Blame for modules/script_02_irregular_diffusion.py

"""Script 2: irregular Laola response with wave train, damping and diffusion proxy.

Standalone point:
- Builds on Script 1 conceptually, but remains executable on its own.
- Individual response delays create cumulative timing error.
- Repeated start impulses show how a packet train becomes noisy and damped.
- Transport-diffusion is a diagnostic coarse-grained reading, not an input law.
"""

from __future__ import annotations

from dataclasses import dataclass
from typing import Tuple

import gradio as gr
import numpy as np

try:
    from modules.common import (
        BASE,
        activity_from_schedule,
        count_visible_packets,
        exact_regular_schedule,
        launch_offsets,
        make_heatmap,
        make_snapshot_comparison,
        make_timing_error_plot,
        time_grid,
        timing_error,
        transport_diffusion_residuals,
    )
except ModuleNotFoundError:  # direct execution as python modules/script_XX.py
    from common import (
        BASE,
        activity_from_schedule,
        count_visible_packets,
        exact_regular_schedule,
        launch_offsets,
        make_heatmap,
        make_snapshot_comparison,
        make_timing_error_plot,
        time_grid,
        timing_error,
        transport_diffusion_residuals,
    )


@dataclass(frozen=True)
class IrregularParams:
    n_people: int = BASE.n_people
    tau_mean: float = BASE.tau
    tau_jitter: float = 0.035
    activation_time: float = BASE.activation_time
    refractory_time: float = BASE.refractory_time
    seat_spacing: float = BASE.seat_spacing
    launch_count: int = BASE.launch_count
    launch_interval: float = BASE.launch_interval
    damping_per_node: float = 0.003
    dt: float = BASE.dt
    t_max: float = BASE.t_max
    snapshot_time: float = BASE.snapshot_time
    seed: int = 7


def activation_schedule(params: IrregularParams) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
    """Return activation times, edge delays and raw jitter ξ.

    Critical detail: if tau_jitter == 0 we construct t_i = i τ directly.
    We do not use cumsum in that zero-noise case, because recursive float addition
    of 0.2 creates ~1e-14 diagnostic artefacts although the model has no jitter.
    """
    n = int(params.n_people)
    tau = float(params.tau_mean)
    sigma = max(float(params.tau_jitter), 0.0)

    if n <= 0:
        return np.zeros(0), np.zeros(0), np.zeros(0)

    if sigma == 0.0:
        activation_times = exact_regular_schedule(n, tau)
        edge_delays = np.full(max(n - 1, 0), tau, dtype=float)
        raw_jitter = np.zeros(max(n - 1, 0), dtype=float)
        return activation_times, edge_delays, raw_jitter

    rng = np.random.default_rng(int(params.seed))
    raw_jitter = rng.normal(0.0, sigma, size=max(n - 1, 0))
    edge_delays = np.clip(tau + raw_jitter, 0.02, None)
    activation_times = np.zeros(n, dtype=float)
    if n > 1:
        activation_times[1:] = np.cumsum(edge_delays)
    return activation_times, edge_delays, raw_jitter


def simulate(params: IrregularParams) -> Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
    n = int(params.n_people)
    times = time_grid(params.t_max, params.dt)
    positions = params.seat_spacing * np.arange(n, dtype=float)
    activation_times, edge_delays, raw_jitter = activation_schedule(params)
    launches = launch_offsets(params.launch_count, params.launch_interval)

    base = activity_from_schedule(times, activation_times, launches, params.activation_time)
    damping = np.exp(-float(params.damping_per_node) * np.arange(n, dtype=float))
    heights = base * damping[None, :]

    reference_activation = exact_regular_schedule(n, params.tau_mean)
    reference = activity_from_schedule(times, reference_activation, launches, params.activation_time)
    return times, positions, activation_times, edge_delays, raw_jitter, launches, heights, reference


def estimate_diffusion(params: IrregularParams, edge_delays: np.ndarray) -> float:
    """Heuristic effective diffusion coefficient.

    If Var(t_i) ≈ i σ_τ² and x ≈ i a, then the spatial variance at time t≈iτ
    scales like a² σ_τ² t / τ³. Matching Var(x) ≈ 2D t yields
    D_eff ≈ a² σ_τ² / (2τ³).

    This is a diagnostic proxy, not a theorem. For zero jitter it is exactly zero.
    """
    if edge_delays.size == 0:
        return 0.0
    tau = max(float(np.mean(edge_delays)), 1e-12)
    sigma2 = float(np.var(edge_delays))
    if sigma2 <= 1e-30:
        return 0.0
    return (params.seat_spacing**2 * sigma2) / (2.0 * tau**3)


