Python math.log1p() Examples

def test_math_subclass(self):
        """verify subtypes of float/long work w/ math functions"""
        import math
        class myfloat(float): pass
        class mylong(long): pass

        mf = myfloat(1)
        ml = mylong(1)

        for x in math.log, math.log10, math.log1p, math.asinh, math.acosh, math.atanh, math.factorial, math.trunc, math.isinf:
            try:
                resf = x(mf)
            except ValueError:
                resf = None
            try:
                resl = x(ml)
            except ValueError:
                resl = None
            self.assertEqual(resf, resl) 

def _isadensalt_m(self, dens_kgpm3):
        """Returns the geopotential alt (m) at which ISA has given density (kg/m^3)"""
        dense_it = np.nditer([dens_kgpm3, None])
        for d_kgpm3, alt_m in dense_it:
            if d_kgpm3 > self._dLevel1:
                # Troposphere
                alt_m[...] = (d_kgpm3 ** (1 / self._L1) - self._I1) / self._J1
            elif d_kgpm3 > self._dLevel2:
                # Lower stratopshere
                alt_m[...] = (1 / self._N2) * math.log1p(dens_kgpm3 / self._M2 - 1)
            elif d_kgpm3 > self._dLevel3:
                # Upper stratopshere
                alt_m[...] = (d_kgpm3 ** (1 / self._L3) - self._I3) / self._J3
            elif dens_kgpm3 > self._dLevel4:
                # Between 32 and 47 km
                alt_m[...] = (d_kgpm3 ** (1 / self._L4) - self._I4) / self._J4
            else:
                # Between 47km and 51km
                alt_m[...] = (1 / self._N5) * math.log1p(d_kgpm3 / self._M5 - 1)
            # Adjust for the temperature offset
        return dense_it.operands[1] 

def include_revision(revision_num, skip_factor=1.1):
  """Decide whether to include a revision.

  If the number of revisions is large, we exclude some revisions to avoid
  a quadratic blowup in runtime, since the article is likely also large.

  We make the ratio between consecutive included revision numbers
  appproximately equal to "factor".

  Args:
    revision_num: an integer
    skip_factor: a floating point number >= 1.0

  Returns:
    a boolean
  """
  if skip_factor <= 1.0:
    return True
  return (int(math.log1p(revision_num) / math.log(skip_factor)) != int(
      math.log(revision_num + 2.0) / math.log(skip_factor))) 

def _isapressalt_m(self, pressure_pa):
        """Returns the geopotential alt (m) at which ISA has the given pressure"""
        press_it = np.nditer([pressure_pa, None])
        for press_pa, alt_m in press_it:
            if press_pa > self._pLevel1:
                # Troposphere
                alt_m[...] = (press_pa ** (1 / self._E1) - self._C1) / self._D1
            elif press_pa > self._pLevel2:
                # Lower stratopshere
                alt_m[...] = (1 / self._G2) * math.log1p(press_pa / self._F2 - 1)
            elif press_pa > self._pLevel3:
                # Upper stratopshere
                alt_m[...] = (press_pa ** (1 / self._E3) - self._C3) / self._D3
            elif press_pa > self._pLevel4:
                # Between 32 and 47 km
                alt_m[...] = (press_pa ** (1 / self._E4) - self._C4) / self._D4
            else:
                # Between 47km and 51km
                alt_m[...] = (1 / self._G5) * math.log1p(press_pa / self._F5 - 1)
        return press_it.operands[1] 

def include_revision(revision_num, skip_factor=1.1):
  """Decide whether to include a revision.

  If the number of revisions is large, we exclude some revisions to avoid
  a quadratic blowup in runtime, since the article is likely also large.

  We make the ratio between consecutive included revision numbers
  appproximately equal to "factor".

  Args:
    revision_num: an integer
    skip_factor: a floating point number >= 1.0

  Returns:
    a boolean
  """
  if skip_factor <= 1.0:
    return True
  return (int(math.log1p(revision_num) / math.log(skip_factor)) != int(
      math.log(revision_num + 2.0) / math.log(skip_factor))) 

