import numpy as np

# -------------------------
# CF Delta-nstep Validation Locked inputs. Gaps N1 significance and N2 robustmness
# -------------------------
delta_win = np.array([1,2,3,4,8,12,20,32,54,90,146], dtype=int)   # Δn_win actual CF delta-nSteps for post-crossing
D_fibo = np.array([1,2,3,5,8,13,21,34,55,89,144], dtype=int) # fibo test
n_cf = np.array([1,2,4,7,11,19,31,51,83,137,227,373],dtype=float)
m = np.arange(1, len(n_cf)+1, dtype=float)

# -------------------------
# Helper: score a window of gaps against targets
# -------------------------
def score_window(delta, targets):
    d = np.min(np.abs(delta[:,None] - targets[None,:]), axis=1)
    S = int(np.sum(d))
    H = int(np.sum(d == 0))
    return S, H, d

# Observed scores
obsS_fibo, obsH_fibo, _ = score_window(delta_win, D_fibo)

# -------------------------
# N1: growth-matched log-like null
# OLS on ln(n_m) = A + B m + eps, compute residual sigma
# -------------------------
lnn = np.log(n_cf)
X = np.vstack([np.ones_like(m), m]).T
A, B = np.linalg.lstsq(X, lnn, rcond=None)[0]
resid = lnn - (A + B*m)
sigma = np.sqrt(np.sum(resid**2)/(len(m)-2))  # df = n-2

# -------------------------
# Null generator N1
# -------------------------
def generate_N1(N, rng, A, B, sigma):
    K = len(n_cf)
    eps = rng.normal(0, sigma, size=(N, K))
    base = A + B*np.arange(1, K+1)
    ntilde = np.exp(base[None,:] + eps)
    n = np.floor(ntilde).astype(int)
    n[n < 1] = 1

    # strict monotonicity "minimal repair"
    for j in range(1, n.shape[1]):
        n[:, j] = np.maximum(n[:, j], n[:, j-1] + 1)

    gaps = n[:, 1:] - n[:, :-1]                 # length 8
    return gaps

# -------------------------
# Null generator N2
# -------------------------
def generate_N2(N, rng, n_cf_int):
    eta = rng.integers(-1, 2, size=(N, len(n_cf_int)))      # uniform in {-1,0,+1}
    n = n_cf_int[None, :] + eta
    n[n < 1] = 1

    # strict monotonicity "minimal repair"
    for j in range(1, n.shape[1]):
        n[:, j] = np.maximum(n[:, j], n[:, j-1] + 1)

    gaps = n[:, 1:] - n[:, :-1]
    return gaps.astype(int)


# -------------------------
# N3: Smooth-growth log-correlated null
# -------------------------
def generate_N3(N, rng, A, B, sigma, n_cf, rho=0.8):
    """
    Generate growth-matched landmark sequences with
    correlated logarithmic residuals.

    rho = adjacent residual correlation
    rho=0 gives approximately N1 behavior.
    rho near 1 produces smoother deviations.

    Returns an array of N simulated gap sequences.
    """
    K = len(n_cf)

    # Stationary AR(1) correlated residuals
    eps = np.empty((N, K), dtype=float)

    eps[:, 0] = rng.normal(0, sigma, size=N)

    innovation_sd = sigma * np.sqrt(1 - rho**2)

    for j in range(1, K):
        eps[:, j] = (
            rho * eps[:, j-1]
            + rng.normal(0, innovation_sd, size=N)
        )

    # Exponential-growth baseline
    base = A + B * np.arange(1, K+1)

    ntilde = np.exp(base[None, :] + eps)

    # Integer landmarks
    n = np.floor(ntilde).astype(int)
    n[n < 1] = 1

    # Enforce strict monotonicity
    for j in range(1, K):
        n[:, j] = np.maximum(
            n[:, j],
            n[:, j-1] + 1
        )

    # Eleven gaps from twelve landmarks
    gaps = np.diff(n, axis=1)

    return gaps


# -------------------------
# N4: Geometric Growth and Phase Null
# -------------------------
def generate_N4(N, rng, A, B, K,
                A_halfwidth=0.25,
                B_halfwidth=0.05):
    """
    Generate deterministic geometric landmark
    sequences with randomized growth rate and
    starting scale.

    Each realization has its own A' and B'.
    No independent landmark noise is added.

    Returns:
        gaps: (N, K-1) integer array
        A_samples, B_samples
    """

    A_samples = A + rng.uniform(
        -A_halfwidth, A_halfwidth, size=N
    )

    B_samples = B + rng.uniform(
        -B_halfwidth, B_halfwidth, size=N
    )

    m_values = np.arange(1, K+1)

    log_landmarks = (
        A_samples[:, None]
        + B_samples[:, None] * m_values[None, :]
    )

    n = np.floor(np.exp(log_landmarks)).astype(int)
    n[n < 1] = 1

    # Strictly increasing integer landmarks
    for j in range(1, K):
        n[:, j] = np.maximum(
            n[:, j], n[:, j-1] + 1
        )

    gaps = np.diff(n, axis=1)

    return gaps, A_samples, B_samples


