Python scipy.special.iv() Examples

def kaiser_bessel_ft(om, npts, alpha, order, d):
    """Computes FT of KB function for scaling in image domain.

    Args:
        om (ndarray): An array of coordinates to interpolate to.
        npts (int): Number of points to use for interpolation in each
            dimension.
        order (ind, default=0): Order of Kaiser-Bessel kernel.
        alpha (double or array of doubles): KB parameter.

    Returns:
        ndarray: The scaling coefficients.
    """
    z = np.sqrt((2 * np.pi * (npts / 2) * om)**2 - alpha**2 + 0j)
    nu = d / 2 + order
    scaling_coef = (2 * np.pi)**(d / 2) * ((npts / 2)**d) * (alpha**order) / \
        special.iv(order, alpha) * special.jv(nu, z) / (z**nu)
    scaling_coef = np.real(scaling_coef)

    return scaling_coef 

def KL_davidson(k, d):
    vmf_entropy = k * ive(d / 2, k) / ive((d / 2) - 1, k) + \
                  (d / 2 - 1) * np.log(k) \
                  - (d / 2) * np.log(2 * np.pi) - np.log(iv(d / 2 - 1, k))

    hyu_ent = np.log(2) + (d / 2) * np.log(np.pi) - sp.loggamma(
        d / 2)

    kl = vmf_entropy + hyu_ent
    return kl
#
# first = k * bessel(d / 2, k) / bessel(d / 2 - 1, k)
# second = (d / 2 - 1) * torch.log(k) - torch.log(bessel(d / 2 - 1, k))
# const = torch.tensor(
#            np.log(3.1415926) * d / 2 + np.log(2) - sp.loggamma(d / 2).real - (d / 2) * np.log(2 * 3.1415926)).to(
#            devic

# for kappa in range(10, 150, 20):
#     for d in range(50, 150, 50):
#         print("Davidson:{}\t\tGuu:{}".format(KL_davidson(kappa, d), KL_guu(kappa, d))) 

def test_ticket_854(self):
        """Real-valued Bessel domains"""
        assert_(isnan(special.jv(0.5, -1)))
        assert_(isnan(special.iv(0.5, -1)))
        assert_(isnan(special.yv(0.5, -1)))
        assert_(isnan(special.yv(1, -1)))
        assert_(isnan(special.kv(0.5, -1)))
        assert_(isnan(special.kv(1, -1)))
        assert_(isnan(special.jve(0.5, -1)))
        assert_(isnan(special.ive(0.5, -1)))
        assert_(isnan(special.yve(0.5, -1)))
        assert_(isnan(special.yve(1, -1)))
        assert_(isnan(special.kve(0.5, -1)))
        assert_(isnan(special.kve(1, -1)))
        assert_(isnan(special.airye(-1)[0:2]).all(), special.airye(-1))
        assert_(not isnan(special.airye(-1)[2:4]).any(), special.airye(-1)) 

def nena(locs, info, callback=None):
    bin_centers, dnfl_ = next_frame_neighbor_distance_histogram(locs, callback)

    def func(d, a, s, ac, dc, sc):
        f = a * (d / s ** 2) * _np.exp(-0.5 * d ** 2 / s ** 2)
        fc = (
            ac
            * (d / sc ** 2)
            * _np.exp(-0.5 * (d ** 2 + dc ** 2) / sc ** 2)
            * _iv(0, d * dc / sc)
        )
        return f + fc

    pdf_model = _lmfit.Model(func)
    params = _lmfit.Parameters()
    area = _np.trapz(dnfl_, bin_centers)
    median_lp = _np.mean([_np.median(locs.lpx), _np.median(locs.lpy)])
    params.add("a", value=area / 2, min=0)
    params.add("s", value=median_lp, min=0)
    params.add("ac", value=area / 2, min=0)
    params.add("dc", value=2 * median_lp, min=0)
    params.add("sc", value=median_lp, min=0)
    result = pdf_model.fit(dnfl_, params, d=bin_centers)
    return result, result.best_values["s"] 

def test_ticket_854(self):
        """Real-valued Bessel domains"""
        assert_(isnan(special.jv(0.5, -1)))
        assert_(isnan(special.iv(0.5, -1)))
        assert_(isnan(special.yv(0.5, -1)))
        assert_(isnan(special.yv(1, -1)))
        assert_(isnan(special.kv(0.5, -1)))
        assert_(isnan(special.kv(1, -1)))
        assert_(isnan(special.jve(0.5, -1)))
        assert_(isnan(special.ive(0.5, -1)))
        assert_(isnan(special.yve(0.5, -1)))
        assert_(isnan(special.yve(1, -1)))
        assert_(isnan(special.kve(0.5, -1)))
        assert_(isnan(special.kve(1, -1)))
        assert_(isnan(special.airye(-1)[0:2]).all(), special.airye(-1))
        assert_(not isnan(special.airye(-1)[2:4]).any(), special.airye(-1)) 

