Python math.erf函数代码示例

 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