#!/usr/bin/env python3
"""DDIA Ch. 9 - The Trouble with Distributed Systems

There is no safe failure-detection timeout.

A healthy node's round-trip time (RTT) is heavy-tailed: fast almost always,
but load / GC pauses / network queueing occasionally produce a spike. A
reader reasons "median RTT is 50 ms, so a 300 ms timeout (6x headroom) never
fires on a live node." This script measures two curves that CANNOT both be
small:

  - FALSE-POSITIVE RATE  = P(healthy RTT > T): a live node whose heartbeat
    exceeded the timeout, so it is wrongly declared dead.
  - DETECTION DELAY      ~= T: for a genuinely dead node you must wait about
    the whole timeout before you are entitled to declare it dead.

The RTT is modelled as a lognormal, so P(RTT > T) is available EXACTLY from
the standard-normal CDF (statistics.NormalDist) -- no sampling error. A
seeded Monte Carlo draw is run alongside purely to confirm the analytic
numbers and to report the distribution's percentiles.

Deterministic: analytic values are exact; the Monte Carlo uses random.Random(42).
"""

import math
import random
import statistics

# ---- RTT distribution -------------------------------------------------------
# Lognormal: log(RTT) ~ Normal(MU, SIGMA).  Median = exp(MU).
# Calibrated so the median is 50 ms and the tail is fat.
MEDIAN_MS = 50.0
MU = math.log(MEDIAN_MS)     # 3.912 -> median exp(MU) = 50 ms exactly
SIGMA = 1.0                  # fat tail: p99 ~ half a second, extreme tail multi-second

STD_NORMAL = statistics.NormalDist(0.0, 1.0)


def fp_rate(timeout_ms):
    """Exact P(healthy RTT > timeout) = 1 - Phi((ln T - MU)/SIGMA).

    This is the irreducible false-positive floor of a heartbeat timeout T:
    the fraction of healthy heartbeats whose RTT exceeds T by chance.
    """
    z = (math.log(timeout_ms) - MU) / SIGMA
    return 1.0 - STD_NORMAL.cdf(z)


def rtt_percentile(p):
    """Exact p-quantile of the lognormal RTT, in ms."""
    z = STD_NORMAL.inv_cdf(p)
    return math.exp(MU + SIGMA * z)


TIMEOUTS_MS = [100, 200, 300, 500, 1000, 2000, 5000]
SAFE_TIMEOUT_MS = 300          # the reader's "6x median = surely safe" choice
SAFE_HEADROOM = SAFE_TIMEOUT_MS / MEDIAN_MS


def line(width=64):
    print("-" * width)


def step1_distribution():
    print("=" * 64)
    print("Step 1 - the healthy node's RTT distribution (calibration)")
    line()
    print("RTT ~ lognormal(mu=%.3f, sigma=%.1f)  =>  median = %.0f ms"
          % (MU, SIGMA, MEDIAN_MS))
    print()
    for label, p in (("p50 ", 0.50), ("p90 ", 0.90), ("p99 ", 0.99),
                     ("p99.9", 0.999), ("p99.99", 0.9999)):
        print("    %-7s = %8.1f ms" % (label, rtt_percentile(p)))
    print()
    print("    A healthy node: median %.0f ms, but the p99.9 is a full %.0fx"
          % (MEDIAN_MS, rtt_percentile(0.999) / MEDIAN_MS))
    print("    the median and the extreme tail runs into seconds.")
    print()


def step2_safe_timeout():
    print("=" * 64)
    print("Step 2 - the 'obviously safe' timeout: 6x the median = 300 ms")
    line()
    fp = fp_rate(SAFE_TIMEOUT_MS)
    print("Timeout T = %d ms  (%.0fx the %.0f ms median headroom)"
          % (SAFE_TIMEOUT_MS, SAFE_HEADROOM, MEDIAN_MS))
    print()
    print("    FALSE-POSITIVE RATE  P(healthy RTT > %d ms) = %.4f = %.2f%%"
          % (SAFE_TIMEOUT_MS, fp, fp * 100))
    print("    i.e. 1 healthy heartbeat in %.0f is wrongly declared DEAD."
          % (1.0 / fp,))
    print()
    print("    NOT zero. 6x the median leaves a fat slice of the tail")
    print("    above the timeout, and every heartbeat that lands there is")
    print("    a live node mistaken for a corpse.")
    print()


def step3_sweep():
    print("=" * 64)
    print("Step 3 - sweep the timeout: false positives vs detection delay")
    line()
    print("%8s   %18s   %20s" % ("T (ms)", "false-positive rate",
                                  "detection delay (dead)"))
    print("%8s   %18s   %20s" % ("------", "-------------------",
                                  "--------------------"))
    for t in TIMEOUTS_MS:
        fp = fp_rate(t)
        mark = "   <- 6x median" if t == SAFE_TIMEOUT_MS else ""
        print("%8d   %17.4f%%   %17d ms%s"
              % (t, fp * 100, t, mark))
    print()
    print("    Read the two columns against each other: every ms you add to")
    print("    T to shrink the false-positive rate is a ms you must WAIT")
    print("    before detecting a genuinely dead node. They move in opposite")
    print("    directions. There is no row where both are small.")
    print()


def step4_chasing_zero():
    print("=" * 64)
    print("Step 4 - how big must T be to drive false positives to ~zero?")
    line()
    for target in (0.01, 0.001, 0.0001):
        # invert: find T with P(RTT > T) = target  =>  T = exp(MU + SIGMA*z_{1-target})
        z = STD_NORMAL.inv_cdf(1.0 - target)
        t = math.exp(MU + SIGMA * z)
        print("    to reach FP = %6.3f%%   need T = %8.0f ms  (%5.0fx median,"
              " wait %.2f s per dead node)"
              % (target * 100, t, t / MEDIAN_MS, t / 1000.0))
    print()
    print("    FP only reaches EXACTLY 0 at T = infinity -- a node you never")
    print("    declare dead. The tail has no upper bound, so no finite timeout")
    print("    buys a zero false-positive rate.")
    print()


def main():
    # Seeded Monte Carlo, used only to CONFIRM the analytic numbers.
    rng = random.Random(42)
    N = 1_000_000
    samples = sorted(rng.lognormvariate(MU, SIGMA) for _ in range(N))

    step1_distribution()
    step2_safe_timeout()
    step3_sweep()
    step4_chasing_zero()

    # Cross-check: empirical FP at the safe timeout vs the exact value.
    print("=" * 64)
    print("Monte Carlo cross-check (n=%s, seed=42)" % f"{N:,}")
    line()
    empirical_fp = sum(1 for s in samples if s > SAFE_TIMEOUT_MS) / N
    print("    empirical P(RTT > %d ms) = %.4f%%   (analytic %.4f%%)"
          % (SAFE_TIMEOUT_MS, empirical_fp * 100, fp_rate(SAFE_TIMEOUT_MS) * 100))

    def emp_pct(p):
        return samples[min(int(p * N), N - 1)]
    print("    empirical p50 = %6.1f ms   (analytic %6.1f ms)"
          % (emp_pct(0.50), rtt_percentile(0.50)))
    print("    empirical p99 = %6.1f ms   (analytic %6.1f ms)"
          % (emp_pct(0.99), rtt_percentile(0.99)))
    print("    empirical max = %6.1f ms   over %s draws"
          % (samples[-1], f"{N:,}"))
    print("=" * 64)


if __name__ == "__main__":
    main()