Chapter 42 — Exponential Growth

# --- Visualization: Linear vs polynomial vs exponential growth ---
import numpy as np
import matplotlib.pyplot as plt

plt.style.use('seaborn-v0_8-whitegrid')

n = np.linspace(0, 20, 300)

fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# --- Panel 1: Linear scale (exponential dominates quickly) ---
ax = axes[0]
ax.plot(n, n,         label='Linear: n',       color='green')
ax.plot(n, n**2,      label='Quadratic: n²',   color='blue')
ax.plot(n, n**3,      label='Cubic: n³',       color='orange')
ax.plot(n, 2**n,      label='Exponential: 2ⁿ', color='red', linewidth=2)
ax.set_ylim(0, 500)
ax.set_title('Growth Comparison (Linear Scale)')
ax.set_xlabel('n')
ax.set_ylabel('f(n)')
ax.legend()

# --- Panel 2: Log scale (straightens exponential to a line) ---
ax = axes[1]
ax.plot(n, n,         label='Linear: n',       color='green')
ax.plot(n, n**2,      label='Quadratic: n²',   color='blue')
ax.plot(n, n**3,      label='Cubic: n³',       color='orange')
ax.plot(n, 2**n,      label='Exponential: 2ⁿ', color='red', linewidth=2)
ax.set_yscale('log')
ax.set_title('Growth Comparison (Log Scale)')
ax.set_xlabel('n')
ax.set_ylabel('f(n)  [log scale]')
ax.legend()

plt.tight_layout()
plt.show()

# --- Visualization: Continuous exponential and decay; doubling time ---
import numpy as np
import matplotlib.pyplot as plt

plt.style.use('seaborn-v0_8-whitegrid')

t = np.linspace(0, 5, 500)
P0 = 1.0          # initial value
k_growth = 0.8    # growth rate
k_decay  = -0.8   # decay rate

P_growth = P0 * np.exp(k_growth * t)
P_decay  = P0 * np.exp(k_decay  * t)
doubling_time = np.log(2) / k_growth  # exact formula

fig, axes = plt.subplots(1, 2, figsize=(13, 5))

ax = axes[0]
ax.plot(t, P_growth, color='red',   label=f'Growth: e^({k_growth}t)')
ax.plot(t, P_decay,  color='blue',  label=f'Decay:  e^({k_decay}t)')
ax.axhline(2 * P0, color='gray', linestyle='--', linewidth=0.8)
ax.axvline(doubling_time, color='red', linestyle='--', linewidth=0.8,
           label=f'Doubling time ≈ {doubling_time:.2f}')
ax.set_title('Continuous Exponential Growth and Decay')
ax.set_xlabel('Time t')
ax.set_ylabel('P(t)')
ax.legend()

# Compound interest: annual compounding vs continuous
ax = axes[1]
years = np.linspace(0, 30, 300)
r = 0.07  # 7% annual rate
for n_per_year, label in [(1, 'Annual (n=1)'), (12, 'Monthly (n=12)'),
                           (365, 'Daily (n=365)'), (np.inf, 'Continuous')]:
    if n_per_year == np.inf:
        vals = np.exp(r * years)
    else:
        vals = (1 + r / n_per_year) ** (n_per_year * years)
    ax.plot(years, vals, label=label)
ax.set_title(f'Compound Interest at {int(r*100)}% (r={r})')
ax.set_xlabel('Years')
ax.set_ylabel('Growth of $1')
ax.legend()

plt.tight_layout()
plt.show()

# --- Verify the limit definition of e ---
import numpy as np

print("Approximating e via (1 + 1/n)^n")
print("-" * 45)
for n in [1, 10, 100, 1_000, 1_000_000, 1_000_000_000]:
    approx = (1 + 1/n) ** n
    print(f"n = {n:>12,}:  (1+1/n)^n = {approx:.10f}")
print(f"\nnp.e           =             {np.e:.10f}")
print("\nNote: the convergence is slow (only ~log10(n) decimal digits per order of magnitude)")

# --- Implementation: Discrete and continuous exponential models ---
import numpy as np

def discrete_exponential(P0, r, n_steps):
    """
    Simulate discrete exponential growth/decay.

    P(n) = P0 * r^n

    Args:
        P0:      initial value
        r:       growth factor per step (r > 1 grows, 0 < r < 1 decays)
        n_steps: number of steps (non-negative integer)

    Returns:
        Array of values P(0), P(1), ..., P(n_steps)
    """
    steps = np.arange(0, n_steps + 1)
    return P0 * r ** steps


def continuous_exponential(P0, k, t_max, n_points=500):
    """
    Evaluate continuous exponential model P(t) = P0 * e^(kt).

    Args:
        P0:       initial value at t=0
        k:        continuous rate (positive=growth, negative=decay)
        t_max:    upper time limit
        n_points: resolution

    Returns:
        t array, P(t) array
    """
    t = np.linspace(0, t_max, n_points)
    P = P0 * np.exp(k * t)
    return t, P