def test_math_subclass(self):
        """verify subtypes of float/long work w/ math functions"""
        import math
        class myfloat(float): pass
        class mylong(long): pass

        mf = myfloat(1)
        ml = mylong(1)

        for x in math.log, math.log10, math.log1p, math.asinh, math.acosh, math.atanh, math.factorial, math.trunc, math.isinf:
            try:
                resf = x(mf)
            except ValueError:
                resf = None
            try:
                resl = x(ml)
            except ValueError:
                resl = None
            self.assertEqual(resf, resl) 

def __rsub__(self, other):
        if other < 0:  # Result will be negative
            return other + -float(self)

        log_other = self._log(other)
        if log_other > self.log_value:
            log_l, log_s = log_other, self.log_value
        elif log_other < self.log_value:  # Result will be negative
            return other + -float(self)
        else:  # Must be equal, so return 0
            return Probability(float("-inf"), log_value=True)

        if log_s == float("-inf"):  # Just return largest value
            return Probability(log_l, log_value=True)

        exp_diff = exp(log_s - log_l)
        if exp_diff == 1:  # Diff too small, so result is effectively zero
            return Probability(float("-inf"), log_value=True)

        return Probability(log_l + log1p(-exp_diff),
                           log_value=True) 

def __sub__(self, other):
        if other < 0:
            return self + -other

        log_other = self._log(other)
        if self.log_value > log_other:
            log_l, log_s = self.log_value, log_other
        elif self.log_value < log_other:  # Result will be negative
            return float(self) - other
        else:  # Must be equal, so return 0
            return Probability(float("-inf"), log_value=True)

        if log_s == float("-inf"):  # Just return largest value
            return Probability(log_l, log_value=True)

        exp_diff = exp(log_s - log_l)
        if exp_diff == 1:  # Diff too small, so result is effectively zero
            return Probability(float("-inf"), log_value=True)

        return Probability(log_l + log1p(-exp_diff),
                           log_value=True) 

def __add__(self, other):
        if other < 0:
            return self - -other

        log_other = self._log(other)
        if self.log_value > log_other:
            log_l, log_s = self.log_value, log_other
        elif self.log_value < log_other:
            log_l, log_s = log_other, self.log_value
        else:  # Must be equal, so just double value
            return self * 2

        if log_s == float("-inf"):  # Just return largest value
            return Probability(log_l, log_value=True)

        return Probability(log_l + log1p(exp(log_s - log_l)),
                           log_value=True) 

def _log1mexp(x):
  """Numerically stable computation of log(1-exp(x))."""
  if x < -1:
    return math.log1p(-math.exp(x))
  elif x < 0:
    return math.log(-math.expm1(x))
  elif x == 0:
    return -np.inf
  else:
    raise ValueError("Argument must be non-positive.") 

def compute_logq_gaussian(counts, sigma):
  """Returns an upper bound on ln Pr[outcome != argmax] for GNMax.

  Implementation of Proposition 7.

  Args:
    counts: A numpy array of scores.
    sigma: The standard deviation of the Gaussian noise in the GNMax mechanism.

  Returns:
    logq: Natural log of the probability that outcome is different from argmax.
  """
  n = len(counts)
  variance = sigma**2
  idx_max = np.argmax(counts)
  counts_normalized = counts[idx_max] - counts
  counts_rest = counts_normalized[np.arange(n) != idx_max]  # exclude one index
  # Upper bound q via a union bound rather than a more precise calculation.
  logq = _logaddexp(
      scipy.stats.norm.logsf(counts_rest, scale=math.sqrt(2 * variance)))

  # A sketch of a more accurate estimate, which is currently disabled for two
  # reasons:
  # 1. Numerical instability;
  # 2. Not covered by smooth sensitivity analysis.
  # covariance = variance * (np.ones((n - 1, n - 1)) + np.identity(n - 1))
  # logq = np.log1p(-statsmodels.sandbox.distributions.extras.mvnormcdf(
  #     counts_rest, np.zeros(n - 1), covariance, maxpts=1e4))

  return min(logq, math.log(1 - (1 / n))) 

def log1p(node): return merge([node], math.log1p) 

