Python special.erfc方法代码示例

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import erfc [as 别名]
def test_erf_complex():
    # need to increase mpmath precision for this test
    old_dps, old_prec = mpmath.mp.dps, mpmath.mp.prec
    try:
        mpmath.mp.dps = 70
        x1, y1 = np.meshgrid(np.linspace(-10, 1, 31), np.linspace(-10, 1, 11))
        x2, y2 = np.meshgrid(np.logspace(-80, .8, 31), np.logspace(-80, .8, 11))
        points = np.r_[x1.ravel(),x2.ravel()] + 1j*np.r_[y1.ravel(),y2.ravel()]

        assert_func_equal(sc.erf, lambda x: complex(mpmath.erf(x)), points,
                          vectorized=False, rtol=1e-13)
        assert_func_equal(sc.erfc, lambda x: complex(mpmath.erfc(x)), points,
                          vectorized=False, rtol=1e-13)
    finally:
        mpmath.mp.dps, mpmath.mp.prec = old_dps, old_prec


#------------------------------------------------------------------------------
# lpmv
#------------------------------------------------------------------------------ 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import erfc [as 别名]
def test_erf_complex():
    # need to increase mpmath precision for this test
    old_dps, old_prec = mpmath.mp.dps, mpmath.mp.prec
    try:
        mpmath.mp.dps = 70
        x1, y1 = np.meshgrid(np.linspace(-10, 1, 31), np.linspace(-10, 1, 11))
        x2, y2 = np.meshgrid(np.logspace(-80, .8, 31), np.logspace(-80, .8, 11))
        points = np.r_[x1.ravel(),x2.ravel()] + 1j*np.r_[y1.ravel(), y2.ravel()]

        assert_func_equal(sc.erf, lambda x: complex(mpmath.erf(x)), points,
                          vectorized=False, rtol=1e-13)
        assert_func_equal(sc.erfc, lambda x: complex(mpmath.erfc(x)), points,
                          vectorized=False, rtol=1e-13)
    finally:
        mpmath.mp.dps, mpmath.mp.prec = old_dps, old_prec


# ------------------------------------------------------------------------------
# lpmv
# ------------------------------------------------------------------------------ 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import erfc [as 别名]
def lempelzivcompressiontest1(binin):
    ''' The focus of this test is the number of cumulatively distinct patterns (words) in the sequence. The purpose of the test is to determine how far the tested sequence can be compressed. The sequence is considered to be non-random if it can be significantly compressed. A random sequence will have a characteristic number of distinct patterns.'''
    i = 1
    j = 0
    n = len(binin)
    mu = 69586.25
    sigma = 70.448718
    words = []
    while (i+j) <= n:
        tmp = binin[i:i+j:]
        if words.count(tmp) > 0:
            j += 1
        else:
            words.append(tmp)
            i += j+1
            j = 0
    wobs = len(words)
    pval = 0.5*spc.erfc((mu-wobs)/np.sqrt(2.0*sigma))
    return pval

# test 2.11 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import erfc [as 别名]
def get_selu_parameters(fixed_point_mean=0.0, fixed_point_var=1.0):
    """ Finding the parameters of the SELU activation function. The function returns alpha and lambda for the desired
    fixed point. """
    aa = Symbol('aa')
    ll = Symbol('ll')
    nu = fixed_point_mean
    tau = fixed_point_var
    mean = 0.5 * ll * (nu + np.exp(-nu ** 2 / (2 * tau)) * np.sqrt(2 / np.pi) * np.sqrt(tau) +
                       nu * erf(nu / (np.sqrt(2 * tau))) - aa * erfc(nu / (np.sqrt(2 * tau))) +
                       np.exp(nu + tau / 2) * aa * erfc((nu + tau) / (np.sqrt(2 * tau))))
    var = 0.5 * ll ** 2 * (np.exp(-nu ** 2 / (2 * tau)) * np.sqrt(2 / np.pi * tau) * nu + (nu ** 2 + tau) *
                           (1 + erf(nu / (np.sqrt(2 * tau)))) + aa ** 2 * erfc(nu / (np.sqrt(2 * tau))) -
                           aa ** 2 * 2 * np.exp(nu + tau / 2) * erfc((nu + tau) / (np.sqrt(2 * tau))) +
                           aa ** 2 * np.exp(2 * (nu + tau)) * erfc((nu + 2 * tau) / (np.sqrt(2 * tau)))) - mean ** 2
    eq1 = mean - nu
    eq2 = var - tau
    res = nsolve((eq2, eq1), (aa, ll), (1.67, 1.05))
    return float(res[0]), float(res[1]) 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import erfc [as 别名]