def run_model(
    n_people: int,
    tau_mean: float,
    tau_jitter: float,
    activation_time: float,
    refractory_time: float,
    seat_spacing: float,
    launch_count: int,
    launch_interval: float,
    damping_per_node: float,
    dt: float,
    t_max: float,
    snapshot_time: float,
    seed: int,
):
    params = IrregularParams(
        n_people=int(n_people),
        tau_mean=float(tau_mean),
        tau_jitter=float(tau_jitter),
        activation_time=float(activation_time),
        refractory_time=float(refractory_time),
        seat_spacing=float(seat_spacing),
        launch_count=int(launch_count),
        launch_interval=float(launch_interval),
        damping_per_node=float(damping_per_node),
        dt=float(dt),
        t_max=float(t_max),
        snapshot_time=float(snapshot_time),
        seed=int(seed),
    )

    times, positions, activation_times, edge_delays, raw_jitter, launches, heights, reference = simulate(params)
    reference_activation = exact_regular_schedule(params.n_people, params.tau_mean)
    errors = timing_error(activation_times, reference_activation)

    if float(params.tau_jitter) == 0.0:
        tau_emp = float(params.tau_mean)
        sigma_emp = 0.0
    else:
        tau_emp = float(np.mean(edge_delays)) if edge_delays.size else float(params.tau_mean)
        sigma_emp = float(np.std(edge_delays)) if edge_delays.size else 0.0
    c_emp = params.seat_spacing / max(tau_emp, 1e-12)
    d_eff = estimate_diffusion(params, edge_delays)
    r_transport, r_diffusion = transport_diffusion_residuals(times, positions, heights, c_emp, d_eff)

    zero_irregularity = (float(params.tau_jitter) == 0.0) and (float(params.damping_per_node) == 0.0)
    heat_title = (
        "Script 2: identisch zu Script 1, weil ξ=0 und α=0"
        if zero_irregularity
        else "Script 2: irregulärer/dämpfender Laola-Wellenzug"
    )
    heatmap = make_heatmap(times, positions, heights, heat_title)
    snapshot = make_snapshot_comparison(
        times,
        positions,
        heights,
        reference,
        params.snapshot_time,
        "Script 2 simuliert",
    )
    schedule = make_timing_error_plot(
        positions,
        errors,
        "Script 2: kumulativer Timing-Fehler relativ zur idealen Script-1-Referenz",
    )

    peak_left = float(np.max(heights[:, 0])) if heights.size else float("nan")
    peak_right = float(np.max(heights[:, -1])) if heights.size else float("nan")
    attenuation_ratio = peak_right / max(peak_left, 1e-12)
    snapshot_idx = int(np.argmin(np.abs(times - params.snapshot_time)))
    visible_packets = count_visible_packets(heights[snapshot_idx, :])
    active_nodes = int(np.sum(heights[snapshot_idx, :] > 0.05))
    launch_spacing_nodes = params.launch_interval / max(tau_emp, 1e-12)
    refractory_warning = params.launch_interval < params.refractory_time
    warning_text = "Ja — Wellen starten enger als die gesetzte Refraktärzeit." if refractory_warning else "Nein."
    max_abs_error = float(np.max(np.abs(errors))) if errors.size else 0.0
    max_abs_delta_snapshot = float(np.max(np.abs(heights[snapshot_idx, :] - reference[snapshot_idx, :]))) if heights.size else 0.0

    zero_note = "Ja — Script 2 fällt exakt auf die Script-1-Referenz zurück." if zero_irregularity else "Nein."

    interpretation = f"""
## Script 2 — Interpretation

**Diskrete Regel mit irregulärem Wellenzug**

```text
t_0 = 0
t_{{i+1}} = t_i + τ_i
τ_i = τ + ξ_i
t_i^(m) = t_i + m · T_launch
h_i(t) = exp(-α i) · max_m p(t - t_i^(m))
```

**v4-Korrektur:** Für `std(ξ)=0` wird die Aktivierungszeit direkt als `t_i=i·τ` konstruiert.
Damit verschwindet die frühere scheinbare kumulative Abweichung von ca. `1e-14 s`.
Das war nur Floating-Point-Rundung aus rekursiver `cumsum`-Addition, kein physikalischer Timing-Fehler.

| Größe | Wert |
|---|---:|
| Knotenanzahl N | {params.n_people} |
| Startimpulse launch_count | {params.launch_count} |
| Startabstand T_launch | {params.launch_interval:.3f} s |
| Startabstand in Knoten T_launch/mean(τ_i) | {launch_spacing_nodes:.2f} Knoten |
| sichtbare Pakete im Snapshot | {visible_packets} |
| aktive Knoten im Snapshot | {active_nodes} |
| mittlere Reaktionszeit mean(τ_i) | {tau_emp:.4f} s |
| Timing-Streuung std(τ_i) | {sigma_emp:.4e} s |
| max. kumulativer Timing-Fehler | {max_abs_error:.4e} s |
| effektive Frontgeschwindigkeit c≈a/mean(τ_i) | {c_emp:.4f} Knotenabstände/s |
| Dämpfung pro Knoten α | {params.damping_per_node:.5f} |
| Peak rechts / Peak links | {attenuation_ratio:.4f} |
| max. Snapshot-Abweichung gegen Script 1 | {max_abs_delta_snapshot:.4e} |
| heuristisches D_eff≈a²Var(ξ)/(2τ³) | {d_eff:.6e} |
| Refraktärwarnung T_launch < T_ref | {warning_text} |
| ξ=0 und α=0 ⇒ exakt Script 1? | {zero_note} |
| normierter Transport-Residualindikator | {r_transport:.4f} |
| normierter Transport-Diffusions-Residualindikator | {r_diffusion:.4f} |

**Lesart:** Lokale Unregelmäßigkeiten erzeugen kumulative Timing-Fehler. Nach coarse-graining
ist die saubere Transportgleichung nicht mehr die beste Sprache; ein Transport-Diffusions-Term
ist die naheliegende effektive Näherung:

```text
∂t h + c ∂x h ≈ D ∂x² h
ω(k) ≈ c k - i D k²
```

**Unsicherheit:** `D_eff` ist hier eine kontrollierte heuristische Diagnose, kein bewiesener Grenzwertsatz.
Für CNNA wäre später entscheidend, ob die effektiven Diffusions-/Dämpfungsterme aus Schur-/DtN-
Randantworten des wachsenden Netzwerks abgeleitet werden können.
"""
    return heatmap, snapshot, schedule, interpretation