# -------------------------
# Gap Score Generator
# -------------------------
def scores_for_gaps(gaps, targets):
    d = np.min(np.abs(gaps[:,:,None] - targets[None,None,:]), axis=2)  # (N,8)
    S = np.sum(d, axis=1)
    H = np.sum(d == 0, axis=1)
    return S, H

# -------------------------
# Tail Probability Score Generator
# -------------------------
def tail_probs(S, H, obsS, obsH):
    pS = np.mean(S <= obsS)
    pH = np.mean(H >= obsH)
    pJ = np.mean((S <= obsS) & (H >= obsH))
    return pS, pH, pJ


# -------------------------
# Run Monte Carlo
# -------------------------
rng = np.random.default_rng(12345)
N = 200000

# N1
gaps1 = generate_N1(N, rng, A, B, sigma)
S1_fibo, H1_fibo = scores_for_gaps(gaps1, D_fibo)

tp1_fibo  = tail_probs(S1_fibo,  H1_fibo,  obsS_fibo,  obsH_fibo)

all4_1 = np.mean((S1_fibo  <= obsS_fibo)  & (H1_fibo  >= obsH_fibo))

# N2
gaps2 = generate_N2(N, rng, n_cf.astype(int))
S2_fibo, H2_fibo = scores_for_gaps(gaps2, D_fibo)

tp2_fibo  = tail_probs(S2_fibo,  H2_fibo,  obsS_fibo,  obsH_fibo)

all4_2 = np.mean((S2_fibo  <= obsS_fibo)  & (H2_fibo  >= obsH_fibo))

# -------------------------
# N3: Smooth-growth Monte Carlo
# -------------------------

rho_N3 = 0.8  # Prespecified correlation

gaps3 = generate_N3(
    N, rng, A, B, sigma, n_cf,
    rho=rho_N3
)

S3_fibo, H3_fibo = scores_for_gaps(
    gaps3, D_fibo
)

tp3_fibo = tail_probs(
    S3_fibo, H3_fibo,
    obsS_fibo, obsH_fibo
)

all4_3 = np.mean(
    (S3_fibo <= obsS_fibo) &
    (H3_fibo >= obsH_fibo)
)


# -------------------------
# N4: Monte Carlo evaluation
# -------------------------

gaps4, A4, B4 = generate_N4(
    N=N,
    rng=rng,
    A=A,
    B=B,
    K=len(n_cf),
    A_halfwidth=0.25,
    B_halfwidth=0.05
)

S4_fibo, H4_fibo = scores_for_gaps(
    gaps4, D_fibo
)

tp4_fibo = tail_probs(
    S4_fibo, H4_fibo,
    obsS_fibo, obsH_fibo
)

print("A,B,sigma =", A, B, sigma)
print("Observed fibo  S,H =", obsS_fibo, obsH_fibo)

print("N1 fibo  (pS,pH,pJ) =", tp1_fibo)
print("N1 joint =", all4_1)

print("N2 fibo (pS,pH,pJ) =", tp2_fibo)
print("N2 joint =", all4_2)


print("N3 rho =", rho_N3)
print("N3 fibo (pS,pH,pJ) =", tp3_fibo)
print("N3 joint =", all4_3)

print("N3 S median =", np.median(S3_fibo))
print("N3 S 95th percentile =",
      np.percentile(S3_fibo, 95))

print("N4 fibo (pS,pH,pJ) =", tp4_fibo)

print("N4 S median =", np.median(S4_fibo))
print("N4 S 5th percentile =",
      np.percentile(S4_fibo, 5))
print("N4 S 95th percentile =",
      np.percentile(S4_fibo, 95))

print("N4 growth-factor range =",
      np.exp(B-0.05), np.exp(B+0.05))

# Best-fitting sampled geometric alternative
best_idx = np.argmin(S4_fibo)

print("N4 best S,H =",
      S4_fibo[best_idx],
      H4_fibo[best_idx])

print("N4 best A,B =",
      A4[best_idx],
      B4[best_idx])

print("N4 best gaps =",
      gaps4[best_idx])