"""Script 3: refractory response gate for the Laola analogue. Standalone point: - Builds on Script 1/2, but adds the first local state dynamic. - A node does not respond automatically when a signal arrives. - It responds only if its local live state has recovered from previous activation. - Therefore repeated waves can pass, weaken, or fail. Didactic status: - Still not a CNNA Schur/DtN model. - First transition from kinematic schedules to stateful response propagation. """ from __future__ import annotations from dataclasses import dataclass from typing import Literal, Tuple import gradio as gr import matplotlib.pyplot as plt 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, response_profile, ) 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, response_profile, ) RefractoryMode = Literal["binary", "partial_recovery"] FailurePolicy = Literal["block", "attenuate"] # Script 3 intentionally uses a stress default so the new refractory gate is visible. # The Script-1/2 baseline remains available by setting T_launch back to BASE.launch_interval. SCRIPT3_DEMO_LAUNCH_INTERVAL = 2.00 SCRIPT3_DEMO_SNAPSHOT_TIME = 18.00 @dataclass(frozen=True) class RefractoryParams: n_people: int = BASE.n_people tau_mean: float = BASE.tau tau_jitter: float = 0.0 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 = SCRIPT3_DEMO_LAUNCH_INTERVAL damping_per_node: float = 0.0 dt: float = BASE.dt t_max: float = BASE.t_max snapshot_time: float = SCRIPT3_DEMO_SNAPSHOT_TIME refractory_mode: RefractoryMode = "binary" failure_policy: FailurePolicy = "block" wrap_ring: bool = False seed: int = 7 def edge_delay_schedule(params: RefractoryParams) -> Tuple[np.ndarray, np.ndarray]: """Return edge delays τ_i and raw jitter ξ_i. Edge i is the delay from node i to node i+1 in the open-chain case. For ring mode the last edge N-1 goes from node N-1 back to node 0. Zero-jitter is constructed exactly as a constant array to avoid turning binary floating point drift into a false model diagnostic. """ n = max(0, int(params.n_people)) edge_count = n if bool(params.wrap_ring) else max(n - 1, 0) tau = float(params.tau_mean) sigma = max(float(params.tau_jitter), 0.0) if edge_count == 0: return np.zeros(0, dtype=float), np.zeros(0, dtype=float) if sigma == 0.0: return np.full(edge_count, tau, dtype=float), np.zeros(edge_count, dtype=float) rng = np.random.default_rng(int(params.seed)) raw_jitter = rng.normal(0.0, sigma, size=edge_count) edge_delays = np.clip(tau + raw_jitter, 0.02, None) return edge_delays, raw_jitter def _recovery_factor(delta_since_last: float, params: RefractoryParams) -> float: """Local live-state gate for a new response attempt. delta_since_last measures the time between the current arrival and the node's previous activation. The node is still active for T_act, then refractory for T_ref. In binary mode it accepts only after the full recovery gap. In partial mode it recovers linearly after the active period. """ if not np.isfinite(delta_since_last): return 1.0 t_act = max(float(params.activation_time), 0.0) t_ref = max(float(params.refractory_time), 1e-12) full_gap = t_act + t_ref if str(params.refractory_mode) == "binary": return 1.0 if delta_since_last >= full_gap else 0.0 # partial_recovery: active phase gives zero response; refractory phase recovers linearly. return float(np.clip((delta_since_last - t_act) / t_ref, 0.0, 1.0)) def propagate_events( params: RefractoryParams, ) -> Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray]: """Propagate launch signals through local refractory gates. Returns ------- event_times: Shape (launch_count, N). Accepted activation time or NaN. response_strength: Shape (launch_count, N). Local response amplitude. 1 = full response, 0 = blocked or unreached. status: Shape (launch_count, N). 1 accepted, 0 blocked at this node, -1 unreached. edge_delays: Edge delay schedule. raw_jitter: Raw ξ values used to build edge_delays. """ n = max(0, int(params.n_people)) launches = launch_offsets(params.launch_count, params.launch_interval) m_count = len(launches) event_times = np.full((m_count, n), np.nan, dtype=float) response_strength = np.zeros((m_count, n), dtype=float) status = np.full((m_count, n), -1.0, dtype=float) last_activation = np.full(n, -np.inf, dtype=float) edge_delays, raw_jitter = edge_delay_schedule(params) for m, launch_time in enumerate(launches): arrival_time = float(launch_time) node = 0 visited_edges = 0 max_edges = n if bool(params.wrap_ring) else max(n - 1, 0) while 0 <= node < n: delta = arrival_time - last_activation[node] recovery = _recovery_factor(delta, params) accepted = False local_strength = 0.0 if str(params.failure_policy) == "block": if recovery > 0.0: accepted = True local_strength = 1.0 if str(params.refractory_mode) == "binary" else recovery else: # attenuate accepted = True local_strength = recovery if not accepted or local_strength <= 0.0: status[m, node] = 0.0 break status[m, node] = 1.0 event_times[m, node] = arrival_time response_strength[m, node] = local_strength last_activation[node] = arrival_time if not bool(params.wrap_ring): if node >= n - 1: break edge_idx = node arrival_time = arrival_time + float(edge_delays[edge_idx]) node = node + 1 else: if visited_edges >= max_edges - 1: break edge_idx = node arrival_time = arrival_time + float(edge_delays[edge_idx]) node = (node + 1) % n visited_edges += 1 return event_times, response_strength, status, edge_delays, raw_jitter def activity_from_events( times: np.ndarray, event_times: np.ndarray, response_strength: np.ndarray, activation_time: float, damping_per_node: float, ) -> np.ndarray: """Build H[t, i] from accepted stateful events.""" if event_times.ndim != 2: raise ValueError("event_times must have shape (launch_count, n_people)") launch_count, n = event_times.shape heights = np.zeros((len(times), n), dtype=float) node_damping = np.exp(-max(float(damping_per_node), 0.0) * np.arange(n, dtype=float)) for m in range(launch_count): for i in range(n): t0 = event_times[m, i] amp = float(response_strength[m, i]) * float(node_damping[i]) if not np.isfinite(t0) or amp <= 0.0: continue local = times - float(t0) heights[:, i] = np.maximum(heights[:, i], amp * response_profile(local, activation_time)) return heights def first_activation_times(event_times: np.ndarray) -> np.ndarray: """First accepted activation per node; NaN if no accepted event reached it.""" if event_times.size == 0: return np.zeros(0, dtype=float) with np.errstate(all="ignore"): first = np.nanmin(event_times, axis=0) first[~np.isfinite(first)] = np.nan return first def first_wave_reference_from_edges(n: int, edge_delays: np.ndarray, tau_mean: float) -> np.ndarray: """Reference schedule for the first launch using the same edge delays. With zero jitter this is exactly i·τ. With jitter it matches Script 2's cumulative edge-delay schedule, so refractory failures are not confused with timing noise. """ if n <= 0: return np.zeros(0, dtype=float) if edge_delays.size == 0: return np.zeros(n, dtype=float) if np.all(edge_delays == float(tau_mean)): return exact_regular_schedule(n, float(tau_mean)) out = np.zeros(n, dtype=float) if n > 1: out[1:] = np.cumsum(edge_delays[: n - 1]) return out def make_response_deficit_plot(positions: np.ndarray, status: np.ndarray): """Plot missing accepted responses relative to the ideal Script-1 launch count. Script 1 assumes every launch reaches every node. Script 3 may block a launch at a gate; downstream nodes are then not reached. This plot is therefore the correct default-visible deviation for the refractory module. Pure timing error of accepted events can remain zero while the response count already differs. """ fig, ax = plt.subplots(figsize=(10, 3.6)) if status.size == 0: ax.set_title("Script 3: fehlende Responses gegen ideale Referenz — keine Daten") return fig accepted_counts = np.sum(status == 1.0, axis=0) launch_count = status.shape[0] deficit = launch_count - accepted_counts blocked_counts = np.sum(status == 0.0, axis=0) unreached_counts = np.sum(status == -1.0, axis=0) ax.plot(positions, deficit, marker="o", linewidth=1.5, markersize=3, label="fehlende akzeptierte Responses") ax.plot(positions, blocked_counts, linestyle="--", linewidth=1.2, label="direkt blockiert") ax.plot(positions, unreached_counts, linestyle=":", linewidth=1.2, label="nicht erreicht") ax.set_xlabel("Knoten / Sitzposition xᵢ") ax.set_ylabel("Anzahl Launches") ax.set_title("Script 3: Response-Defizit relativ zum idealen Script-1-Wellenzug") ax.set_ylim(-0.05, max(1.0, float(launch_count)) + 0.2) ax.grid(True, alpha=0.3) ax.legend(loc="best") fig.tight_layout() return fig def make_refractory_map(positions: np.ndarray, status: np.ndarray): """Plot launch × node acceptance map. Values: -1 = signal did not reach this node because an earlier gate stopped the wave. 0 = signal reached this node but was blocked by local refractory state. 1 = node accepted and responded. """ fig, ax = plt.subplots(figsize=(10, 3.8)) if status.size == 0: ax.set_title("Refraktär-/Failure-Map: keine Daten") return fig extent = [positions[0], positions[-1], -0.5, status.shape[0] - 0.5] im = ax.imshow(status, aspect="auto", origin="lower", extent=extent, vmin=-1.0, vmax=1.0) cbar = fig.colorbar(im, ax=ax) cbar.set_label("-1 nicht erreicht · 0 blockiert · 1 akzeptiert") ax.set_xlabel("Knoten / Sitzposition xᵢ") ax.set_ylabel("Startimpuls m") ax.set_title("Script 3: Refraktär-/Failure-Map pro Launch und Knoten") fig.tight_layout() return fig def simulate( params: RefractoryParams, ) -> Tuple[np.ndarray, np.ndarray, 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) launches = launch_offsets(params.launch_count, params.launch_interval) event_times, strengths, status, edge_delays, raw_jitter = propagate_events(params) heights = activity_from_events(times, event_times, strengths, params.activation_time, params.damping_per_node) reference_activation = exact_regular_schedule(n, params.tau_mean) reference = activity_from_schedule(times, reference_activation, launches, params.activation_time) return times, positions, launches, event_times, strengths, status, edge_delays, raw_jitter, heights, reference 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, refractory_mode: str, failure_policy: str, wrap_ring: bool, dt: float, t_max: float, snapshot_time: float, seed: int, ): params = RefractoryParams( 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), refractory_mode=str(refractory_mode), # type: ignore[arg-type] failure_policy=str(failure_policy), # type: ignore[arg-type] wrap_ring=bool(wrap_ring), dt=float(dt), t_max=float(t_max), snapshot_time=float(snapshot_time), seed=int(seed), ) ( times, positions, launches, event_times, strengths, status, edge_delays, raw_jitter, heights, reference, ) = simulate(params) first_actual = first_activation_times(event_times) edge_reference = first_wave_reference_from_edges(params.n_people, edge_delays, params.tau_mean) errors = timing_error(first_actual, edge_reference) errors[~np.isfinite(first_actual)] = np.nan tau_emp = float(np.mean(edge_delays)) if edge_delays.size else float(params.tau_mean) sigma_emp = 0.0 if float(params.tau_jitter) == 0.0 else float(np.std(edge_delays)) c_emp = params.seat_spacing / max(tau_emp, 1e-12) sigma2 = float(np.var(edge_delays)) if edge_delays.size else 0.0 d_eff = 0.0 if sigma2 <= 1e-30 else (params.seat_spacing**2 * sigma2) / (2.0 * max(tau_emp, 1e-12) ** 3) r_transport, r_diffusion = transport_diffusion_residuals(times, positions, heights, c_emp, d_eff) no_extra_irregularity = (float(params.tau_jitter) == 0.0) and (float(params.damping_per_node) == 0.0) full_recovery_gap = max(float(params.activation_time), 0.0) + max(float(params.refractory_time), 0.0) gate_should_matter = float(params.launch_interval) < full_recovery_gap heat_title = ( "Script 3: Refraktär-Gate aktiv — Wellen können lokal scheitern" if gate_should_matter else "Script 3: Refraktär-Gate aktiv — Default lässt vollständige Erholung zu" ) heatmap = make_heatmap(times, positions, heights, heat_title) snapshot = make_snapshot_comparison( times, positions, heights, reference, params.snapshot_time, "Script 3 simuliert", ) schedule = make_timing_error_plot( positions, errors, "Script 3: erster akzeptierter Launch — Timing-Fehler relativ zur edge-basierten Referenz", ) response_deficit = make_response_deficit_plot(positions, status) refractory_map = make_refractory_map(positions, status) 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)) accepted = int(np.sum(status == 1.0)) blocked = int(np.sum(status == 0.0)) unreached = int(np.sum(status == -1.0)) total_attempt_slots = int(status.size) pass_rate = accepted / max(total_attempt_slots, 1) full_launches = int(np.sum(np.all(status == 1.0, axis=1))) if status.size else 0 failed_launches = int(params.launch_count) - full_launches first_block_per_launch = [] for m in range(status.shape[0]): blocked_nodes = np.where(status[m, :] == 0.0)[0] first_block_per_launch.append(int(blocked_nodes[0]) if blocked_nodes.size else None) first_blocks_text = ", ".join("–" if x is None else str(x) for x in first_block_per_launch) max_abs_error = float(np.nanmax(np.abs(errors))) if np.any(np.isfinite(errors)) else float("nan") max_abs_delta_snapshot = float(np.max(np.abs(heights[snapshot_idx, :] - reference[snapshot_idx, :]))) if heights.size else 0.0 launch_spacing_nodes = params.launch_interval / max(tau_emp, 1e-12) refractory_warning = "Ja — spätere Launches können in noch nicht erholte Knoten laufen." if gate_should_matter else "Nein — Default-Abstand lässt volle Erholung zu." zero_note = "Ja — bis auf das neue Refraktär-Gate." if no_extra_irregularity else "Nein — Timing-Streuung oder Dämpfung ist aktiv." interpretation = f""" ## Script 3 — Interpretation **Neue lokale Regel: Refraktär-Gate / endliche Responsefähigkeit** Script 1 und 2 zeichnen Aktivierungen im Wesentlichen aus einem Schedule. Script 3 macht daraus zum ersten Mal eine lokale Zustandsdynamik: ```text ansprechbar → aktiv → refraktär → ansprechbar Signal kommt bei Knoten i an. Knoten i antwortet nur, wenn seit seiner letzten Aktivierung mindestens T_act + T_ref vergangen ist. ``` Minimal-binär: ```text response_i = 1, falls Δt_i ≥ T_act + T_ref response_i = 0, sonst ``` Bei `partial_recovery` wird die Erholung linear angenähert: ```text response_i = clamp((Δt_i - T_act) / T_ref, 0, 1) ``` | Größe | Wert | |---|---:| | Knotenanzahl N | {params.n_people} | | Startimpulse launch_count | {params.launch_count} | | Startabstand T_launch | {params.launch_interval:.3f} s | | Erholungsschwelle T_act + T_ref | {full_recovery_gap:.3f} s | | Startabstand in Knoten T_launch/mean(τ_i) | {launch_spacing_nodes:.2f} Knoten | | Refraktärwarnung T_launch < T_act+T_ref | {refractory_warning} | | Refraktärmodus | {params.refractory_mode} | | Failure-Policy | {params.failure_policy} | | Ringmodus | {params.wrap_ring} | | sichtbare Pakete im Snapshot | {visible_packets} | | aktive Knoten im Snapshot | {active_nodes} | | akzeptierte Knoten-Launch-Ereignisse | {accepted} | | blockierte Ereignisse | {blocked} | | nicht erreichte Ereignisse | {unreached} | | Pass-Rate akzeptiert / Slots | {pass_rate:.3f} | | vollständig durchgelaufene Launches | {full_launches} | | gescheiterte Launches | {failed_launches} | | erster Block pro Launch | {first_blocks_text} | | mean(τ_i) | {tau_emp:.4f} s | | std(τ_i) | {sigma_emp:.4e} s | | Dämpfung pro Knoten α | {params.damping_per_node:.5f} | | max. Timing-Fehler des ersten akzeptierten Launchs | {max_abs_error:.4e} s | | max. Snapshot-Abweichung gegen Script 1 | {max_abs_delta_snapshot:.4e} | | heuristisches D_eff aus Timing-Streuung | {d_eff:.6e} | | normierter Transport-Residualindikator | {r_transport:.4f} | | normierter Transport-Diffusions-Residualindikator | {r_diffusion:.4f} | | keine Script-2-Irregularität aktiv? | {zero_note} | **Lesart:** Eine Welle ist jetzt nicht mehr garantiert. Die Propagation hängt vom lokalen Live-Zustand ab. Damit wird aus einem vorgegebenen Wellenzug eine Response-Kette mit Gedächtnis. Das ist der erste relevante Schritt in Richtung CNNA-Sprache: lokale Antwortfähigkeit, Relaxation, Record/live-Unterscheidung und Scheitern einer Response-Propagation. **Wichtiger Default-Hinweis:** Script 3 verwendet jetzt als Demonstrations-Default `T_launch = 2.0 s`. Da `T_act + T_ref = 3.2 s`, laufen einige spätere Launches in noch nicht erholte Knoten. Die gemeinsame Script-1/2-Baseline erhältst du durch `T_launch = 4.0 s`; dann bleibt das Gate korrekt, aber im Plot weitgehend unsichtbar. **Warum der Timing-Plot trotzdem null sein kann:** Bei `std(ξ)=0` verschiebt das Refraktär-Gate akzeptierte Ereignisse nicht zeitlich; es akzeptiert oder blockiert sie. Die relevante Abweichung ist daher nicht Timing-Drift, sondern Response-Defizit und Failure-Map. **Unsicherheit / Modellgrenze:** Das ist noch kein Schur-/DtN-Netzwerk. Die Refraktärregel ist eine kontrollierte lokale Automat-Regel. Im späteren CNNA-Modul müsste die Antwortfähigkeit aus der live-Randantwort beziehungsweise einer lokalen Schur-/DtN-Response abgeleitet werden. """ return heatmap, snapshot, schedule, response_deficit, refractory_map, interpretation def build_tab() -> None: gr.Markdown( r""" ## Script 3 — Refraktär-Gate: lokale Response kann scheitern Dieses Modul baut auf Script 1/2 auf, fügt aber ein neues wesentliches Detail hinzu: Ein Knoten reagiert nicht automatisch. Er reagiert nur, wenn seine lokale live-Response wieder frei ist. Der Script-3-Default ist absichtlich ein Refraktär-Stresstest (`T_launch = 2.0 s`), damit das neue Gate sofort sichtbar wird. Die gemeinsame Script-1/2-Baseline ist weiterhin verfügbar: `T_launch = 4.0 s` setzen. """ ) 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.0, 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=SCRIPT3_DEMO_LAUNCH_INTERVAL, step=0.10, label="T_launch Abstand der Startimpulse [s]") damping_per_node = gr.Slider(0.00, 0.04, value=0.0, step=0.001, label="α Dämpfung pro Knoten") refractory_mode = gr.Radio( choices=["binary", "partial_recovery"], value="binary", label="Refraktärmodus", ) failure_policy = gr.Radio( choices=["block", "attenuate"], value="block", label="Failure-Policy", ) wrap_ring = gr.Checkbox(value=False, label="Ringmodus statt offener Reihe") 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=SCRIPT3_DEMO_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="Timing-Fehler akzeptierter erster Durchlauf") response_deficit = gr.Plot(label="Response-Defizit gegen ideale Referenz") refractory_map = gr.Plot(label="Refraktär-/Failure-Map") interpretation = gr.Markdown() inputs = [ n_people, tau_mean, tau_jitter, activation_time, refractory_time, seat_spacing, launch_count, launch_interval, damping_per_node, refractory_mode, failure_policy, wrap_ring, dt, t_max, snapshot_time, seed, ] outputs = [heatmap, snapshot, schedule, response_deficit, refractory_map, interpretation] run.click(run_model, inputs=inputs, outputs=outputs) if __name__ == "__main__": params = RefractoryParams() t, x, launches, events, strengths, status, delays, jitter, h, ref = simulate(params) print( t.shape, x.shape, launches.shape, events.shape, strengths.shape, status.shape, delays.shape, jitter.shape, h.shape, ref.shape, )