Python mpmath.inf() Examples

def compute_b_mp(self, sigma, q, order, verbose=False):
    """Compute B_lambda for arbitrary lambda by numerical integration."""
    mu0, _, mu = self._distributions_mp(sigma, q)
    b_lambda_fn = lambda z: mu0(z) * (mu0(z) / mu(z)) ** order
    b_numeric = self._integral_inf_mp(b_lambda_fn)

    if verbose:
      _, z1 = rdp_accountant._compute_zs(sigma, q)
      print("z1 = ", z1)
      print("x in the Taylor series = ", q / (1 - q) * np.exp(
          (2 * z1 - 1) / (2 * sigma ** 2)))

      b0_numeric = self._integral_bounded_mp(b_lambda_fn, -np.inf, z1)
      b1_numeric = self._integral_bounded_mp(b_lambda_fn, z1, +np.inf)

      print("B: numerically {} = {} + {}".format(b_numeric, b0_numeric,
                                                 b1_numeric))
    return float(b_numeric) 

def test_wrightomega_branch():
    x = -np.logspace(10, 0, 25)
    picut_above = [np.nextafter(np.pi, np.inf)]
    picut_below = [np.nextafter(np.pi, -np.inf)]
    npicut_above = [np.nextafter(-np.pi, np.inf)]
    npicut_below = [np.nextafter(-np.pi, -np.inf)]
    for i in range(50):
        picut_above.append(np.nextafter(picut_above[-1], np.inf))
        picut_below.append(np.nextafter(picut_below[-1], -np.inf))
        npicut_above.append(np.nextafter(npicut_above[-1], np.inf))
        npicut_below.append(np.nextafter(npicut_below[-1], -np.inf))
    y = np.hstack((picut_above, picut_below, npicut_above, npicut_below))
    x, y = np.meshgrid(x, y)
    z = (x + 1j*y).flatten()

    dataset = []
    for z0 in z:
        dataset.append((z0, complex(_mpmath_wrightomega(z0, 25))))
    dataset = np.asarray(dataset)

    FuncData(sc.wrightomega, dataset, 0, 1, rtol=1e-8).check() 

def test_e1_complex(self):
        # E_1 oscillates as Im[z] -> +- inf, so limit range
        assert_mpmath_equal(sc.exp1,
                            mpmath.e1,
                            [ComplexArg(complex(-np.inf, -1e8), complex(np.inf, 1e8))],
                            rtol=1e-11)

        # Check cross-over region
        assert_mpmath_equal(sc.exp1,
                            mpmath.e1,
                            (np.linspace(-50, 50, 171)[:, None] +
                             np.r_[0, np.logspace(-3, 2, 61),
                                   -np.logspace(-3, 2, 11)]*1j).ravel(),
                            rtol=1e-11)
        assert_mpmath_equal(sc.exp1,
                            mpmath.e1,
                            (np.linspace(-50, -35, 10000) + 0j),
                            rtol=1e-11) 

def test_lanczos_sum_expg_scaled(self):
        maxgamma = 171.624376956302725
        e = np.exp(1)
        g = 6.024680040776729583740234375

        def gamma(x):
            with np.errstate(over='ignore'):
                fac = ((x + g - 0.5)/e)**(x - 0.5)
                if fac != np.inf:
                    res = fac*_lanczos_sum_expg_scaled(x)
                else:
                    fac = ((x + g - 0.5)/e)**(0.5*(x - 0.5))
                    res = fac*_lanczos_sum_expg_scaled(x)
                    res *= fac
            return res

        assert_mpmath_equal(gamma,
                            mpmath.gamma,
                            [Arg(0, maxgamma, inclusive_a=False)],
                            rtol=1e-13) 

