def output_branches(self):
     outfile = open(self.fout + ".line.txt", "w")
     for node in self.tree.traverse(strategy="postorder"):
         if not node.is_leaf():
             cn = self.name_coords[node.name]
             childs = node.get_children()
             cl = self.name_coords[childs[0].name]
             cr = self.name_coords[childs[1].name]
             l_mds_dis = self.innernode_dis_mds_matrix[node.name + ":" + childs[0].name]
             l_tree_dis = self.innernode_dis_tree_matrix[node.name + ":" + childs[0].name] 
             r_mds_dis = self.innernode_dis_mds_matrix[node.name + ":" + childs[1].name]
             r_tree_dis = self.innernode_dis_tree_matrix[node.name + ":" + childs[1].name] 
             l_ratio = (l_tree_dis - l_mds_dis)/l_mds_dis
             r_ratio = (r_tree_dis - r_mds_dis)/r_mds_dis
             lstroke = math.erf(l_ratio)
             rstroke = math.erf(r_ratio)
             if lstroke < 0:
                 lstroke = 2.0
             else:
                 lstroke = 4.0*lstroke + 2.0
             
             if rstroke < 0:
                 rstroke = 2.0
             else:
                 rstroke = 4.0*rstroke + 2.0 
             
             outfile.write(repr(cn[0])+","+repr(cn[1])+","+repr(cl[0])+","+repr(cl[1])+","+repr(lstroke)+"\n")
             outfile.write(repr(cn[0])+","+repr(cn[1])+","+repr(cr[0])+","+repr(cr[1])+","+repr(rstroke)+"\n")
     outfile.close()

def response( t, sig, maxamp, maxt, fastconst, slowconst ):
    """
    t=0 is assumed to be max of distribution
    """
    # simplistic
    #farg = t/sig
    #f = 0.5*maxamp*np.exp( -0.5*farg*farg )#/sig/np.sqrt(2*3.14159)
    #if t<0:
    #    s = 0.0
    #else:
    #    s = 0.5*maxamp*np.exp( -t/config.slowconst )

    # slow component shape: expo convolved with gaus
    t_smax = 95.0 # peak of only slow component. numerically solved for det. smearing=3.5*15.625 ns, decay time const= 1500 ns
    t_fmax = 105.0 # numerically solved for det. smearing=3.5*15.625 ns, decay time const= 6 ns
    #dt_smax = -10.0 # expect slow comp peak to be 10 ns earlier than fast component peak
    smax = np.exp( sig*sig/(2*slowconst*slowconst) - t_fmax/slowconst )*(1 - math.erf( (sig*sig - slowconst*t_fmax )/(np.sqrt(2)*sig*slowconst ) ) )
    # normalize max at fast component peak
    As = 0.3*maxamp/smax
    s = As*np.exp( sig*sig/(2*slowconst*slowconst) - t/slowconst )*(1 - math.erf( (sig*sig - slowconst*t )/(np.sqrt(2)*sig*slowconst ) ) )
    #s = np.exp( sig*sig/(2*slowconst*slowconst))*(1-math.erf( (sig*sig)/(np.sqrt(2)*sig*slowconst ) ) )
    #s = maxamp*np.exp( sig*sig/(2*slowconst*slowconst) - t/slowconst )*(1 - math.erf( (sig*sig - slowconst*t )/(np.sqrt(2)*sig*slowconst ) ) )

    # fast component: since time const is smaller than spe response, we model as simple gaussian
    #
    farg = t/sig
    fmax = np.exp( -0.5*farg*farg )
    Af = 0.8*maxamp
    f = Af*np.exp( -0.5*farg*farg )

    #return fastfraction*f + slowfraction*s
    #print t, f, s
    return f+s