def build_tab() -> None:
    gr.Markdown(
        r"""
## Script 2 — Unregelmäßige Menschen → Dämpfung/Diffusion

Dieses Modul ist eigenständig, übernimmt aber die Idee von Script 1. Neu ist:
individuelle Reaktionszeit `τ_i = τ + ξ_i`, kumulative Frontunschärfe und optionaler Response-Verlust.

Ab v4 hat Script 2 dieselben Vergleichsplots wie Script 1. Bei `std(ξ)=0` und `α=0` muss es exakt
auf Script 1 zurückfallen.
"""
    )
    with gr.Row():
        with gr.Column(scale=1):
            n_people = gr.Slider(20, 240, value=BASE.n_people, step=1, label="N Menschen / Knoten")
            tau_mean = gr.Slider(0.05, 0.80, value=BASE.tau, step=0.01, label="mittlere τ [s]")
            tau_jitter = gr.Slider(0.00, 0.20, value=0.035, step=0.005, label="std(ξ) Timing-Streuung [s]")
            activation_time = gr.Slider(0.20, 4.00, value=BASE.activation_time, step=0.05, label="T_act Aktivitätsdauer [s]")
            refractory_time = gr.Slider(0.20, 8.00, value=BASE.refractory_time, step=0.05, label="T_ref Refraktärzeit [s]")
            seat_spacing = gr.Slider(0.20, 3.00, value=BASE.seat_spacing, step=0.05, label="a Knotenabstand")
            launch_count = gr.Slider(1, 12, value=BASE.launch_count, step=1, label="launch_count Startimpulse")
            launch_interval = gr.Slider(0.50, 12.00, value=BASE.launch_interval, step=0.10, label="T_launch Abstand der Startimpulse [s]")
            damping_per_node = gr.Slider(0.00, 0.04, value=0.003, step=0.001, label="α Dämpfung pro Knoten")
            dt = gr.Slider(0.01, 0.10, value=BASE.dt, step=0.01, label="dt Zeitauflösung [s]")
            t_max = gr.Slider(4.00, 80.00, value=BASE.t_max, step=0.50, label="t_max [s]")
            snapshot_time = gr.Slider(0.00, 80.00, value=BASE.snapshot_time, step=0.10, label="Momentaufnahme t [s]")
            seed = gr.Number(value=7, precision=0, label="Seed für reproduzierbare ξᵢ")
            run = gr.Button("Simulation ausführen", variant="primary")
        with gr.Column(scale=2):
            heatmap = gr.Plot(label="Raum-Zeit-Aktivität")
            snapshot = gr.Plot(label="Momentaufnahme gegen ideale Referenz")
            schedule = gr.Plot(label="kumulativer Timing-Fehler")
            interpretation = gr.Markdown()

    inputs = [
        n_people,
        tau_mean,
        tau_jitter,
        activation_time,
        refractory_time,
        seat_spacing,
        launch_count,
        launch_interval,
        damping_per_node,
        dt,
        t_max,
        snapshot_time,
        seed,
    ]
    outputs = [heatmap, snapshot, schedule, interpretation]
    run.click(run_model, inputs=inputs, outputs=outputs)


if __name__ == "__main__":
    params = IrregularParams()
    t, x, activations, delays, raw_jitter, launches, h, ref = simulate(params)
    print(t.shape, x.shape, activations.shape, delays.shape, raw_jitter.shape, launches.shape, h.shape, ref.shape)