def rdp_pure_eps(logq, pure_eps, orders):
  """Computes the RDP value given logq and pure privacy eps.

  Implementation of https://arxiv.org/abs/1610.05755, Theorem 3.

  The bound used is the min of three terms. The first term is from
  https://arxiv.org/pdf/1605.02065.pdf.
  The second term is based on the fact that when event has probability (1-q) for
  q close to zero, q can only change by exp(eps), which corresponds to a
  much smaller multiplicative change in (1-q)
  The third term comes directly from the privacy guarantee.

  Args:
    logq: Natural logarithm of the probability of a non-optimal outcome.
    pure_eps: eps parameter for DP
    orders: array_like list of moments to compute.

  Returns:
    Array of upper bounds on rdp (a scalar if orders is a scalar).
  """
  orders_vec = np.atleast_1d(orders)
  q = math.exp(logq)
  log_t = np.full_like(orders_vec, np.inf)
  if q <= 1 / (math.exp(pure_eps) + 1):
    logt_one = math.log1p(-q) + (
        math.log1p(-q) - _log1mexp(pure_eps + logq)) * (
            orders_vec - 1)
    logt_two = logq + pure_eps * (orders_vec - 1)
    log_t = np.logaddexp(logt_one, logt_two)

  ret = np.minimum(
      np.minimum(0.5 * pure_eps * pure_eps * orders_vec,
                 log_t / (orders_vec - 1)), pure_eps)
  if np.isscalar(orders):
    return np.asscalar(ret)
  else:
    return ret 

def compute_logq_gaussian(counts, sigma):
  """Returns an upper bound on ln Pr[outcome != argmax] for GNMax.

  Implementation of Proposition 7.

  Args:
    counts: A numpy array of scores.
    sigma: The standard deviation of the Gaussian noise in the GNMax mechanism.

  Returns:
    logq: Natural log of the probability that outcome is different from argmax.
  """
  n = len(counts)
  variance = sigma**2
  idx_max = np.argmax(counts)
  counts_normalized = counts[idx_max] - counts
  counts_rest = counts_normalized[np.arange(n) != idx_max]  # exclude one index
  # Upper bound q via a union bound rather than a more precise calculation.
  logq = _logaddexp(
      scipy.stats.norm.logsf(counts_rest, scale=math.sqrt(2 * variance)))

  # A sketch of a more accurate estimate, which is currently disabled for two
  # reasons:
  # 1. Numerical instability;
  # 2. Not covered by smooth sensitivity analysis.
  # covariance = variance * (np.ones((n - 1, n - 1)) + np.identity(n - 1))
  # logq = np.log1p(-statsmodels.sandbox.distributions.extras.mvnormcdf(
  #     counts_rest, np.zeros(n - 1), covariance, maxpts=1e4))

  return min(logq, math.log(1 - (1 / n))) 

def ZettaYear(probability):
	trials = 2*365*10 ** 21
	return 1-exp(trials * log1p(-probability)) 

def _log_add(logx, logy):
  """Add two numbers in the log space."""
  a, b = min(logx, logy), max(logx, logy)
  if a == -np.inf:  # adding 0
    return b
  # Use exp(a) + exp(b) = (exp(a - b) + 1) * exp(b)
  return math.log1p(math.exp(a - b)) + b  # log1p(x) = log(x + 1) 

def testLog1p(self):
        self.assertRaises(TypeError, math.log1p)
        self.ftest('log1p(1/e -1)', math.log1p(1/math.e-1), -1)
        self.ftest('log1p(0)', math.log1p(0), 0)
        self.ftest('log1p(e-1)', math.log1p(math.e-1), 1)
        self.ftest('log1p(1)', math.log1p(1), math.log(2))
        self.assertEqual(math.log1p(INF), INF)
        self.assertRaises(ValueError, math.log1p, NINF)
        self.assertTrue(math.isnan(math.log1p(NAN)))
        n= 2**90
        self.assertAlmostEqual(math.log1p(n), 62.383246250395075)
        self.assertAlmostEqual(math.log1p(n), math.log1p(float(n))) 

def testLog1p(self):
        self.assertRaises(TypeError, math.log1p)
        for n in [2, 2**90, 2**300]:
            self.assertAlmostEqual(math.log1p(n), math.log1p(float(n)))
        self.assertRaises(ValueError, math.log1p, -1)
        self.assertEqual(math.log1p(INF), INF) 

def _log1mexp(x):
  """Numerically stable computation of log(1-exp(x))."""
  if x < -1:
    return math.log1p(-math.exp(x))
  elif x < 0:
    return math.log(-math.expm1(x))
  elif x == 0:
    return -np.inf
  else:
    raise ValueError("Argument must be non-positive.") 