def dynamical_friction_sis(x, y, z, vx, vy, vz, M_sat):
    x = x * units.kpc
    y = y * units.kpc
    z = z * units.kpc
    r = np.sqrt(x**2 + y**2 + z**2)
    vx = vx * units.km / units.s
    vy = vy * units.km / units.s
    vz = vz * units.km / units.s
    v = np.sqrt(vx**2 + vy**2 + vz**2)
    a = 10 # concentration parameter
    # Density at distance r and a velocity v
    rho = dens_sis(10, r.value, v.value) # a, r, v
    v = v.to(units.kpc / units.s)
    M_sat = M_sat * units.Msun
    factor = - 4 * np.pi * G**2  
    Coulomb = coulomb_log(r)
    sigma = v / np.sqrt(2)
    X = v / ( np.sqrt(2) * sigma ) #Check this
    F_dfx = factor * M_sat * rho * Coulomb / v**3 * (  erf(X) - 2*X/(np.sqrt(np.pi) * np.exp(-X**2))  ) * vx
    F_dfy = factor * M_sat * rho * Coulomb / v**3 * (  erf(X) - 2*X/(np.sqrt(np.pi) * np.exp(-X**2))  ) * vy
    F_dfz = factor * M_sat * rho * Coulomb / v**3 * (  erf(X) - 2*X/(np.sqrt(np.pi) * np.exp(-X**2))  ) * vz
    F_dfx = F_dfx.to(units.kpc / units.Gyr**2)
    F_dfy = F_dfy.to(units.kpc / units.Gyr**2)
    F_dfz = F_dfz.to(units.kpc / units.Gyr**2)
    return F_dfx.value, F_dfy.value, F_dfz.value

  def __call__(self, val):
    import numpy
    from math import sqrt, log, pi, erf
    sigma_m = val
    #sigma_m = abs(val[0])

    dphi2 = self.dphi / 2
    den =  sqrt(2) * sigma_m / abs(self.zeta)# / 0.670

    a = (self.tau + dphi2) / den
    b = (self.tau - dphi2) / den

    a = numpy.array([erf(a[i]) for i in range(len(a))])
    b = numpy.array([erf(b[i]) for i in range(len(b))])

    r = (a - b) / 2.0#(2 * self.dphi)

#        ind = numpy.where(r > 0)[0]
#        ret = r
#        ret[ind] = numpy.log(r[ind])
    ind = numpy.where(r <= 0)
    if len(ind[0]) > 0:
      r[ind] = 1e-7
    #ret[ind] = 0
    #print ind, r[ind], self.zeta[ind], self.tau[ind]
    return numpy.log(r)

def p_y(n,y,M):
    t_1 = m.exp(-0.5*n*y**2)
    t_2 = m.erf((0.5*n)**(0.5)*(y-0.5*M)) 
    t_3 = m.erf((0.5*n)**(0.5)*(y+0.5*M))
    m_f1 = (n/(8*m.pi))**2
    m_f2 = (n/(2*m.pi))**2/M
    return m_f1*(t_1 - m_f2*(t_2 - t_3))

def black_calc():
    stage4 = fmv * math.exp(-div * exp_term)
    stage5 = (1.0 + math.erf(d1 / math.sqrt(2.0))) / 2.0
    stage6 = exp_price * math.exp(-rate * exp_term)
    stage7 = (1.0 + math.erf(d2 / math.sqrt(2.0))) / 2.0
    c = (stage4*stage5)-(stage6*stage7)
    return c

