Python scipy.special() Examples

def fit_transform( self, X ):

        i = np.argsort( X, axis = 0 )
        j = np.argsort( i, axis = 0 )

        assert ( j.min() == 0 ).all()
        assert ( j.max() == len( j ) - 1 ).all()

        j_range = len( j ) - 1
        self.divider = j_range / self.range

        transformed = j / self.divider
        transformed = transformed - self.upper
        transformed = scipy.special.erfinv( transformed )
        ############
        # transformed = transformed - np.mean(transformed)

        return transformed 

def alpha_abs(self, x, y, n_sersic, r_eff, k_eff, center_x=0, center_y=0):
        """

        :param x:
        :param y:
        :param n_sersic:
        :param r_eff:
        :param k_eff:
        :param center_x:
        :param center_y:
        :return:
        """
        n = n_sersic
        x_red = self._x_reduced(x, y, n_sersic, r_eff, center_x, center_y)
        b = self.b_n(n_sersic)
        a_eff = self._alpha_eff(r_eff, n_sersic, k_eff)
        alpha = 2. * a_eff * x_red ** (-n) * (special.gammainc(2 * n, b * x_red))
        return alpha 

def frequency(self, w, s=1.0):
        """Frequency representation of derivative of Gaussian.

        Parameters
        ----------
        w : float
            Angular frequency. If `s` is not specified, i.e. set to 1,
            this can be used as the non-dimensional angular
            frequency w * s.
        s : float
            Scaling factor. Default is 1.

        Returns
        -------
        out : complex
            Value of the derivative of Gaussian wavelet at the
            given time
        """
        m = self.m
        x = s * w
        gamma = scipy.special.gamma
        const = -1j ** m / gamma(m + 0.5) ** .5
        function = x ** m * np.exp(-x ** 2 / 2)
        return const * function 

def test_UnconvertibleValueErrorIsUncatchable(self):
        import scipy.special
        def f():
            try:
                return scipy.special.airy(0)
            except:
                return 0

        try:
            self.evaluateWithExecutor(f)
            self.assertTrue(False)
        except pyfora.ComputationError as e:
            self.assertIsInstance(
                e.remoteException,
                Exceptions.UnconvertibleValueError
                ) 

def test_UnconvertibleValueErrorIsUncatchable(self):
        import scipy.special
        def f():
            try:
                return scipy.special.airy(0)
            except:
                return 0

        try:
            self.evaluateWithExecutor(f)
            self.assertTrue(False)
        except pyfora.ComputationError as e:
            self.assertIsInstance(
                e.remoteException,
                Exceptions.UnconvertibleValueError
                ) 

def test_typeWeCantTranslateYet_raise_4(self):
        import scipy
        def f():
            x = scipy.special
            return 0

        try:
            self.evaluateWithExecutor(f)
            self.assertTrue(False)
        except pyfora.ComputationError as e:
            self.assertTrue(
                str(e).startswith("Pyfora didn't know how to convert scipy.special")
                )
            self.assertIsInstance(
                e.remoteException,
                Exceptions.UnconvertibleValueError
                ) 

def frequency(self, w, s=1.0):
        """Frequency representation of derivative of Gaussian.

        Parameters
        ----------
        w : float
            Angular frequency. If `s` is not specified, i.e. set to 1,
            this can be used as the non-dimensional angular
            frequency w * s.
        s : float
            Scaling factor. Default is 1.

        Returns
        -------
        out : complex
            Value of the derivative of Gaussian wavelet at the
            given time
        """
        m = self.m
        x = s * w
        gamma = scipy.special.gamma
        const = -1j ** m / gamma(m + 0.5) ** .5
        function = x ** m * np.exp(-x ** 2 / 2)
        return const * function 

def frequency(self, w, s=1.0):
        """Frequency representation of derivative of Gaussian.

        Parameters
        ----------
        w : float
            Angular frequency. If `s` is not specified, i.e. set to 1,
            this can be used as the non-dimensional angular
            frequency w * s.
        s : float
            Scaling factor. Default is 1.

        Returns
        -------
        out : complex
            Value of the derivative of Gaussian wavelet at the
            given time
        """
        m = self.m
        x = s * w
        gamma = scipy.special.gamma
        const = -1j ** m / gamma(m + 0.5) ** .5
        function = x ** m * np.exp(-x ** 2 / 2)
        return const * function 

def frequency(self, w, s=1.0):
        """Frequency representation of derivative of Gaussian.

        Parameters
        ----------
        w : float
            Angular frequency. If `s` is not specified, i.e. set to 1,
            this can be used as the non-dimensional angular
            frequency w * s.
        s : float
            Scaling factor. Default is 1.

        Returns
        -------
        out : complex
            Value of the derivative of Gaussian wavelet at the
            given time
        """
        m = self.m
        x = s * w
        gamma = scipy.special.gamma
        const = -1j ** m / gamma(m + 0.5) ** .5
        function = x ** m * np.exp(-x ** 2 / 2)
        return const * function 