def _cephes_vs_amos_points(self):
        """Yield points at which to compare Cephes implementation to AMOS"""
        # check several points, including large-amplitude ones
        for v in [-120, -100.3, -20., -10., -1., -.5,
                  0., 1., 12.49, 120., 301]:
            for z in [-1300, -11, -10, -1, 1., 10., 200.5, 401., 600.5,
                      700.6, 1300, 10003]:
                yield v, z

        # check half-integers; these are problematic points at least
        # for cephes/iv
        for v in 0.5 + arange(-60, 60):
            yield v, 3.5 

def _pdf(self, x, df, nc):
        output = 0.5*numpy.e**(-0.5*(x+nc))
        output *= (x/nc)**(0.25*df-0.5)
        output *= special.iv(0.5*df-1, (nc*x)**0.5)
        return output 

def test_iv_cephes_vs_amos(self):
        olderr = np.seterr(all='ignore')
        try:
            self.check_cephes_vs_amos(special.iv, special.iv, rtol=5e-9, atol=1e-305)
        finally:
            np.seterr(**olderr) 

def test_iv_cephes_vs_amos_mass_test(self):
        N = 1000000
        np.random.seed(1)
        v = np.random.pareto(0.5, N) * (-1)**np.random.randint(2, size=N)
        x = np.random.pareto(0.2, N) * (-1)**np.random.randint(2, size=N)

        imsk = (np.random.randint(8, size=N) == 0)
        v[imsk] = v[imsk].astype(int)

        old_err = np.seterr(all='ignore')
        try:
            c1 = special.iv(v, x)
            c2 = special.iv(v, x+0j)

            # deal with differences in the inf and zero cutoffs
            c1[abs(c1) > 1e300] = np.inf
            c2[abs(c2) > 1e300] = np.inf
            c1[abs(c1) < 1e-300] = 0
            c2[abs(c2) < 1e-300] = 0

            dc = abs(c1/c2 - 1)
            dc[np.isnan(dc)] = 0
        finally:
            np.seterr(**old_err)

        k = np.argmax(dc)

        # Most error apparently comes from AMOS and not our implementation;
        # there are some problems near integer orders there
        assert_(dc[k] < 2e-7, (v[k], x[k], special.iv(v[k], x[k]), special.iv(v[k], x[k]+0j))) 

def test_ticket_503(self):
        """Real-valued Bessel I overflow"""
        assert_allclose(special.iv(1, 700), 1.528500390233901e302)
        assert_allclose(special.iv(1000, 1120), 1.301564549405821e301) 

def test_iv_hyperg_poles(self):
        assert_allclose(special.iv(-0.5, 1), 1.231200214592967) 

def test_iv_series(self):
        for v in [-20., -10., -1., 0., 1., 12.49, 120.]:
            for z in [1., 10., 200.5, -1+2j]:
                value, err = self.iv_series(v, z)
                assert_allclose(special.iv(v, z), value, atol=err, err_msg=(v, z)) 

def test_iv(self):
        iv1 = special.iv(0,.1)*exp(-.1)
        assert_almost_equal(iv1,0.90710092578230106,10) 

def test_ivp0(self):
        assert_almost_equal(special.iv(1,2), special.ivp(0,2), 10) 

def test_ivp(self):
        y = (special.iv(0,2) + special.iv(2,2))/2
        x = special.ivp(1,2)
        assert_almost_equal(x,y,10) 

def test_error_raising():
    assert_raises(special.SpecialFunctionError, special.iv, 1, 1e99j) 

def test_besseli(self):
        assert_mpmath_equal(sc.iv,
                            exception_to_nan(lambda v, z: mpmath.besseli(v, z, **HYPERKW)),
                            [Arg(-1e100, 1e100), Arg()],
                            atol=1e-270) 

def test_besseli_complex(self):
        assert_mpmath_equal(lambda v, z: sc.iv(v.real, z),
                            exception_to_nan(lambda v, z: mpmath.besseli(v, z, **HYPERKW)),
                            [Arg(-1e100, 1e100), ComplexArg()]) 

def _chebyshev_gaussian(self, A, a, order=10, mu=0.5, s=0.5):
        # NE Enhancement via Spectral Propagation
        print("Chebyshev Series -----------------")
        t1 = time.time()

        if order == 1:
            return a

        A = sp.eye(self.num_node) + A
        DA = preprocessing.normalize(A, norm="l1")
        L = sp.eye(self.num_node) - DA

        M = L - mu * sp.eye(self.num_node)

        Lx0 = a
        Lx1 = M.dot(a)
        Lx1 = 0.5 * M.dot(Lx1) - a

        conv = iv(0, s) * Lx0
        conv -= 2 * iv(1, s) * Lx1
        for i in range(2, order):
            Lx2 = M.dot(Lx1)
            Lx2 = (M.dot(Lx2) - 2 * Lx1) - Lx0
            #         Lx2 = 2*L.dot(Lx1) - Lx0
            if i % 2 == 0:
                conv += 2 * iv(i, s) * Lx2
            else:
                conv -= 2 * iv(i, s) * Lx2
            Lx0 = Lx1
            Lx1 = Lx2
            del Lx2
            print("Bessell time", i, time.time() - t1)
        mm = A.dot(a - conv)
        emb = self._get_embedding_dense(mm, self.dimension)
        return emb 

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

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 iv(*args, **kwargs):
        import mpmath
        return mpmath.besseli(*args, **kwargs) 