def calc_charge_loss_fraction_in_line(accelerator, **kwargs):
    """Calculate charge loss in a line

    Keyword arguments:
    twiss_at_entrance -- Twiss parameters at the start of first element
    global_coupling   -- Global coupling
    energy_spread     -- Relative energy spread
    emittance         -- [m·rad]
    delta_rx          -- [m]
    delta_angle       -- [rad]
    hmax              -- [m]
    hmin              -- [m]
    vmax              -- [m]
    vmin              -- [m]
    """
    init_twiss, energy_spread, emittance, hmax, hmin, vmax, vmin = _process_loss_fraction_args(accelerator, **kwargs)
    coupling = kwargs['global_coupling']

    try:
        twiss, m66 = pyaccel.optics.calc_twiss(accelerator, init_twiss = init_twiss, indices ='open')
        betax, etax, betay, etay = twiss.betax, twiss.etax, twiss.betay, twiss.etay
        if math.isnan(betax[-1]):
            loss_fraction = 1.0
            return (loss_fraction, None, None)
    except (numpy.linalg.linalg.LinAlgError, pyaccel.optics.OpticsException, pyaccel.tracking.TrackingException):
        loss_fraction = 1.0
        return (loss_fraction, None, None)

    emitx = emittance * 1 / (1 + coupling)
    emity = emittance * coupling / (1 + coupling)
    sigmax = numpy.sqrt(betax * emitx + (etax * energy_spread)**2)
    sigmay = numpy.sqrt(betay * emity + (etax * energy_spread)**2)
    h_vc = hmax - hmin
    v_vc = vmax - vmin
    co = twiss.co
    rx, ry = co[0,:], co[2,:]
    xlim_inf, xlim_sup = rx - hmin, hmax - rx
    ylim_inf, ylim_sup = ry - vmin, vmax - ry
    xlim_inf[xlim_inf < 0] = 0
    xlim_sup[xlim_sup < 0] = 0
    ylim_inf[ylim_inf < 0] = 0
    ylim_sup[ylim_sup < 0] = 0
    xlim_inf[xlim_inf > h_vc] = 0
    xlim_sup[xlim_sup > h_vc] = 0
    ylim_inf[ylim_inf > v_vc] = 0
    ylim_sup[ylim_sup > v_vc] = 0
    min_xfrac_inf = numpy.amin(xlim_inf/sigmax)
    min_xfrac_sup = numpy.amin(xlim_sup/sigmax)
    min_yfrac_inf = numpy.amin(ylim_inf/sigmay)
    min_yfrac_sup = numpy.amin(ylim_sup/sigmay)
    sqrt2 = math.sqrt(2)

    x_surviving_fraction = 0.5*math.erf(min_xfrac_inf/sqrt2) + \
                           0.5*math.erf(min_xfrac_sup/sqrt2)
    y_surviving_fraction = 0.5*math.erf(min_yfrac_inf/sqrt2) + \
                           0.5*math.erf(min_yfrac_sup/sqrt2)
    surviving_fraction = x_surviving_fraction * y_surviving_fraction
    loss_fraction = 1.0 - surviving_fraction
    return loss_fraction, twiss, m66

 def marg_rating(self):
     if self._marg_rating is None:
         s = self.sigmac2
         sp4 = s+4
         self._marg_rating = 0.5*log(2*pi*s/sp4)*exp(-9/(8*sp4))
         self._marg_rating += log(erf(0.75*sqrt(0.5*s/sp4))
                                  - erf(-0.25*(17*s+80)/sqrt(2*s*sp4)))
     return self._marg_rating

def alpha2(a, N, Nmax, Nmin=1):
    y = sqrt(pi*Nmin*Nmax)/(2.0*a) * exp((a * log2(sqrt(Nmax/Nmin)))**2.0)
    y = y * exp((log(2.0)/(2.0*a))**2.0)
    y = y * erf(a * log2(sqrt(Nmax/Nmin)) - log(2.0)/(2.0*a))
    y += erf(a * log2(sqrt(Nmax/Nmin)) + log(2.0)/(2.0*a))
    y -= N

    return y # find alpha