def __call__(self, n, k):
        if not isinstance(n, int):
            n = int(n)
        if not isinstance(k, int):
            k = int(k)

        res = 1.0
        if n < 0 or k < 0:
            return 0.0

        if k == 0:
            return res

        if (n - k) < k:
            return scipy.special.comb(n, n - k)

        for ix in xrange(k):
            res = (res * (n - ix)) / (k - ix)

        return res 

def test_function(name, specialization):
    if (name, specialization) in SKIP_LIST:
        pytest.xfail()

    scipy_func = getattr(sc, name)

    @numba.njit
    def numba_func(*args):
        return scipy_func(*args)

    args = itertools.product(*(
        NUMBA_TYPES_TO_TEST_POINTS[numba_type] for numba_type in specialization
    ))
    with warnings.catch_warnings():
        # Ignore warnings about unsafe casts generated by SciPy.
        warnings.filterwarnings(
            action='ignore',
            message='floating point number truncated to an integer',
            category=RuntimeWarning,
        )
        compare_functions(args, scipy_func, numba_func) 

def quadratic(x):
    """A quadratic B-spline.

    This is a special case of `bspline`, and equivalent to ``bspline(x, 2)``.
    """
    ax = abs(asarray(x))
    res = zeros_like(ax)
    cond1 = less(ax, 0.5)
    if cond1.any():
        ax1 = ax[cond1]
        res[cond1] = 0.75 - ax1 ** 2
    cond2 = ~cond1 & less(ax, 1.5)
    if cond2.any():
        ax2 = ax[cond2]
        res[cond2] = (ax2 - 1.5) ** 2 / 2.0
    return res 

def cubic(x):
    """A cubic B-spline.

    This is a special case of `bspline`, and equivalent to ``bspline(x, 3)``.
    """
    ax = abs(asarray(x))
    res = zeros_like(ax)
    cond1 = less(ax, 1)
    if cond1.any():
        ax1 = ax[cond1]
        res[cond1] = 2.0 / 3 - 1.0 / 2 * ax1 ** 2 * (2 - ax1)
    cond2 = ~cond1 & less(ax, 2)
    if cond2.any():
        ax2 = ax[cond2]
        res[cond2] = 1.0 / 6 * (2 - ax2) ** 3
    return res 

def frequency(self, w, s=1.0):
        """Frequency representation of derivative of Gaussian.

        Parameters
        ----------
        w : float
            Angular frequency. If `s` is not specified, i.e. set to 1,
            this can be used as the non-dimensional angular
            frequency w * s.
        s : float
            Scaling factor. Default is 1.

        Returns
        -------
        out : complex
            Value of the derivative of Gaussian wavelet at the
            given time
        """
        m = self.m
        x = s * w
        gamma = scipy.special.gamma
        const = -1j ** m / gamma(m + 0.5) ** .5
        function = x ** m * np.exp(-x ** 2 / 2)
        return const * function 

def df_slice(df, elements, levels, drop_level=True, reset_index=False, return_none=False):
    if df is None:
        return None
    elements, levels = ensure_iterable(elements), ensure_iterable(levels)
    if not len(levels):
        return None
    if len(elements) != len(levels) and len(levels) > 1:
        raise ValueError('Number of elements ' + str(len(elements)) + ' must match the number of levels ' + str(len(levels)))

    # special case where we use a different method to handle multiple elements
    if len(levels) == 1 and len(elements) > 1:
        df =  df.reset_index().loc[df.reset_index()[levels[0]].isin(elements)].set_index(df.index.names)
    else:
        # remove levels if they are not in the df
        elements, levels = zip(*[(e, l) for e, l in zip(elements, levels) if l in df.index.names])
        result = df.xs(elements, level=levels, drop_level=drop_level)
        df = result.reset_index().set_index(result.index.names) if reset_index else result
    if not len(df) and return_none:
        return None
    else:
        return df 

def _convert(self, string):
        """
        Converts string in query to either a Parameter instance, a Numpy function, a scipy.special function or
        a float number.

        Args:
            string(str): string to convert

        Returns:
            Parameter instance, function or float
        """
        if string in self.names:
            param = [p for p in self.parameters if p.name == string][0]
            return param.copy()
        elif isinstance(string, str) and (string in dir(np)) and callable(getattr(np, string)):
            return getattr(np, string)
        elif isinstance(string, str) and (string in dir(sp)) and callable(getattr(sp, string)):
            return getattr(sp, string)
        else:
            return float(string) 

def rel_entr(self, other):
        try:
            naked_other = naked(other)
        except TypeError:  # pragma: no cover
            return NotImplemented

        xp = get_array_module(self.spmatrix)

        if xp is np:
            from scipy.special import rel_entr
        else:  # pragma: no cover
            from cupyx.scipy.special import rel_entr

        if get_array_module(naked_other).isscalar(naked_other):  # pragma: no cover
            return call_sparse_binary_scalar(rel_entr, self, naked_other)
        else:
            if issparse(naked_other):  # pragma: no cover
                naked_other = other.toarray()
            x = get_sparse_module(self.spmatrix).csr_matrix(
                rel_entr(self.toarray(), naked_other))
        if issparse(x):
            return SparseNDArray(x, shape=self.shape)
        return get_array_module(x).asarray(x) 