def _log_erfc(x):
  """Compute log(erfc(x)) with high accuracy for large x."""
  try:
    return math.log(2) + special.log_ndtr(-x * 2**.5)
  except NameError:
    # If log_ndtr is not available, approximate as follows:
    r = special.erfc(x)
    if r == 0.0:
      # Using the Laurent series at infinity for the tail of the erfc function:
      #     erfc(x) ~ exp(-x^2-.5/x^2+.625/x^4)/(x*pi^.5)
      # To verify in Mathematica:
      #     Series[Log[Erfc[x]] + Log[x] + Log[Pi]/2 + x^2, {x, Infinity, 6}]
      return (-math.log(math.pi) / 2 - math.log(x) - x**2 - .5 * x**-2 +
              .625 * x**-4 - 37. / 24. * x**-6 + 353. / 64. * x**-8)
    else:
      return math.log(r) 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import erfc [as 别名]
def get_quantiles(acquisition_par, fmin, m, s):
    '''
    Quantiles of the Gaussian distribution useful to determine the acquisition function values
    :param acquisition_par: parameter of the acquisition function
    :param fmin: current minimum.
    :param m: vector of means.
    :param s: vector of standard deviations.
    '''
    if isinstance(s, np.ndarray):
        s[s<1e-10] = 1e-10
    elif s< 1e-10:
        s = 1e-10
    u = (fmin - m - acquisition_par)/s

    phi = np.exp(-0.5 * u**2) / np.sqrt(2*np.pi)
    # vectorized version of erfc to not depend on scipy
    Phi = 0.5 * erfc(-u / np.sqrt(2))
    return (phi, Phi, u) 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import erfc [as 别名]
def get_ewald_params(cell, precision=INTEGRAL_PRECISION, mesh=None):
    r'''Choose a reasonable value of Ewald 'eta' and 'cut' parameters.
    eta^2 is the exponent coefficient of the model Gaussian charge for nucleus
    at R:  \frac{eta^3}{pi^1.5} e^{-eta^2 (r-R)^2}

    Choice is based on largest G vector and desired relative precision.

    The relative error in the G-space sum is given by

        precision ~ 4\pi Gmax^2 e^{(-Gmax^2)/(4 \eta^2)}

    which determines eta. Then, real-space cutoff is determined by (exp.
    factors only)

        precision ~ erfc(eta*rcut) / rcut ~ e^{(-eta**2 rcut*2)}

    Returns:
        ew_eta, ew_cut : float
            The Ewald 'eta' and 'cut' parameters.
    '''
    if cell.natm == 0:
        return 0, 0
    elif (cell.dimension < 2 or
          (cell.dimension == 2 and cell.low_dim_ft_type == 'inf_vacuum')):