def test_genz_gaussian_exact():
    u = np.array([1, 21, 2], dtype=float)
    a = np.array([1/10, 1/100, 1/500], dtype=float)
    val = ti.genz_gaussian_exact(u, a)
    exact = 62500*pow(np.pi, 3/2)*math.erf(1/10)*(math.erf(1/250) -
            math.erf(1/500))*(math.erf(21/100) - math.erf(1/5))

    assert np.allclose([val], [exact])

 def transitionProbability(self, T_total, rho0, rhoBar,sigmaRhoSquared,tauEff):
     
     # argument for the Error Function
     x1 =  -(self.rhoStar - rhoBar + (rhoBar - rho0)*np.exp(-T_total/tauEff))/(np.sqrt(sigmaRhoSquared*(1.-np.exp(-2.*T_total/tauEff))))
     # transition probability
     if rho0 == 0.:
         return (1. + math.erf(x1))/2.
     else:
         return (1. - math.erf(x1))/2.

 def _cdf(self, x):
     """
     Calculate cumulative distribution function in a certain point
     """
     if isinstance(x, float):
         return 1.0 / 2.0 * (1 + math.erf(x / np.sqrt(2)))
     else:
         return (
             1.0 / 2.0 *
             (1 + np.array([math.erf(n / np.sqrt(2)) for n in x]))
         )

    def alpha(self, a, Nmax, Nmin=1):

        """Numerically solve for Preston's a. Needed to estimate S using the lognormal"""

        y = sqrt(pi*Nmin*Nmax)/(2.0*a) * exp((a * log2(sqrt(Nmax/Nmin)))**2.0)
        y = y * exp((log(2.0)/(2.0*a))**2.0)
        y = y * erf(a * log2(sqrt(Nmax/Nmin)) - log(2.0)/(2.0*a))
        y += erf(a * log2(sqrt(Nmax/Nmin)) + log(2.0)/(2.0*a))
        y -= self.N

        return y # find alpha

 def _phit(self, z, zstar):
     """
     Model component of hitting target
     :param z: Reading, in cm
     :param zstar: Expected reading, from ray casting
     :return p: Probability of 'hit'
     """
     if z < self._max_z:
         N = 1.0/math.sqrt(2*math.pi*self._sigmahit**2)*math.e**(-0.5*(z-zstar)**2/self._sigmahit**2)
         eta = 0.5*(math.erf((self._max_z-zstar)/(self._sigmahit*math.sqrt(2))) + math.erf(zstar/(math.sqrt(2)*self._sigmahit)))
         return N*eta
     else:
         return 0

def erfContribution(dx, sigma, N):
    	"""
    	Solves for the distance to the closest periodic replica of x1 to x2. Integrates a gaussian of stdev sigma.
    	x1 and x2 are in reduced units, and N is the number of reduced units in the cubic system.
	Quantities in dx must be in [-N, N].
	Uses the same distance computation as GaussianContribution.
    	"""
	half = float(N)/2.0
    	dist_vec = [ abs(abs(abs(x)-half)-half) for x in dx]
	half_box = .5/sigma
    	my_erf = lambda x : - math.erf(-half_box - 2*half_box*x) + math.erf(half_box - 2*half_box*x)
    	contribution = reduce( mul, [ my_erf(x) for x in dist_vec] ) / 8
    	return contribution

  def run(self):
    from dials.algorithms.profile_model.gaussian_rs import \
      PartialityCalculator3D
    from dials.array_family import flex

    calculator = PartialityCalculator3D(
      self.experiment.beam,
      self.experiment.goniometer,
      self.experiment.scan,
      self.sigma_m)

    predicted = flex.reflection_table.from_predictions_multi(self.experiments)
    predicted['bbox'] = predicted.compute_bbox(self.experiments)

    # Remove any touching edges of scan to get only fully recorded
    x0, x1, y0, y1, z0, z1 = predicted['bbox'].parts()
    predicted = predicted.select((z0 > 0) & (z1 < 100))
    assert(len(predicted) > 0)

    # Compute partiality
    partiality = calculator(
      predicted['s1'],
      predicted['xyzcal.px'].parts()[2],
      predicted['bbox'])

    # Should have all fully recorded
    assert(len(partiality) == len(predicted))
    from math import sqrt, erf
    three_sigma = 0.5 * (erf(3.0 / sqrt(2.0)) - erf(-3.0 / sqrt(2.0)))
    assert(partiality.all_gt(three_sigma))

    # Trim bounding boxes
    x0, x1, y0, y1, z0, z1 = predicted['bbox'].parts()
    z0 = z0 + 1
    z1 = z1 - 1
    predicted['bbox'] = flex.int6(x0, x1, y0, y1, z0, z1)
    predicted = predicted.select(z1 > z0)
    assert(len(predicted) > 0)

    # Compute partiality
    partiality = calculator(
      predicted['s1'],
      predicted['xyzcal.px'].parts()[2],
      predicted['bbox'])

    # Should have all partials
    assert(len(partiality) == len(predicted))
    assert(partiality.all_le(1.0) and partiality.all_gt(0))

    print 'OK'