def iv(*args, **kwargs):
        import mpmath
        return mpmath.besseli(*args, **kwargs) 

def _kaiser(n, beta):
    """Independant Kaiser window

    For the definition of the Kaiser window, see A. V. Oppenheim & R. W. Schafer, "Discrete-Time Signal Processing".

    The continuous version of width n centered about x=0 is:

    .. note:: 2 times slower than scipy.kaiser
    """
    from scipy.special import iv as besselI
    m = n - 1
    k = arange(0, m)
    k = 2. * beta / m * sqrt (k * (m - k))
    w = besselI (0, k) / besselI (0, beta)
    return w 

def rt1d_coda(r, t, c, g0):
    from scipy.special import iv
    l = 1 / g0
    arg0 = np.sqrt(c ** 2 * t ** 2 - r ** 2)
    arg1 = arg0 / 2 / l
    t1 = 1 / 4 / l * np.exp(-c * t / 2 / l)
    t2 = iv(0, arg1) + c * t * iv(1, arg1) / arg0
    return t1 * t2  # * Heaviside(c * t - r) 

def chebyshev_gaussian(self, A, a, order=10, mu=0.5, s=0.5):
		# NE Enhancement via Spectral Propagation
		print('Chebyshev Series -----------------')
		t1 = time.time()

		if order == 1:
			return a

		A = sp.eye(self.node_number) + A
		DA = preprocessing.normalize(A, norm='l1')
		L = sp.eye(self.node_number) - DA

		M = L - mu * sp.eye(self.node_number)

		Lx0 = a
		Lx1 = M.dot(a)
		Lx1 = 0.5 * M.dot(Lx1) - a

		conv = iv(0, s) * Lx0
		conv -= 2 * iv(1, s) * Lx1
		for i in range(2, order):
			Lx2 = M.dot(Lx1)
			Lx2 = (M.dot(Lx2) - 2 * Lx1) - Lx0
			#         Lx2 = 2*L.dot(Lx1) - Lx0
			if i % 2 == 0:
				conv += 2 * iv(i, s) * Lx2
			else:
				conv -= 2 * iv(i, s) * Lx2
			Lx0 = Lx1
			Lx1 = Lx2
			del Lx2
			print('Bessell time', i, time.time() - t1)
		mm = A.dot(a - conv)
		emb = self.get_embedding_dense(mm, self.dimension)
		return emb 

def _chebyshev_gaussian(self, A, a, order=5, mu=0.5, s=0.2, plus=False, nn=False):
        # NE Enhancement via Spectral Propagation
        print("Chebyshev Series -----------------")
        t1 = time.time()
        num_node = a.shape[0]

        if order == 1:
            return a

        A = sp.eye(num_node) + A
        DA = preprocessing.normalize(A, norm="l1")
        L = sp.eye(num_node) - DA

        M = L - mu * sp.eye(num_node)

        Lx0 = a
        Lx1 = M.dot(a)
        Lx1 = 0.5 * M.dot(Lx1) - a

        conv = iv(0, s) * Lx0
        conv -= 2 * iv(1, s) * Lx1
        for i in range(2, order):
            Lx2 = M.dot(Lx1)
            Lx2 = (M.dot(Lx2) - 2 * Lx1) - Lx0
            #         Lx2 = 2*L.dot(Lx1) - Lx0
            if i % 2 == 0:
                conv += 2 * iv(i, s) * Lx2
            else:
                conv -= 2 * iv(i, s) * Lx2
            Lx0 = Lx1
            Lx1 = Lx2
            del Lx2
            print("Bessell time", i, time.time() - t1)
        emb = mm = conv
        if not plus:
            mm = A.dot(a - conv)
        if not nn:
            emb = self._get_embedding_dense(mm, self.dimension)
        return emb 

def _G_integrand(xip, B):
    try:
        from scipy.special import iv as besseli
    except:
        """ [PHIDL] To run this function you need scipy, please install it with
        pip install scipy """
    return besseli(0, B*sqrt(1-xip**2)) 

def test_iv_cephes_vs_amos(self):
        olderr = np.seterr(all='ignore')
        try:
            self.check_cephes_vs_amos(special.iv, special.iv, rtol=5e-9, atol=1e-305)
        finally:
            np.seterr(**olderr) 

def forward(ctx, dim, kappa):
        """
        In the forward pass we receive a Tensor containing the input and return
        a Tensor containing the output. ctx is a context object that can be used
        to stash information for backward computation. You can cache arbitrary
        objects for use in the backward pass using the ctx.save_for_backward method.
        """
        ctx.save_for_backward(dim, kappa)

        kappa_copy = kappa.clone()
        m = sp.iv(dim, kappa_copy)
        x = torch.tensor(m).to(device)
        return x.clone()