def _compute_delta(log_moments, eps):
  """Compute delta for given log_moments and eps.
  Args:
    log_moments: the log moments of privacy loss, in the form of pairs
      of (moment_order, log_moment)
    eps: the target epsilon.
  Returns:
    delta
  """
  min_delta = 1.0
  for moment_order, log_moment in log_moments:
    if moment_order == 0:
      continue
    if math.isinf(log_moment) or math.isnan(log_moment):
      sys.stderr.write("The %d-th order is inf or Nan\n" % moment_order)
      continue
    if log_moment < moment_order * eps:
      min_delta = min(min_delta,
                      math.exp(log_moment - moment_order * eps))
  return min_delta 

def _compute_eps(log_moments, delta):
  """Compute epsilon for given log_moments and delta.
  Args:
    log_moments: the log moments of privacy loss, in the form of pairs
      of (moment_order, log_moment)
    delta: the target delta.
  Returns:
    epsilon
  """
  min_eps = float("inf")
  for moment_order, log_moment in log_moments:
    if moment_order == 0:
      continue
    if math.isinf(log_moment) or math.isnan(log_moment):
      sys.stderr.write("The %d-th order is inf or Nan\n" % moment_order)
      continue
    min_eps = min(min_eps, (log_moment - math.log(delta)) / moment_order)
  return min_eps 

def test_compute_rdp_sequence(self):
    rdp_vec = rdp_accountant.compute_rdp(0.01, 2.5, 50,
                                         [1.5, 2.5, 5, 50, 100, np.inf])
    self.assertSequenceAlmostEqual(
        rdp_vec, [0.00065, 0.001085, 0.00218075, 0.023846, 167.416307, np.inf],
        delta=1e-5) 

def _to_np_float64(v):
  if math.isnan(v) or math.isinf(v):
    return np.inf
  return np.float64(v)


######################
# FLOAT64 ARITHMETIC #
###################### 

def integral_inf(fn):
  integral, _ = integrate.quad(fn, -np.inf, np.inf)
  return integral 

def integral_inf_mp(fn):
  integral, _ = mp.quad(fn, [-mp.inf, mp.inf], error=True)
  return integral 

def _integral_mp(self, fn, bounds=(-inf, inf)):
    integral, _ = quad(fn, bounds, error=True, maxdegree=8)
    return integral 

def test_si_complex(self):
        def si(z):
            return sc.sici(z)[0]
        # si oscillates as Re[z] -> +- inf, so limit range
        assert_mpmath_equal(si,
                            mpmath.si,
                            [ComplexArg(complex(-1e8, -np.inf), complex(1e8, np.inf))],
                            rtol=1e-12) 

def _integral_mp(self, fn, bounds=(-mp.inf, mp.inf)):
    integral, _ = mp.quad(fn, bounds, error=True, maxdegree=8)
    return integral 

def test_compute_rdp_sequence(self):
    rdp_vec = rdp_accountant.compute_rdp(0.01, 2.5, 50,
                                         [1.5, 2.5, 5, 50, 100, np.inf])
    self.assertSequenceAlmostEqual(
        rdp_vec, [0.00065, 0.001085, 0.00218075, 0.023846, 167.416307, np.inf],
        delta=1e-5) 

def test_nstr():
    m = matrix([[0.75, 0.190940654, -0.0299195971],
                [0.190940654, 0.65625, 0.205663228],
                [-0.0299195971, 0.205663228, 0.64453125e-20]])
    assert nstr(m, 4, min_fixed=-inf) == \
    '''[    0.75  0.1909                    -0.02992]
[  0.1909  0.6563                      0.2057]
[-0.02992  0.2057  0.000000000000000000006445]'''
    assert nstr(m, 4) == \
    '''[    0.75  0.1909   -0.02992]
[  0.1909  0.6563     0.2057]
[-0.02992  0.2057  6.445e-21]''' 

def _integral_inf_mp(self, fn):
    integral, _ = mp.quad(
        fn, [-mp.inf, mp.inf], error=True,
        maxdegree=7)  # maxdegree = 6 is not enough
    return integral 