def pvalue(x):
    '''Cumulative distribution function for the standard normal distribution
    python
    import math
    x=2.33
    x=1.96
    1-(1.0 + math.erf(x / math.sqrt(2.0))) / 2.0
    x=-2.33
    (1.0 + math.erf(x / math.sqrt(2.0))) / 2.0
    '''
    if x < 0:
        return (1.0 + math.erf(x / math.sqrt(2.0))) / 2.0
    else:
        return 1-(1.0 + math.erf(x / math.sqrt(2.0))) / 2.0

def genDocument(cur):

	mu = random.randrange(0, g_wordsNum)
	sigma = random.randrange(1, 10*precision)/precision
	docSize = random.randrange(10, 100)
	cur.execute("INSERT INTO Documents (docName, docSize, mu, sigma)\
		VALUES(\"\",\"" + str(docSize)+ "\",\"" + str(mu)+ "\", \"" + str(sigma) + "\" )")
	docID = cur.lastrowid
	yFrom = (1 + math.erf((0 - mu)/ (sigma * math.sqrt(2)))) / 2
	yTo   = (1 + math.erf((g_wordsNum - mu)/ (sigma * math.sqrt(2)))) / 2

	for i in range(1, docSize):
		addInvertedID(docID, randomWord(mu, sigma, yFrom, yTo), cur)
	
	return docID

    def sort_key(self):
        """Returns 1.0 for a ratio of 1, falling to 0.0 for extremely small
        or large values."""

        score = math.log10(self.score)
        cdf = (1.0 + math.erf(score / math.sqrt(2.0))) / 2.0
        return 1 - 2*math.fabs(0.5 - cdf)

    def calc_stats(self, count, target, total):
        num = max(len(target), 1)
        r = {}
        for t in ('pgs', 'objects', 'bytes'):
            if total[t] == 0:
                r[t] = {
                    'avg': 0,
                    'stddev': 0,
                    'sum_weight': 0,
                    'score': 0,
                }
                continue

            avg = float(total[t]) / float(num)
            dev = 0.0

            # score is a measure of how uneven the data distribution is.
            # score lies between [0, 1), 0 means perfect distribution.
            score = 0.0
            sum_weight = 0.0

            for k, v in six.iteritems(count[t]):
                # adjust/normalize by weight
                if target[k]:
                    adjusted = float(v) / target[k] / float(num)
                else:
                    adjusted = 0.0

                # Overweighted devices and their weights are factors to calculate reweight_urgency.
                # One 10% underfilled device with 5 2% overfilled devices, is arguably a better
                # situation than one 10% overfilled with 5 2% underfilled devices
                if adjusted > avg:
                    '''
                    F(x) = 2*phi(x) - 1, where phi(x) = cdf of standard normal distribution
                    x = (adjusted - avg)/avg.
                    Since, we're considering only over-weighted devices, x >= 0, and so phi(x) lies in [0.5, 1).
                    To bring range of F(x) in range [0, 1), we need to make the above modification.

                    In general, we need to use a function F(x), where x = (adjusted - avg)/avg
                    1. which is bounded between 0 and 1, so that ultimately reweight_urgency will also be bounded.
                    2. A larger value of x, should imply more urgency to reweight.
                    3. Also, the difference between F(x) when x is large, should be minimal.
                    4. The value of F(x) should get close to 1 (highest urgency to reweight) with steeply.

                    Could have used F(x) = (1 - e^(-x)). But that had slower convergence to 1, compared to the one currently in use.

                    cdf of standard normal distribution: https://stackoverflow.com/a/29273201
                    '''
                    score += target[k] * (math.erf(((adjusted - avg)/avg) / math.sqrt(2.0)))
                    sum_weight += target[k]
                dev += (avg - adjusted) * (avg - adjusted)
            stddev = math.sqrt(dev / float(max(num - 1, 1)))
            score = score / max(sum_weight, 1)
            r[t] = {
                'avg': avg,
                'stddev': stddev,
                'sum_weight': sum_weight,
                'score': score,
            }
        return r