# Non-uniform PW grids are used for low-dimensional ewald summation.  The cutoff
# estimation for long range part based on exp(G^2/(4*eta^2)) does not work for
# non-uniform grids.  Smooth model density is preferred.
        ew_cut = cell.rcut
        ew_eta = np.sqrt(max(np.log(4*np.pi*ew_cut**2/precision)/ew_cut**2, .1))
    else:
        if mesh is None:
            mesh = cell.mesh
        mesh = _cut_mesh_for_ewald(cell, mesh)
        Gmax = min(np.asarray(mesh)//2 * lib.norm(cell.reciprocal_vectors(), axis=1))
        log_precision = np.log(precision/(4*np.pi*(Gmax+1e-100)**2))
        ew_eta = np.sqrt(-Gmax**2/(4*log_precision)) + 1e-100
        ew_cut = _estimate_rcut(ew_eta**2, 0, 1., precision)
    return ew_eta, ew_cut 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import erfc [as 别名]
def ewald_ion(self):
        r"""
        Compute ion contribution to Ewald sums.  Since the ions don't move in our calculations, the ion-ion term only needs to be computed once.

        Note: We ignore the constant term :math:`\frac{1}{2} \sum_{I} Z_I^2 C_{\rm self\ image}` in the real-space ion-ion sum corresponding to the interaction of an ion with its own image in other cells.

        The real-space part:

        .. math:: E_{\rm real\ space}^{\text{ion-ion}} = \sum_{\vec{n}} \sum_{I<J}^{N_{ion}} Z_I Z_J \frac{{\rm erfc}(\alpha |\vec{x}_{IJ}+\vec{n}|)}{|\vec{x}_{IJ}+\vec{n}|} 

        The reciprocal-space part:

        .. math:: E_{\rm reciprocal\ space}^{\text{ion-ion}} = \sum_{\vec{G} > 0 } W_G \left| \sum_{I=1}^{N_{ion}} Z_I e^{-i\vec{G}\cdot\vec{x}_I} \right|^2

        Returns:
            ion_ion: float, ion-ion component of Ewald sum
        """
        # Real space part
        if len(self.atom_charges) == 1:
            ion_ion_real = 0
        else:
            dist = pyqmc.distance.MinimalImageDistance(self.latvec)
            ion_distances, ion_inds = dist.dist_matrix(self.atom_coords[np.newaxis])
            rvec = ion_distances[:, :, np.newaxis, :] + self.lattice_displacements
            r = np.linalg.norm(rvec, axis=-1)
            charge_ij = np.prod(self.atom_charges[np.asarray(ion_inds)], axis=1)
            ion_ion_real = np.einsum("j,ijk->", charge_ij, erfc(self.alpha * r) / r)

        # Reciprocal space part
        GdotR = np.dot(self.gpoints, self.atom_coords.T)
        self.ion_exp = np.dot(np.exp(1j * GdotR), self.atom_charges)
        ion_ion_rec = np.dot(self.gweight, np.abs(self.ion_exp) ** 2)

        ion_ion = ion_ion_real + ion_ion_rec
        return ion_ion 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import erfc [as 别名]
def _real_cij(self, dists):
        r = np.zeros(dists.shape[:-1])
        cij = np.zeros(r.shape)
        for ld in self.lattice_displacements:
            r[:] = np.linalg.norm(dists + ld, axis=-1)
            cij += erfc(self.alpha * r) / r
        return cij 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import erfc [as 别名]
def energy_with_test_pos(self, configs, epos):
        """
        Compute Coulomb energy of an additional test electron with a set of configs

        Inputs:
            configs: pyqmc PeriodicConfigs object of shape (nconf, nelec, ndim)
            epos: pyqmc PeriodicConfigs object of shape (nconf, ndim)
        Returns: 
            Vtest: (nconf, nelec+1) array. The first nelec columns are Coulomb energies between the test electron and each electron; the last column is the contribution from all the ions.
        """
        nconf, nelec, ndim = configs.configs.shape
        Vtest = np.zeros((nconf, nelec + 1)) + self.ijconst
        Vtest[:, -1] = self.e_single_test

        # Real space electron-ion part
        # ei_distances shape (conf, atom, dim)
        ei_distances = configs.dist.dist_i(self.atom_coords, epos.configs)
        rvec = ei_distances[:, :, np.newaxis, :] + self.lattice_displacements
        r = np.linalg.norm(rvec, axis=-1)
        Vtest[:, -1] += np.einsum(
            "k,jkl->j", -self.atom_charges, erfc(self.alpha * r) / r
        )

        # Real space electron-electron part
        ee_distances = configs.dist.dist_i(configs.configs, epos.configs)
        rvec = ee_distances[:, :, np.newaxis, :] + self.lattice_displacements
        r = np.linalg.norm(rvec, axis=-1)
        Vtest[:, :-1] += np.sum(erfc(self.alpha * r) / r, axis=-1)

        # Reciprocal space electron-electron part
        e_expGdotR = np.exp(1j * np.dot(configs.configs, self.gpoints.T))
        test_exp = np.exp(1j * np.dot(epos.configs, self.gpoints.T))
        ee_recip_separated = np.dot(np.real(test_exp.conj() * e_expGdotR), self.gweight)
        Vtest[:, :-1] += 2 * ee_recip_separated

        # Reciprocal space electrin-ion part
        coscos_sinsin = np.real(-self.ion_exp.conj() * test_exp)
        ei_recip_separated = np.dot(coscos_sinsin + 0.5, self.gweight)
        Vtest[:, -1] += 2 * ei_recip_separated

        return Vtest 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import erfc [as 别名]
def log_relative_gauss(z):
    if z < -6:
        return 1, -1.0e12, -1
    if z > 6:
        return 0, 0, 1
    else:
        logphi = -0.5 * (z * z + l2p)
        logPhi = np.log(.5 * special.erfc(-z / sq2))
        e = np.exp(logphi - logPhi)
        return e, logPhi, 0 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import erfc [as 别名]
def _pdf(self, x, K):
        # exponnorm.pdf(x, K) =
        #     1/(2*K) exp(1/(2 * K**2)) exp(-x / K) * erfc-(x - 1/K) / sqrt(2))
        invK = 1.0 / K
        exparg = 0.5 * invK**2 - invK * x
        # Avoid overflows; setting np.exp(exparg) to the max float works
        #  all right here
        expval = _lazywhere(exparg < _LOGXMAX, (exparg,), np.exp, _XMAX)
        return 0.5 * invK * expval * sc.erfc(-(x - invK) / np.sqrt(2)) 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import erfc [as 别名]
def _logpdf(self, x, K):
        invK = 1.0 / K
        exparg = 0.5 * invK**2 - invK * x
        return exparg + np.log(0.5 * invK * sc.erfc(-(x - invK) / np.sqrt(2))) 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import erfc [as 别名]
def _cdf(self, x):
        # Equivalent to 2*norm.sf(np.sqrt(1/x))
        return sc.erfc(np.sqrt(0.5 / x)) 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import erfc [as 别名]
def test_erfc(self):
        assert_mpmath_equal(sc.erfc,
                            _exception_to_nan(lambda z: mpmath.erfc(z)),
                            [Arg()])