def doubling_time(r=None, k=None):
    """
    Compute doubling time from either a discrete factor r or continuous rate k.
    Uses: T_double = ln(2) / k  (continuous)  or  log(2)/log(r) (discrete)

    Args:
        r: discrete growth factor (> 1)
        k: continuous growth rate (> 0)

    Returns:
        Doubling time (in same units as the growth rate)
    """
    if k is not None:
        return np.log(2) / k
    elif r is not None:
        return np.log(2) / np.log(r)
    else:
        raise ValueError("Provide either r or k")


# Demonstrations
print("Discrete model: bacterial colony doubling every hour")
colony = discrete_exponential(P0=100, r=2, n_steps=10)
for i, v in enumerate(colony):
    print(f"  Hour {i:2d}: {v:>10,.0f} bacteria")

print(f"\nDoubling time check: {doubling_time(r=2):.4f} steps  (expected 1.0)")

print("\nRule of 72 validation (7% growth):")
exact = doubling_time(k=0.07)
rule_of_72 = 72 / 7
print(f"  Exact doubling time: {exact:.4f} years")
print(f"  Rule of 72 estimate: {rule_of_72:.4f} years")
print(f"  Error: {abs(exact - rule_of_72)/exact * 100:.2f}%")

# --- Experiment 1: The chessboard problem — intuition calibration ---
# Hypothesis: 2^63 is incomprehensibly large.
# Let's compute it and compare to real quantities.
# Try changing: GRAIN_MASS_GRAMS to compare to other large quantities

import numpy as np

GRAIN_MASS_GRAMS = 0.03     # mass of one grain of rice in grams  <-- modify this
EARTH_MASS_KG    = 5.97e24  # mass of Earth in kg

total_grains = sum(2**i for i in range(64))
total_mass_g = total_grains * GRAIN_MASS_GRAMS
total_mass_kg = total_mass_g / 1000

print(f"Total grains on chessboard: {total_grains:.3e}")
print(f"= 2^63 + 2^62 + ... + 1 = 2^64 - 1")
print(f"\nTotal mass: {total_mass_g:.3e} grams = {total_mass_kg:.3e} kg")
print(f"Earth's mass:             {EARTH_MASS_KG:.3e} kg")
print(f"Ratio (rice / Earth):     {total_mass_kg / EARTH_MASS_KG:.3e}")
print(f"\nRice mass is {total_mass_kg / EARTH_MASS_KG:.1e} × Earth's mass.")

# --- Experiment 2: Algorithm complexity — when does O(2^n) become infeasible? ---
# Hypothesis: An O(2^n) algorithm becomes unusable far before an O(n^3) one.
# Assume 10^9 operations per second.
# Try changing: OPS_PER_SEC to see how hardware improvements help (or don't)

import numpy as np

OPS_PER_SEC = 1e9        # <-- modify this (try 1e12 for a supercomputer)
SECONDS_IN_YEAR = 3.15e7

print(f"Assumed speed: {OPS_PER_SEC:.0e} operations/second")
print(f"{'n':>5} {'O(n^2)':>15} {'O(n^3)':>15} {'O(2^n)':>20}")
print("-" * 60)

def format_time(seconds):
    if seconds < 1e-6:   return f"{seconds*1e9:.2f} ns"
    if seconds < 1e-3:   return f"{seconds*1e6:.2f} μs"
    if seconds < 1:      return f"{seconds*1e3:.2f} ms"
    if seconds < 60:     return f"{seconds:.2f} s"
    if seconds < 3600:   return f"{seconds/60:.2f} min"
    if seconds < 86400:  return f"{seconds/3600:.2f} hr"
    if seconds < SECONDS_IN_YEAR: return f"{seconds/86400:.2f} days"
    years = seconds / SECONDS_IN_YEAR
    if years < 1e6:      return f"{years:.2e} years"
    return f"{years:.2e} years"

for n in [10, 20, 30, 40, 50, 60, 100]:
    t_n2 = n**2    / OPS_PER_SEC
    t_n3 = n**3    / OPS_PER_SEC
    t_2n = 2.0**n  / OPS_PER_SEC
    print(f"{n:>5} {format_time(t_n2):>15} {format_time(t_n3):>15} {format_time(t_2n):>20}")

# --- Experiment 3: Compounding frequency and convergence to e^r ---
# Hypothesis: More frequent compounding converges to continuous compounding,
# but with diminishing returns.
# Try changing: ANNUAL_RATE and YEARS

import numpy as np
import matplotlib.pyplot as plt

plt.style.use('seaborn-v0_8-whitegrid')

ANNUAL_RATE = 0.10   # <-- modify this (try 0.05, 0.20)
YEARS       = 10     # <-- modify this
P0          = 1000   # initial principal