def detectThresholdToNoiseSignalToNoisepD(SignalToNoise, pD):
    """ Solve for the threshold to noise ratio, given the signal to noise ratio and
    probability of detection.

    References:

    "Electro-optics handbook," Tech. Rep. EOH-11, RCA, 1974. RCA Technical Series Publication.

    R. D. Hippenstiel, Detection Theory: Applications and Digital Signal Pro-cessing, CRC Press, 2002

    Args:
        | SignalToNoise (float): the signal to noise ratio [-]
        | pD (float): the probability of detection [-]

    Returns:
        | range (float): signal to noise ratio

    Raises:
        | No exception is raised.
    """

    import scipy.special

    ThresholdToNoise = SignalToNoise - np.sqrt(2) * scipy.special.erfinv(2 * pD -1)

    return ThresholdToNoise


##############################################################################
## 

def mean_weibul_factor(beta):
    """ beta is shape parameter of weibul
    http://reliawiki.org/index.php/The_Weibull_Distribution
    """
    return scipy.special.gamma(1 + 1. / beta) 

def std_weibul_factor(beta):
    """ beta is shape parameter of weibul
    http://reliawiki.org/index.php/The_Weibull_Distribution
    """
    return ((scipy.special.gamma(1 + 2. / beta)) - (scipy.special.gamma(1 + 1. / beta) ** 2)) ** .5 

def CDF(x,mu,sig, amp=1, lg=0, **kargs):
    if 'Amp' in kargs:
        from warnings import warn
        warn("Parameter Amp is deprecated. Please use amp in order to set the amplitude!")
        amp = kargs.pop('Amp')
    from scipy.special import erf
    g = sig*np.sqrt(2*np.log(2))
    return amp*lg*(.5+np.arctan2(x-mu,g)/np.pi)+(1-lg)*amp*.5*(1+erf((x-mu)/(sig*np.sqrt(2)))) 

def detectSignalToNoiseThresholdToNoisePd(ThresholdToNoise, pD):
    """ Solve for the signal to noise ratio, given the threshold to noise ratio and
    probability of detection.

    Using the theory of matched filter design, calculate the
    signal to noise ratio, to achieve a required probability of detection.

    References:

    "Electro-optics handbook," Tech. Rep. EOH-11, RCA, 1974. RCA Technical Series Publication.

    R. D. Hippenstiel, Detection Theory: Applications and Digital Signal Pro-cessing, CRC Press, 2002

    Args:
        | ThresholdToNoise (float): the threshold to noise ratio [-]
        | pD (float): the probability of detection [-]

    Returns:
        | range (float): signal to noise ratio

    Raises:
        | No exception is raised.
    """

    import scipy.special

    SignalToNoise = np.sqrt(2) * scipy.special.erfinv(2 * pD -1) + ThresholdToNoise

    return SignalToNoise

##############################################################################
## 

def _bas(ord, angle, COS, TRI):
    """Define Basis vectors for a given polynomial order "order" and a
       given projection angle "angle".

    """

    basis_vec = scipy.special.eval_legendre(ord, angle) *\
                scipy.special.eval_legendre(ord, COS) * TRI
    return basis_vec 

def iv(*args, **kwargs):
        from scipy.special import iv
        return iv(*args, **kwargs) 

def __cdf_percentile(self, params, percentile, backscatter_model):
        mode = params[0]
        k = params[1]
        offset = params[2]
        if backscatter_model == RadarHistogram.BACKSCATTER_MODEL_GAUSSIAN:
            return scipy.special.erfinv(percentile / 0.5 - 1) * k * math.sqrt(2) + mode
        theta = (mode - offset) / (k - 1)
        v = scipy.special.gammaincinv(k, percentile) * theta + offset
        return v 

def __cdf(self, params, x, backscatter_model):
        mode = params[0]
        k = params[1]
        offset = params[2]
        if backscatter_model == RadarHistogram.BACKSCATTER_MODEL_GAUSSIAN:
            return 0.5 * (1 + scipy.special.erf((x - mode) / (k * math.sqrt(2))))
        theta = (mode - offset) / (k - 1)
        return scipy.special.gammainc(k, (x - offset) / theta)

    # find x where __cdf(params, offset, x) = percentile 

def erf(*args, **kwargs):
            from scipy.special import erf
            return erf(*args, **kwargs)


#    from scipy.special import lambertw, ellipe, gammaincc, gamma # fluids
#    from scipy.special import i1, i0, k1, k0, iv # ht
#    from scipy.special import hyp2f1    
#    if erf is None:
#        from scipy.special import erf 

def ellipkinc(phi, m):
        from scipy.special import ellipkinc
        return ellipkinc(phi, m) 

def hyp2f1(*args, **kwargs):
        from scipy.special import hyp2f1
        return hyp2f1(*args, **kwargs) 

def k0(*args, **kwargs):
        from scipy.special import k0
        return k0(*args, **kwargs)