def test_compute_rdp(self):
    rdp_scalar = rdp_accountant.compute_rdp(0.1, 2, 10, 5)
    self.assertAlmostEqual(rdp_scalar, 0.07737, places=5)

    rdp_vec = rdp_accountant.compute_rdp(0.01, 2.5, 50, [1.5, 2.5, 5, 50, 100,
                                                         np.inf])

    correct = [0.00065, 0.001085, 0.00218075, 0.023846, 167.416307, np.inf]
    for i in range(len(rdp_vec)):
      self.assertAlmostEqual(rdp_vec[i], correct[i], places=5) 

def test_spence(self):
        # mpmath uses a different convention for the dilogarithm
        def dilog(x):
            return mpmath.polylog(2, 1 - x)
        # Spence has a branch cut on the negative real axis
        assert_mpmath_equal(sc.spence,
                            exception_to_nan(dilog),
                            [Arg(0, np.inf)], rtol=1e-14) 

def test_rgamma(self):
        def rgamma(x):
            if x < -8000:
                return np.inf
            else:
                v = mpmath.rgamma(x)
            return v
        # n=500 (non-xslow default) fails for one bad point
        assert_mpmath_equal(sc.rgamma,
                            rgamma,
                            [Arg()],
                            n=5000,
                            ignore_inf_sign=True) 

def test_lambertw_real(self):
        assert_mpmath_equal(lambda x, k: sc.lambertw(x, int(k.real)),
                            lambda x, k: mpmath.lambertw(x, int(k.real)),
                            [ComplexArg(-np.inf, np.inf), IntArg(0, 10)],
                            rtol=1e-13, nan_ok=False) 

def test_ei_complex(self):
        # Ei oscillates as Im[z] -> +- inf, so limit range
        assert_mpmath_equal(sc.expi,
                            mpmath.ei,
                            [ComplexArg(complex(-np.inf, -1e8), complex(np.inf, 1e8))],
                            rtol=1e-9) 

def test_expm1_complex(self):
        # Oscillates as a function of Im[z], so limit range to avoid loss of precision
        assert_mpmath_equal(sc.expm1,
                            mpmath.expm1,
                            [ComplexArg(complex(-np.inf, -1e7), complex(np.inf, 1e7))]) 

def test_exprel(self):
        assert_mpmath_equal(sc.exprel,
                            lambda x: mpmath.expm1(x)/x if x != 0 else mpmath.mpf('1.0'),
                            [Arg(a=-np.log(np.finfo(np.double).max), b=np.log(np.finfo(np.double).max))])
        assert_mpmath_equal(sc.exprel,
                            lambda x: mpmath.expm1(x)/x if x != 0 else mpmath.mpf('1.0'),
                            np.array([1e-12, 1e-24, 0, 1e12, 1e24, np.inf]), rtol=1e-11)
        assert_(np.isinf(sc.exprel(np.inf)))
        assert_(sc.exprel(-np.inf) == 0) 

def test_ci_complex(self):
        def ci(z):
            return sc.sici(z)[1]
        # ci oscillates as Re[z] -> +- inf, so limit range
        assert_mpmath_equal(ci,
                            mpmath.ci,
                            [ComplexArg(complex(-1e8, -np.inf), complex(1e8, np.inf))],
                            rtol=1e-8) 

def test_bessely_complex(self):
        def mpbessely(v, x):
            r = complex(mpmath.bessely(v, x, **HYPERKW))
            if abs(r) > 1e305:
                # overflowing to inf a bit earlier is OK
                olderr = np.seterr(invalid='ignore')
                try:
                    r = np.inf * np.sign(r)
                finally:
                    np.seterr(**olderr)
            return r
        assert_mpmath_equal(lambda v, z: sc.yv(v.real, z),
                            exception_to_nan(mpbessely),
                            [Arg(), ComplexArg()],
                            n=15000) 

def test_bessely(self):
        def mpbessely(v, x):
            r = float(mpmath.bessely(v, x, **HYPERKW))
            if abs(r) > 1e305:
                # overflowing to inf a bit earlier is OK
                r = np.inf * np.sign(r)
            if abs(r) == 0 and x == 0:
                # invalid result from mpmath, point x=0 is a divergence
                return np.nan
            return r
        assert_mpmath_equal(sc.yv,
                            exception_to_nan(mpbessely),
                            [Arg(-1e100, 1e100), Arg(-1e8, 1e8)],
                            n=5000) 

