Python mpmath.log() Examples

def test_boxcox(self):

        def mp_boxcox(x, lmbda):
            x = mpmath.mp.mpf(x)
            lmbda = mpmath.mp.mpf(lmbda)
            if lmbda == 0:
                return mpmath.mp.log(x)
            else:
                return mpmath.mp.powm1(x, lmbda) / lmbda

        assert_mpmath_equal(sc.boxcox,
                            exception_to_nan(mp_boxcox),
                            [Arg(a=0, inclusive_a=False), Arg()],
                            n=200,
                            dps=60,
                            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 Swamee_Jain_1976(Re, eD):
    r'''Calculates Darcy friction factor using the method in Swamee and
    Jain (1976) [2]_ as shown in [1]_.

    .. math::
        \frac{1}{\sqrt{f_f}} = -4\log\left[\left(\frac{6.97}{Re}\right)^{0.9}
        + (\frac{\epsilon}{3.7D})\right]

    Parameters
    ----------
    Re : float
        Reynolds number, [-]
    eD : float
        Relative roughness, [-]

    Returns
    -------
    fd : float
        Darcy friction factor [-]

    Notes
    -----
    Range is 5E3 <= Re <= 1E8;  1E-6 <= eD <= 5E-2.

    Examples
    --------
    >>> Swamee_Jain_1976(1E5, 1E-4)
    0.018452424431901808

    References
    ----------
    .. [1] Winning, H. and T. Coole. "Explicit Friction Factor Accuracy and
       Computational Efficiency for Turbulent Flow in Pipes." Flow, Turbulence
       and Combustion 90, no. 1 (January 1, 2013): 1-27.
       doi:10.1007/s10494-012-9419-7
    .. [2] Swamee, Prabhata K., and Akalank K. Jain."Explicit Equations for
       Pipe-Flow Problems." Journal of the Hydraulics Division 102, no. 5
       (May 1976): 657-664.
    '''
    ff = (-4*log10((6.97/Re)**0.9 + eD/3.7))**-2
    return 4*ff 

def _compute_rdp_mp(self, q, sigma, order):
    log_a_mp = float(mp.log(self.compute_a_mp(sigma, q, order)))
    log_b_mp = float(mp.log(self.compute_b_mp(sigma, q, order)))
    return log_a_mp, log_b_mp

  # TEST ROUTINES 

def test_compute_log_a_equals_mp(self, q, sigma, order):
    # Compare the cheap computation of log(A) with an expensive, multi-precision
    # computation.
    log_a = rdp_accountant._compute_log_a(q, sigma, order)
    log_a_mp = self._log_float_mp(self._compute_a_mp(sigma, q, order))
    np.testing.assert_allclose(log_a, log_a_mp, rtol=1e-4) 

def _log_float_mp(self, x):
    # Convert multi-precision input to float log space.
    if x >= sys.float_info.min:
      return float(mp.log(x))
    else:
      return -np.inf 

def eta(lam):
    """Function from DLMF 8.12.1 shifted to be centered at 0."""
    if lam > 0:
        return mp.sqrt(2*(lam - mp.log(lam + 1)))
    elif lam < 0:
        return -mp.sqrt(2*(lam - mp.log(lam + 1)))
    else:
        return 0 

def test_compute_log_a_equals_mp(self, q, sigma, order):
    # Compare the cheap computation of log(A) with an expensive, multi-precision
    # computation.
    log_a = rdp_accountant._compute_log_a(q, sigma, order)
    log_a_mp = self._log_float_mp(self._compute_a_mp(sigma, q, order))
    np.testing.assert_allclose(log_a, log_a_mp, rtol=1e-4) 

def _log_float_mp(self, x):
    # Convert multi-precision input to float log space.
    if x >= sys.float_info.min:
      return float(log(x))
    else:
      return -np.inf 

def Avci_Karagoz_2009(Re, eD):
    r'''Calculates Darcy friction factor using the method in Avci and Karagoz
    (2009) [2]_ as shown in [1]_.

    .. math::
        f_D = \frac{6.4} {\left\{\ln(Re) - \ln\left[
        1 + 0.01Re\frac{\epsilon}{D}\left(1 + 10(\frac{\epsilon}{D})^{0.5}
        \right)\right]\right\}^{2.4}}

    Parameters
    ----------
    Re : float
        Reynolds number, [-]
    eD : float
        Relative roughness, [-]

    Returns
    -------
    fd : float
        Darcy friction factor [-]

    Notes
    -----
    No range of validity specified for this equation.

    Examples
    --------
    >>> Avci_Karagoz_2009(1E5, 1E-4)
    0.01857058061066499

    References
    ----------
    .. [1] Winning, H. and T. Coole. "Explicit Friction Factor Accuracy and
       Computational Efficiency for Turbulent Flow in Pipes." Flow, Turbulence
       and Combustion 90, no. 1 (January 1, 2013): 1-27.
       doi:10.1007/s10494-012-9419-7
    .. [2]	Avci, Atakan, and Irfan Karagoz."A Novel Explicit Equation for
       Friction Factor in Smooth and Rough Pipes." Journal of Fluids
       Engineering 131, no. 6 (2009): 061203. doi:10.1115/1.3129132.
    '''
    return 6.4*(log(Re) - log(1 + 0.01*Re*eD*(1+10*eD**0.5)))**-2.4 

def eta(lam):
    """Function from DLMF 8.12.1 shifted to be centered at 0."""
    if lam > 0:
        return mp.sqrt(2*(lam - mp.log(lam + 1)))
    elif lam < 0:
        return -mp.sqrt(2*(lam - mp.log(lam + 1)))
    else:
        return 0 

def Jain_1976(Re, eD):
    r'''Calculates Darcy friction factor using the method in Jain (1976)
    [2]_ as shown in [1]_.

    .. math::
        \frac{1}{\sqrt{f_f}} = 2.28 - 4\log\left[ \frac{\epsilon}{D} +
        \left(\frac{29.843}{Re}\right)^{0.9}\right]

    Parameters
    ----------
    Re : float
        Reynolds number, [-]
    eD : float
        Relative roughness, [-]

    Returns
    -------
    fd : float
        Darcy friction factor [-]

    Notes
    -----
    Range is 5E3 <= Re <= 1E7;  4E-5 <= eD <= 5E-2.

    Examples
    --------
    >>> Jain_1976(1E5, 1E-4)
    0.018436560312693327

    References
    ----------
    .. [1] Winning, H. and T. Coole. "Explicit Friction Factor Accuracy and
       Computational Efficiency for Turbulent Flow in Pipes." Flow, Turbulence
       and Combustion 90, no. 1 (January 1, 2013): 1-27.
       doi:10.1007/s10494-012-9419-7
    .. [2] 	Jain, Akalank K."Accurate Explicit Equation for Friction Factor."
       Journal of the Hydraulics Division 102, no. 5 (May 1976): 674-77.
    '''
    ff = (2.28-4*log10(eD+(29.843/Re)**0.9))**-2
    return 4*ff 

def Eck_1973(Re, eD):
    r'''Calculates Darcy friction factor using the method in Eck (1973)
    [2]_ as shown in [1]_.

    .. math::
        \frac{1}{\sqrt{f_d}} = -2\log\left[\frac{\epsilon}{3.715D}
        + \frac{15}{Re}\right]

    Parameters
    ----------
    Re : float
        Reynolds number, [-]
    eD : float
        Relative roughness, [-]

    Returns
    -------
    fd : float
        Darcy friction factor [-]

    Notes
    -----
    No range of validity specified for this equation.

    Examples
    --------
    >>> Eck_1973(1E5, 1E-4)
    0.01775666973488564

    References
    ----------
    .. [1] Winning, H. and T. Coole. "Explicit Friction Factor Accuracy and
       Computational Efficiency for Turbulent Flow in Pipes." Flow, Turbulence
       and Combustion 90, no. 1 (January 1, 2013): 1-27.
       doi:10.1007/s10494-012-9419-7
    .. [2] Eck, B.: Technische Stromungslehre. Springer, New York (1973)
    '''
    return (-2*log10(eD/3.715 + 15/Re))**-2 

def Churchill_1973(Re, eD):
    r'''Calculates Darcy friction factor using the method in Churchill (1973)
    [2]_ as shown in [1]_

    .. math::
        \frac{1}{\sqrt{f_d}} = -2\log\left[\frac{\epsilon}{3.7D} +
        (\frac{7}{Re})^{0.9}\right]

    Parameters
    ----------
    Re : float
        Reynolds number, [-]
    eD : float
        Relative roughness, [-]

    Returns
    -------
    fd : float
        Darcy friction factor [-]

    Notes
    -----
    No range of validity specified for this equation.

    Examples
    --------
    >>> Churchill_1973(1E5, 1E-4)
    0.01846708694482294

    References
    ----------
    .. [1] Winning, H. and T. Coole. "Explicit Friction Factor Accuracy and
       Computational Efficiency for Turbulent Flow in Pipes." Flow, Turbulence
       and Combustion 90, no. 1 (January 1, 2013): 1-27.
       doi:10.1007/s10494-012-9419-7
    .. [2] Churchill, Stuart W. "Empirical Expressions for the Shear
       Stress in Turbulent Flow in Commercial Pipe." AIChE Journal 19, no. 2
       (March 1, 1973): 375-76. doi:10.1002/aic.690190228.
    '''
    return (-2*log10(eD/3.7 + (7./Re)**0.9))**-2 

def entropy(prob_vector, logbase=2.):
    """Entropy of a probability vector that sums too one."""
    prob_vector = np.array(prob_vector)
    pos_prob_vector = prob_vector[prob_vector > 0]
    return - np.sum(pos_prob_vector * np.log(pos_prob_vector)/np.log(logbase)) 

def compute_b_mp(sigma, q, lmbd, verbose=False):
  lmbd_int = int(math.ceil(lmbd))
  if lmbd_int == 0:
    return 1.0

  mu0, _, mu = distributions_mp(sigma, q)

  b_lambda_fn = lambda z: mu0(z) * (mu0(z) / mu(z)) ** lmbd_int
  b_lambda = integral_inf_mp(b_lambda_fn)

  m = sigma ** 2 * (mp.log((2 - q) / (1 - q)) + 1 / (2 * (sigma ** 2)))
  b_fn = lambda z: ((mu0(z) / mu(z)) ** lmbd_int -
                    (mu(-z) / mu0(z)) ** lmbd_int)
  if verbose:
    print "M =", m
    print "f(-M) = {} f(M) = {}".format(b_fn(-m), b_fn(m))
    assert b_fn(-m) < 0 and b_fn(m) < 0

  b_lambda_int1_fn = lambda z: mu0(z) * (mu0(z) / mu(z)) ** lmbd_int
  b_lambda_int2_fn = lambda z: mu0(z) * (mu(z) / mu0(z)) ** lmbd_int
  b_int1 = integral_bounded_mp(b_lambda_int1_fn, -m, m)
  b_int2 = integral_bounded_mp(b_lambda_int2_fn, -m, m)

  a_lambda_m1 = compute_a_mp(sigma, q, lmbd - 1)
  b_bound = a_lambda_m1 + b_int1 - b_int2

  if verbose:
    print "B by numerical integration", b_lambda
    print "B must be no more than    ", b_bound
  assert b_lambda < b_bound + 1e-5
  return _to_np_float64(b_lambda) 

def test_beta():
    np.random.seed(1234)

    b = np.r_[np.logspace(-200, 200, 4),
              np.logspace(-10, 10, 4),
              np.logspace(-1, 1, 4),
              np.arange(-10, 11, 1),
              np.arange(-10, 11, 1) + 0.5,
              -1, -2.3, -3, -100.3, -10003.4]
    a = b

    ab = np.array(np.broadcast_arrays(a[:,None], b[None,:])).reshape(2, -1).T

    old_dps, old_prec = mpmath.mp.dps, mpmath.mp.prec
    try:
        mpmath.mp.dps = 400

        assert_func_equal(sc.beta,
                          lambda a, b: float(mpmath.beta(a, b)),
                          ab,
                          vectorized=False,
                          rtol=1e-10,
                          ignore_inf_sign=True)

        assert_func_equal(
            sc.betaln,
            lambda a, b: float(mpmath.log(abs(mpmath.beta(a, b)))),
            ab,
            vectorized=False,
            rtol=1e-10)
    finally:
        mpmath.mp.dps, mpmath.mp.prec = old_dps, old_prec


#------------------------------------------------------------------------------
# Machinery for systematic tests
#------------------------------------------------------------------------------ 

def eta(lam):
    """Function from DLMF 8.12.1 shifted to be centered at 0."""
    if lam > 0:
        return mp.sqrt(2*(lam - mp.log(lam + 1)))
    elif lam < 0:
        return -mp.sqrt(2*(lam - mp.log(lam + 1)))
    else:
        return 0 

def test_log(self):
        with mp.workdps(30):
            logcoeffs = mp.taylor(lambda x: mp.log(1 + x), 0, 10)
            expcoeffs = mp.taylor(lambda x: mp.exp(x) - 1, 0, 10)
            invlogcoeffs = lagrange_inversion(logcoeffs)
            mp_assert_allclose(invlogcoeffs, expcoeffs) 

def test_beta():
    np.random.seed(1234)

    b = np.r_[np.logspace(-200, 200, 4),
              np.logspace(-10, 10, 4),
              np.logspace(-1, 1, 4),
              np.arange(-10, 11, 1),
              np.arange(-10, 11, 1) + 0.5,
              -1, -2.3, -3, -100.3, -10003.4]
    a = b

    ab = np.array(np.broadcast_arrays(a[:,None], b[None,:])).reshape(2, -1).T

    old_dps, old_prec = mpmath.mp.dps, mpmath.mp.prec
    try:
        mpmath.mp.dps = 400

        assert_func_equal(sc.beta,
                          lambda a, b: float(mpmath.beta(a, b)),
                          ab,
                          vectorized=False,
                          rtol=1e-10,
                          ignore_inf_sign=True)

        assert_func_equal(
            sc.betaln,
            lambda a, b: float(mpmath.log(abs(mpmath.beta(a, b)))),
            ab,
            vectorized=False,
            rtol=1e-10)
    finally:
        mpmath.mp.dps, mpmath.mp.prec = old_dps, old_prec


# ------------------------------------------------------------------------------
# loggamma
# ------------------------------------------------------------------------------ 

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_log1p_complex(self):
        assert_mpmath_equal(sc.log1p,
                            lambda x: mpmath.log(x+1),
                            [ComplexArg()], dps=60) 

def test_log_ndtr(self):
        assert_mpmath_equal(sc.log_ndtr,
                            exception_to_nan(lambda z: mpmath.log(mpmath.ncdf(z))),
                            [Arg()], n=600, dps=300) 

def compute_log_moment(q, sigma, steps, lmbd, verify=False, verbose=False):
  """Compute the log moment of Gaussian mechanism for given parameters.
  Args:
    q: the sampling ratio.
    sigma: the noise sigma.
    steps: the number of steps.
    lmbd: the moment order.
    verify: if False, only compute the symbolic version. If True, computes
      both symbolic and numerical solutions and verifies the results match.
    verbose: if True, print out debug information.
  Returns:
    the log moment with type np.float64, could be np.inf.
  """
  moment = compute_a(sigma, q, lmbd, verbose=verbose)
  if verify:
    mp.dps = 50
    moment_a_mp = compute_a_mp(sigma, q, lmbd, verbose=verbose)
    moment_b_mp = compute_b_mp(sigma, q, lmbd, verbose=verbose)
    np.testing.assert_allclose(moment, moment_a_mp, rtol=1e-10)
    if not np.isinf(moment_a_mp):
      # The following test fails for (1, np.inf)!
      np.testing.assert_array_less(moment_b_mp, moment_a_mp)
  if np.isinf(moment):
    return np.inf
  else:
    return np.log(moment) * steps 

def compute_b(sigma, q, lmbd, verbose=False):
  mu0, _, mu = distributions(sigma, q)

  b_lambda_fn = lambda z: mu0(z) * np.power(cropped_ratio(mu0(z), mu(z)), lmbd)
  b_lambda = integral_inf(b_lambda_fn)
  m = sigma ** 2 * (np.log((2. - q) / (1. - q)) + 1. / (2 * sigma ** 2))

  b_fn = lambda z: (np.power(mu0(z) / mu(z), lmbd) -
                    np.power(mu(-z) / mu0(z), lmbd))
  if verbose:
    print "M =", m
    print "f(-M) = {} f(M) = {}".format(b_fn(-m), b_fn(m))
    assert b_fn(-m) < 0 and b_fn(m) < 0

  b_lambda_int1_fn = lambda z: (mu0(z) *
                                np.power(cropped_ratio(mu0(z), mu(z)), lmbd))
  b_lambda_int2_fn = lambda z: (mu0(z) *
                                np.power(cropped_ratio(mu(z), mu0(z)), lmbd))
  b_int1 = integral_bounded(b_lambda_int1_fn, -m, m)
  b_int2 = integral_bounded(b_lambda_int2_fn, -m, m)

  a_lambda_m1 = compute_a(sigma, q, lmbd - 1)
  b_bound = a_lambda_m1 + b_int1 - b_int2

  if verbose:
    print "B: by numerical integration", b_lambda
    print "B must be no more than     ", b_bound
  print b_lambda, b_bound
  return _to_np_float64(b_lambda)


###########################
# MULTIPRECISION ROUTINES #
########################### 

def test_lgam1p(self):
        def param_filter(x):
            # Filter the poles
            return np.where((np.floor(x) == x) & (x <= 0), False, True)

        def mp_lgam1p(z):
            # The real part of loggamma is log(|gamma(z)|)
            return mpmath.loggamma(1 + z).real

        assert_mpmath_equal(_lgam1p,
                            mp_lgam1p,
                            [Arg()], rtol=1e-13, dps=100,
                            param_filter=param_filter) 

def test_gammaln(self):
        # The real part of loggamma is log(|gamma(z)|).
        def f(z):
            return mpmath.loggamma(z).real

        assert_mpmath_equal(sc.gammaln, exception_to_nan(f), [Arg()]) 

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 Brkic_2011_1(Re, eD):
    r'''Calculates Darcy friction factor using the method in Brkic
    (2011) [2]_ as shown in [1]_.

    .. math::
        f_d = [-2\log(10^{-0.4343\beta} + \frac{\epsilon}{3.71D})]^{-2}

    .. math::
        \beta = \ln \frac{Re}{1.816\ln\left(\frac{1.1Re}{\ln(1+1.1Re)}\right)}

    Parameters
    ----------
    Re : float
        Reynolds number, [-]
    eD : float
        Relative roughness, [-]

    Returns
    -------
    fd : float
        Darcy friction factor [-]

    Notes
    -----
    No range of validity specified for this equation.

    Examples
    --------
    >>> Brkic_2011_1(1E5, 1E-4)
    0.01812455874141297

    References
    ----------
    .. [1] Winning, H. and T. Coole. "Explicit Friction Factor Accuracy and
       Computational Efficiency for Turbulent Flow in Pipes." Flow, Turbulence
       and Combustion 90, no. 1 (January 1, 2013): 1-27.
       doi:10.1007/s10494-012-9419-7
    .. [2] 	Brkic, Dejan."Review of Explicit Approximations to the Colebrook
       Relation for Flow Friction." Journal of Petroleum Science and
       Engineering 77, no. 1 (April 2011): 34-48.
       doi:10.1016/j.petrol.2011.02.006.
    '''
    beta = log(Re/(1.816*log(1.1*Re/log(1+1.1*Re))))
    return (-2*log10(10**(-0.4343*beta)+eD/3.71))**-2