def _log1mexp(x):
  """Numerically stable computation of log(1-exp(x))."""
  if x < -1:
    return math.log1p(-math.exp(x))
  elif x < 0:
    return math.log(-math.expm1(x))
  elif x == 0:
    return -np.inf
  else:
    raise ValueError("Argument must be non-positive.") 

def testLog1p(self):
        self.assertRaises(TypeError, math.log1p)
        n= 2**90
        self.assertAlmostEqual(math.log1p(n), math.log1p(float(n))) 

def MegaYear(probability):
	trials = 2*365*10 ** 6
	return 1-exp(trials * log1p(-probability))



##We evaluate the function pcalc(10,5,3,4)
##print pcalc(10,5,3,4)
##as a test vector
##The answer would be [binom(3,4)*binom(2,6)+(binom(4,4)*binom(1,6)]/binom(10,5)
##[4*15+1*6]/252 = 66/252
##print float(66)/252

##quorum size for ChainLocks 

def log1p(x):
    """log(1 + x) accurate for small x (missing from python 2.5.2)"""

    if sys.version_info > (2, 6):
      return math.log1p(x)

    y = 1 + x
    z = y - 1
    # Here's the explanation for this magic: y = 1 + z, exactly, and z
    # approx x, thus log(y)/z (which is nearly constant near z = 0) returns
    # a good approximation to the true log(1 + x)/x.  The multiplication x *
    # (log(y)/z) introduces little additional error.
    return x if z == 0 else x * math.log(y) / z 

def log_add(log_x, log_y):
    """Given log x and log y, returns log(x + y)."""
    # Swap variables so log_y is larger.
    if log_x > log_y:
        log_x, log_y = log_y, log_x

    # Use the log(1 + e^p) trick to compute this efficiently
    # If the difference is large enough, this is effectively log y.
    delta = log_y - log_x
    return math.log1p(math.exp(delta)) + log_x if delta <= 50.0 else log_y 

def geometric(p):
    if p <= 0 or p > 1:
        raise ValueError('p must be in the interval (0.0, 1.0]')
    if p == 1:
        return 1
    return int(math.log1p(-random.random()) / math.log1p(-p)) + 1 

def log1m_exp(val):
    """Numerically stable implementation of `log(1 - exp(val))`."""
    if val >= 0.:
        return nan
    elif val > LOG_2:
        return log(-expm1(val))
    else:
        return log1p(-exp(val)) 

def log1p_exp(val):
    """Numerically stable implementation of `log(1 + exp(val))`."""
    if val > 0.:
        return val + log1p(exp(-val))
    else:
        return log1p(exp(val)) 

def testLog1p(self):
        self.assertRaises(TypeError, math.log1p)
        self.ftest('log1p(1/e -1)', math.log1p(1/math.e-1), -1)
        self.ftest('log1p(0)', math.log1p(0), 0)
        self.ftest('log1p(e-1)', math.log1p(math.e-1), 1)
        self.ftest('log1p(1)', math.log1p(1), math.log(2))
        self.assertEqual(math.log1p(INF), INF)
        self.assertRaises(ValueError, math.log1p, NINF)
        self.assertTrue(math.isnan(math.log1p(NAN)))
        n= 2**90
        self.assertAlmostEqual(math.log1p(n), 62.383246250395075)
        self.assertAlmostEqual(math.log1p(n), math.log1p(float(n))) 

def update(self, item: str, date: datetime.date) -> None:
        """Add 'like' for item.

        Args:
            item: An item in the list that is being ranked.
            date: The date on which the item has been liked.
        """
        score = self.get(item)
        time = date.toordinal()
        higher = max(score, time * self.rate)
        lower = min(score, time * self.rate)
        self.scores[item] = higher + math.log1p(math.exp(lower - higher)) 

def _log_add(logx, logy):
  """Add two numbers in the log space."""
  a, b = min(logx, logy), max(logx, logy)
  if a == -np.inf:  # adding 0
    return b
  # Use exp(a) + exp(b) = (exp(a - b) + 1) * exp(b)
  return math.log1p(math.exp(a - b)) + b  # log1p(x) = log(x + 1)