n_values = np.logspace(0, 6, 200)  # 1 to 10^6 compounding periods per year
final_values = P0 * (1 + ANNUAL_RATE / n_values) ** (n_values * YEARS)
continuous   = P0 * np.exp(ANNUAL_RATE * YEARS)

fig, ax = plt.subplots(figsize=(9, 5))
ax.semilogx(n_values, final_values, color='steelblue', label='Discrete compounding')
ax.axhline(continuous, color='red', linestyle='--', label=f'Continuous: ${continuous:.2f}')

for n, label in [(1, 'Annual'), (12, 'Monthly'), (365, 'Daily')]:
    val = P0 * (1 + ANNUAL_RATE / n) ** (n * YEARS)
    ax.scatter([n], [val], zorder=5)
    ax.annotate(f'{label}\n${val:.2f}', (n, val), textcoords='offset points',
                xytext=(10, -15), fontsize=8)

ax.set_title(f'Compounding ${P0} at {ANNUAL_RATE*100:.0f}%/yr for {YEARS} years')
ax.set_xlabel('Compounding periods per year (log scale)')
ax.set_ylabel(f'Final value ($)')
ax.legend()
plt.tight_layout()
plt.show()

# --- Mini Project: Early-phase epidemic as exponential growth ---
#
# Problem:
#   In the early phase of an epidemic, when the infected fraction is small,
#   the number of infected individuals grows approximately exponentially:
#   I(t) ≈ I₀ × e^(β t)
#   where β is the transmission rate.
#
#   Your task: fit an exponential model to simulated epidemic data and
#   estimate the growth rate and implied doubling time.
#
# Note: ch079 (Project: Epidemic Simulation) will build the full SIR model.
#       Here we focus on the exponential approximation only.

import numpy as np
import matplotlib.pyplot as plt

plt.style.use('seaborn-v0_8-whitegrid')
np.random.seed(42)

# --- Simulated "observed" data (exponential + noise) ---
TRUE_RATE = 0.25          # true growth rate β
I0        = 10            # initial infected
NOISE_STD = 5             # observation noise

t_obs = np.arange(0, 21)  # days 0 to 20
I_true = I0 * np.exp(TRUE_RATE * t_obs)
I_obs  = I_true + np.random.normal(0, NOISE_STD * I_true / I_true.max(), len(t_obs))
I_obs  = np.maximum(I_obs, 1)  # can't have negative infected

# TODO 1: Log-transform the data to linearize the exponential
#         ln(I) = ln(I0) + β × t  →  linear in t
log_I_obs = np.log(I_obs)

# TODO 2: Fit a line to (t_obs, log_I_obs) using least squares
#         β_hat = slope, ln(I0_hat) = intercept
coeffs = np.polyfit(t_obs, log_I_obs, 1)  # degree-1 polynomial
beta_hat = coeffs[0]
I0_hat   = np.exp(coeffs[1])

# TODO 3: Reconstruct the fitted exponential and compare to true
t_fine   = np.linspace(0, 20, 300)
I_fitted = I0_hat * np.exp(beta_hat * t_fine)
I_truth  = I0 * np.exp(TRUE_RATE * t_fine)

# TODO 4: Compute doubling time from estimated rate
T_double_hat  = np.log(2) / beta_hat
T_double_true = np.log(2) / TRUE_RATE

# --- Plot ---
fig, axes = plt.subplots(1, 2, figsize=(13, 5))

ax = axes[0]
ax.scatter(t_obs, I_obs, label='Observed data', zorder=5, color='gray')
ax.plot(t_fine, I_truth,  label=f'True (β={TRUE_RATE})',      color='green', linestyle='--')
ax.plot(t_fine, I_fitted, label=f'Fitted (β̂={beta_hat:.3f})', color='red')
ax.set_title('Epidemic Early Phase: Exponential Fit')
ax.set_xlabel('Day')
ax.set_ylabel('Infected count')
ax.legend()

ax = axes[1]
ax.scatter(t_obs, log_I_obs, label='log(Observed)', zorder=5, color='gray')
ax.plot(t_fine, np.log(I_truth),  label='log(True)',  color='green', linestyle='--')
ax.plot(t_fine, np.log(I_fitted), label='log(Fitted)', color='red')
ax.set_title('Log-transformed: Exponential → Linear')
ax.set_xlabel('Day')
ax.set_ylabel('log(Infected)')
ax.legend()

plt.tight_layout()
plt.show()

print(f"\nEstimation Results:")
print(f"  True  β: {TRUE_RATE:.4f}   Estimated β̂: {beta_hat:.4f}")
print(f"  True  I₀: {I0:.1f}       Estimated I₀: {I0_hat:.1f}")
print(f"  True  doubling time: {T_double_true:.2f} days")
print(f"  Estimated doubling time: {T_double_hat:.2f} days")