"""Verify the Levy flight-length distribution and its inverse-transform sampler.

Source: viswanathan1999levy, equation (1): P(l) ~ l^(-mu) for 1 < mu <= 3, where
l is the flight length.  The hifuku Levy jump draws l by inverse transform, then
applies round(l) SPR moves.

This script confirms:
  1. the normalized pdf p(l) = (mu - 1) * l^(-mu) on [1, inf) integrates to 1,
  2. the CDF is F(l) = 1 - l^(-(mu - 1)),
  3. inverting U = F(l) gives the sampler l = V^(-1/(mu - 1)) with V = 1 - U,
  4. the survival function of the sampled length has the tail exponent mu - 1,
     so the density tail is l^(-mu),
  5. a numeric check of the empirical tail exponent for mu = 2.
"""

import numpy as np
import sympy as sp

l, mu, U, V, x = sp.symbols("l mu U V x", positive=True)

# --- 1. Normalization of p(l) = (mu - 1) l^(-mu) on [1, inf) ---
# The improper integral converges only for mu > 1, so SymPy returns a Piecewise.
pdf = (mu - 1) * l ** (-mu)
total = sp.integrate(pdf, (l, 1, sp.oo))
if isinstance(total, sp.Piecewise):
    total = total.args[0].expr      # the mu > 1 branch
total = sp.simplify(total)
assert total == 1, f"pdf does not normalize: {total}"

# --- 2. CDF F(l) = integral of pdf from 1 to l ---
cdf = sp.integrate(pdf, (l, 1, x))
if isinstance(cdf, sp.Piecewise):
    cdf = cdf.args[0].expr          # the mu > 1 (or mu < 1) branch
cdf = sp.simplify(cdf)
expected_cdf = 1 - x ** (-(mu - 1))
assert sp.simplify(cdf - expected_cdf) == 0, f"CDF mismatch: {cdf}"

# --- 3. Inverse transform: solve U = F(l) for l ---
sol = sp.solve(sp.Eq(U, 1 - l ** (-(mu - 1))), l)
inv = sol[0]
# With V = 1 - U, the sampler is l = V^(-1/(mu - 1)).
sampler = V ** (-1 / (mu - 1))
check = sp.simplify(inv.subs(U, 1 - V) - sampler)
assert check == 0, f"inverse-transform sampler mismatch: {inv}"

# --- 4. Survival of the sampled length: P(L > x) = x^(-(mu - 1)) ---
# L = V^(-1/(mu-1)) > x  <=>  V < x^(-(mu-1)), and V ~ Uniform(0,1), so
# P(L > x) = x^(-(mu-1)); the density is its negative derivative ~ x^(-mu).
survival = x ** (-(mu - 1))
density_tail = sp.simplify(-sp.diff(survival, x))
assert sp.simplify(density_tail - (mu - 1) * x ** (-mu)) == 0, "tail exponent wrong"

# --- 5. Numeric empirical tail exponent for mu = 2 ---
rng = np.random.default_rng(0)
mu_val = 2.0
v = rng.random(2_000_000)
samples = v ** (-1.0 / (mu_val - 1.0))
# Fit the complementary CDF slope on a log-log grid in the tail.
grid = np.logspace(0.5, 2.5, 40)
surv = np.array([(samples > g).mean() for g in grid])
mask = surv > 0
slope = np.polyfit(np.log(grid[mask]), np.log(surv[mask]), 1)[0]
# Survival exponent should be -(mu - 1) = -1 for mu = 2.
assert abs(slope - (-(mu_val - 1.0))) < 0.05, f"empirical slope {slope}"

print("All Levy flight-length checks passed.")
print(f"  pdf normalizes to {total}")
print(f"  CDF F(l) = {expected_cdf.subs(x, l)}")
print(f"  sampler l = V**(-1/(mu-1))")
print(f"  empirical survival slope for mu=2: {slope:.4f} (expected -1.0)")