def test_besselk_int(self):
        assert_mpmath_equal(sc.kn,
                            mpmath.besselk,
                            [IntArg(-200, 200), Arg(0, np.inf)],
                            nan_ok=False, rtol=1e-12) 

def test_besselk(self):
        assert_mpmath_equal(sc.kv,
                            mpmath.besselk,
                            [Arg(-200, 200), Arg(0, np.inf)],
                            nan_ok=False, rtol=1e-12) 

def test_legenp(self):
        def lpnm(n, m, z):
            try:
                v = sc.lpmn(m, n, z)[0][-1,-1]
            except ValueError:
                return np.nan
            if abs(v) > 1e306:
                # harmonize overflow to inf
                v = np.inf * np.sign(v.real)
            return v

        def lpnm_2(n, m, z):
            v = sc.lpmv(m, n, z)
            if abs(v) > 1e306:
                # harmonize overflow to inf
                v = np.inf * np.sign(v.real)
            return v

        def legenp(n, m, z):
            if (z == 1 or z == -1) and int(n) == n:
                # Special case (mpmath may give inf, we take the limit by
                # continuity)
                if m == 0:
                    if n < 0:
                        n = -n - 1
                    return mpmath.power(mpmath.sign(z), n)
                else:
                    return 0

            if abs(z) < 1e-15:
                # mpmath has bad performance here
                return np.nan

            typ = 2 if abs(z) < 1 else 3
            v = exception_to_nan(mpmath.legenp)(n, m, z, type=typ)

            if abs(v) > 1e306:
                # harmonize overflow to inf
                v = mpmath.inf * mpmath.sign(v.real)

            return v

        assert_mpmath_equal(lpnm,
                            legenp,
                            [IntArg(-100, 100), IntArg(-100, 100), Arg()])

        assert_mpmath_equal(lpnm_2,
                            legenp,
                            [IntArg(-100, 100), Arg(-100, 100), Arg(-1, 1)],
                            atol=1e-10) 

def test_gegenbauer_int(self):
        # Redefine functions to deal with numerical + mpmath issues
        def gegenbauer(n, a, x):
            # Avoid overflow at large `a` (mpmath would need an even larger
            # dps to handle this correctly, so just skip this region)
            if abs(a) > 1e100:
                return np.nan

            # Deal with n=0, n=1 correctly; mpmath 0.17 doesn't do these
            # always correctly
            if n == 0:
                r = 1.0
            elif n == 1:
                r = 2*a*x
            else:
                r = mpmath.gegenbauer(n, a, x)

            # Mpmath 0.17 gives wrong results (spurious zero) in some cases, so
            # compute the value by perturbing the result
            if float(r) == 0 and a < -1 and float(a) == int(float(a)):
                r = mpmath.gegenbauer(n, a + mpmath.mpf('1e-50'), x)
                if abs(r) < mpmath.mpf('1e-50'):
                    r = mpmath.mpf('0.0')

            # Differing overflow thresholds in scipy vs. mpmath
            if abs(r) > 1e270:
                return np.inf
            return r

        def sc_gegenbauer(n, a, x):
            r = sc.eval_gegenbauer(int(n), a, x)
            # Differing overflow thresholds in scipy vs. mpmath
            if abs(r) > 1e270:
                return np.inf
            return r
        assert_mpmath_equal(sc_gegenbauer,
                            exception_to_nan(gegenbauer),
                            [IntArg(0, 100), Arg(-1e9, 1e9), Arg()],
                            n=40000, dps=100,
                            ignore_inf_sign=True, rtol=1e-6)

        # Check the small-x expansion
        assert_mpmath_equal(sc_gegenbauer,
                            exception_to_nan(gegenbauer),
                            [IntArg(0, 100), Arg(), FixedArg(np.logspace(-30, -4, 30))],
                            dps=100,
                            ignore_inf_sign=True)