Python scipy.special.erfc() Examples

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) 

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) 

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
# ------------------------------------------------------------------------------ 

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) 

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]) 

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
#------------------------------------------------------------------------------ 

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 

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 

def test_erfcx_consistent(self):
        self._check_variant_func(
            cephes.erfcx,
            lambda z: np.exp(z*z) * cephes.erfc(z),
            rtol=1e-12
            ) 

def test_erfc_nan_inf(self):
        vals = [np.nan, -np.inf, np.inf]
        expected = [np.nan, 2, 0]
        assert_allclose(special.erfc(vals), expected, rtol=1e-15) 

def test_infs():
    h = 1
    res = 0.5*sc.erfc(h/np.sqrt(2))
    assert_allclose(sc.owens_t(h, np.inf), res, rtol=5e-14)
    assert_allclose(sc.owens_t(h, -np.inf), -res, rtol=5e-14)

    assert_equal(sc.owens_t(np.inf, 1), 0)
    assert_equal(sc.owens_t(-np.inf, 1), 0)

    assert_equal(sc.owens_t(np.inf, np.inf), 0)
    assert_equal(sc.owens_t(-np.inf, np.inf), 0)
    assert_equal(sc.owens_t(np.inf, -np.inf), -0.0)
    assert_equal(sc.owens_t(-np.inf, -np.inf), -0.0) 

def test_erfc_complex(self):
        assert_mpmath_equal(sc.erfc,
                            exception_to_nan(lambda z: mpmath.erfc(z)),
                            [ComplexArg()], n=200) 

def test_ndtr_complex(self):
        assert_mpmath_equal(sc.ndtr,
                            lambda z: mpmath.erfc(-z/np.sqrt(2.))/2.,
                            [ComplexArg(a=complex(-10000, -10000), b=complex(10000, 10000))], n=400) 

def test_log_ndtr_complex(self):
        assert_mpmath_equal(sc.log_ndtr,
                            exception_to_nan(lambda z: mpmath.log(mpmath.erfc(-z/np.sqrt(2.))/2.)),
                            [ComplexArg(a=complex(-10000, -100),
                                        b=complex(10000, 100))], n=200, dps=300) 

def monobitfrequencytest(binin):
    ''' The focus of the test is the proportion of zeroes and ones for the entire sequence. The purpose of this test is to determine whether that number of ones and zeros in a sequence are approximately the same as would be expected for a truly random sequence. The test assesses the closeness of the fraction of ones to 1/2, that is, the number of ones and zeroes in a sequence should be about the same.'''

    ss = [int(el) for el in binin]
    sc = list(map(sumi, ss))
    sn = reduce(su, sc)
    sobs = np.abs(sn) / np.sqrt(len(binin))
    pval = spc.erfc(sobs / np.sqrt(2))
    return pval 

def runstest(binin):
    ''' The focus of this test is the total number of zero and one runs in the entire sequence, where a run is an uninterrupted sequence of identical bits. A run of length k means that a run consists of exactly k identical bits and is bounded before and after with a bit of the opposite value. The purpose of the runs test is to determine whether the number of runs of ones and zeros of various lengths is as expected for a random sequence. In particular, this test determines whether the oscillation between such substrings is too fast or too slow.'''
    ss = [int(el) for el in binin]
    n = len(binin)
    pi = 1.0 * reduce(su, ss) / n
    vobs = len(binin.replace('0', ' ').split()) + \
        len(binin.replace('1', ' ').split())
    pval = spc.erfc(abs(vobs-2*n*pi*(1-pi)) /
                    (2 * pi * (1 - pi) * np.sqrt(2*n)))
    return pval 

def spectraltest(binin):
    '''The focus of this test is the peak heights in the discrete Fast Fourier Transform. The purpose of this test is to detect periodic features (i.e., repetitive patterns that are near each other) in the tested sequence that would indicate a deviation from the assumption of randomness. '''

    n = len(binin)
    ss = [int(el) for el in binin]
    sc = list(map(sumi, ss))
    ft = sff.fft(sc)
    af = abs(ft)[1:floor(n/2)+1:]
    t = np.sqrt(np.log(1/0.05)*n)
    n0 = 0.95*n/2
    n1 = len(np.where(af < t)[0])
    d = (n1 - n0)/np.sqrt(n*0.95*0.05/4)
    pval = spc.erfc(abs(d)/np.sqrt(2))
    return pval 

def maurersuniversalstatistictest(binin, l=6, q=640):
    ''' The focus of this test is the number of bits between matching patterns. The purpose of the test is to detect whether or not the sequence can be significantly compressed without loss of information. An overly compressible sequence is considered to be non-random.'''
    ru = [
        [0.7326495, 0.690],
        [1.5374383, 1.338],
        [2.4016068, 1.901],
        [3.3112247, 2.358],
        [4.2534266, 2.705],
        [5.2177052, 2.954],
        [6.1962507, 3.125],
        [7.1836656, 3.238],
        [8.1764248, 3.311],
        [9.1723243, 3.356],
        [10.170032, 3.384],
        [11.168765, 3.401],
        [12.168070, 3.410],
        [13.167693, 3.416],
        [14.167488, 3.419],
        [15.167379, 3.421],
    ]
    blocks = [int(li, 2) + 1 for li in stringpart(binin, l)]
    k = len(blocks) - q
    states = [0 for x in range(2**l)]
    for x in range(q):
        states[blocks[x]-1] = x+1
    sumi = 0.0
    for x in range(q, len(blocks)):
        sumi += np.log2((x+1)-states[blocks[x]-1])
        states[blocks[x]-1] = x+1
    fn = sumi / k
    c = 0.7-(0.8/l)+(4+(32.0/l))*((k**(-3.0/l))/15)
    sigma = c*np.sqrt((ru[l-1][1])/k)
    pval = spc.erfc(abs(fn-ru[l-1][0]) / (np.sqrt(2)*sigma))
    return pval 

def Q_fctn(x):
    """
    Gaussian Q-function
    """
    return 1./2*erfc(x/np.sqrt(2.)) 

def log_relative_gauss(z):
    """
    log_relative_gauss
    """
    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 

def _pdf(self, x, K):
        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)) 

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))) 

def _cdf(self, x):
        # Equivalent to 2*norm.sf(np.sqrt(1/x))
        return sc.erfc(np.sqrt(0.5 / x)) 

def _log_erfc(x):
  # Can be replaced with a single call to log_ntdr if available:
  # return np.log(2.) + special.log_ntdr(-x * 2**.5)
  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) 

def test_erfc(self):
        assert_mpmath_equal(sc.erfc,
                            _exception_to_nan(lambda z: mpmath.erfc(z)),
                            [Arg()]) 

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 

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 

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 

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 

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))