#!/usr/bin/env python3

import math import argparse import json import copy import time import numpy as np from statistics import mean, stdev from typing import Dict, Any, List, Tuple

def parse_args(): parser = argparse.ArgumentParser(description="Fused ZSF adaptive injection + microcosm composite operator model") # Core toggles parser.add_argument("--debug", action="store_true", help="Enable debug prints") parser.add_argument("--high-noise", action="store_true", help="Increase noise for easier transitions") parser.add_argument("--lower-threshold", action="store_true", help="Lower TRIT thresholds") parser.add_argument("--stabilize", action="store_true", help="Enable stabilizing composite operator") parser.add_argument("--exact-cancel", action="store_true", help="Exact cancel on collision for opposite signs") parser.add_argument("--search-candidates", type=int, default=0, help="If >0, search N candidate syntheses and select best") parser.add_argument("--regression-test", action="store_true", help="Run short regression test comparing baseline single-point behavior") parser.add_argument("--seed", type=int, default=20251204, help="Global seed") # Injection/coupling & controller parser.add_argument("--inj-mult", type=float, default=3.0, help="Injection amplitude multiplier (3-5)") parser.add_argument("--coupling-mult", type=float, default=3.0, help="Spatial coupling multiplier (2-3)") parser.add_argument("--Kp", type=float, default=0.05, help="Adaptive injection proportional gain") parser.add_argument("--Imax", type=float, default=0.5, help="Adaptive injection max absolute value") parser.add_argument("--theta-target-offset", type=float, default=-0.02, help="Target threshold offset relative to thetaplus") parser.add_argument("--inj-hold", type=float, default=0.05, help="Minimum adaptive injection hold after lock") parser.add_argument("--steps", type=int, default=200, help="Single-point evolution steps") # Grid/microcosm params parser.add_argument("--grid-Nx", type=int, default=128, help="Grid points (1D)") parser.add_argument("--Lx", type=float, default=1.0, help="Domain half-length") parser.add_argument("--T-total", type=float, default=0.01, help="Total simulation time") parser.add_argument("--dt-micro", type=float, default=1e-5, help="Micro time step") parser.add_argument("--dt-macro", type=float, default=1e-3, help="Macro operator sync interval") parser.add_argument("--v0", type=float, default=0.1, help="Base advection velocity") parser.add_argument("--alpha-s", type=float, default=0.5, help="Velocity dependence on state") parser.add_argument("--nu", type=float, default=1e-4, help="Diffusion coefficient") parser.add_argument("--gamma-env", type=float, default=0.05, help="Environmental damping") parser.add_argument("--kappa-proj", type=float, default=0.05, help="Projection coupling for residual collection") parser.add_argument("--Nmcsample", type=int, default=50, help="Monte Carlo samples for residual collision") parser.add_argument("--sigma-noise", type=float, default=1e-6, help="Noise (grid)") parser.add_argument("--Z-threshold", type=float, default=1e-2, help="Energy drift threshold") parser.add_argument("--safetydamp-coef", type=float, default=0.9, help="Safety damping multiplier") parser.add_argument("--cfl-max", type=float, default=0.4, help="Max CFL number for adaptive dt") parser.add_argument("--maxGtotscalar", type=float, default=0.1, help="Clamp total source scalar") parser.add_argument("--maxinjectionfrac", type=float, default=0.05, help="Clamp per-step fractions") parser.add_argument("--min-denominator", type=float, default=1e-12, help="Min denominator to avoid div-by-zero") parser.add_argument("--maxabsorbenergy-frac", type=float, default=0.1, help="Max energy absorbed to correct drift") parser.add_argument("--result-save-path", type=str, default="./fused_result.json", help="Where to save final result JSON") parser.add_argument("--nruns", type=int, default=10, help="Comparison runs") parser.add_argument("--use-grid", action="store_true", help="Enable grid microcosm evolution (original)") return parser.parse_args()

args = parse_args() RNG = np.random.default_rng(args.seed)

def clamp_val(x: float, a: float, b: float) -> float: lo, hi = (a, b) if a <= b else (b, a) return max(lo, min(hi, x))

def rng_default(seed=None) -> np.random.Generator: return np.random.default_rng(seed if seed is not None else 42)

def energy_of_state(S: np.ndarray, dx: float = 1.0, dy: float = None, maxclip: float = 1e6) -> float: Sc = np.clip(S, -maxclip, maxclip) if dy is None: return 0.5 * float(np.sum(Sc * Sc)) * dx return 0.5 * float(np.sum(Sc * Sc)) * dx * dy

def xi_of(t: float, xi0: float = 1.0, amp: float = 0.5, freq: float = 0.01, mode: str = "sin") -> float: if mode == "sin": return xi0 * (1.0 + amp * math.sin(2.0 * math.pi * freq * t)) if mode == "cos": return xi0 * (1.0 + amp * math.cos(2.0 * math.pi * freq * t)) if mode == "constant": return xi0 raise ValueError(f"Unsupported xi mode: {mode}")

def laplacian_1d(field: np.ndarray, dx: float) -> np.ndarray: return (np.roll(field, -1) - 2.0 * field + np.roll(field, 1)) / (dx * dx)

def adjust_dt_by_cfl_1d(V: np.ndarray, dx: float, currentdt: float, cflmax: float = 0.4, dtmin: float = 1e-10, dtmax: float = 1e-2) -> float: try: vmax = float(np.max(np.abs(V))) except Exception: vmax = 1e-16 vmax = max(vmax, 1e-16) cfl = vmax * currentdt / dx if cfl > cflmax: newdt = max(dtmin, currentdt * 0.5) while vmax * newdt / dx > cflmax and newdt > dtmin: newdt = max(dtmin, newdt * 0.5) return newdt else: newdt = min(dtmax, currentdt * 1.05) if vmax * newdt / dx <= cflmax: return newdt return currentdt

def sample_collision_residual(residual_scalar: float, Nsamples: int, rng: np.random.Generator = None) -> np.ndarray: if rng is None: rng = rng_default() if residual_scalar <= 0 or Nsamples <= 0: return np.zeros(0, dtype=float) fracs = rng.beta(2.0, 5.0, size=int(Nsamples)) samples = fracs * max(0.0, residual_scalar) samples += rng.normal(scale=1e-12, size=samples.shape) return np.maximum(0.0, samples)

SOURCE_A = { "name": "A_core", "params": { "Gamma": 1e-3, "delta": 1e-3, "etaamp": 1e-5, "thetaplus": 0.5, "thetaminus": -0.5, "kappa": 1e-3, "xi": 0.021, "lambda": 0.12, "coefxi": -0.5, "coeflambda": -0.3 }, "operators": { "linear": {"coefxi": -0.5, "coeflambda": -0.3}, "dissipative": {"Gamma": 1e-3}, "nonlinear": {"delta": 1e-3}, "noise": {"etaamp": 1e-5} } }

SOURCE_B = { "name": "B_evolved", "params": { "projectionclip": 1.023, "mixalpha": 0.42, "mixnoisesigma": 9.6e-7, "injectionamp": 9.8e-5, "diffusionnu": 9.87e-4, "advection_v": 0.097, "xi": 0.01055, "lambda": 0.121, "mix_gain": 0.63 }, "operators": { "projection": {"clip": 1.023}, "mix": {"alpha": 0.42}, "injection": {"amp": 9.8e-5}, "diffusion": {"nu": 9.87e-4}, "advection": {"v": 0.097} } }

SOURCE_C = { "name": "C_final2", "params": { "etaamp": 5e-5, "thetaplus": 0.42, "thetaminus": -0.42, "injectionamp": 0.06, "mixgain": 10.0, "hold_strength": 0.12, "xi": 0.021, "lambda": 0.121 }, "operators": { "stabilizer": {"injectionamp": 0.06, "mixgain": 10.0, "hold_strength": 0.12} } }

def numericcollide(a, b, rng: np.random.Generator, exactcancel: bool = False) -> float: if a is None and b is None: return 0.0 if a is None: return float(b) if b is None: return float(a) a = float(a) b = float(b) if a * b < 0: if exactcancel: return 0.0 mag = (abs(a) + abs(b)) * 0.5 return float(rng.normal(0.0, max(mag * 0.02, 1e-12))) avg = 0.5 * (a + b) perturb = rng.normal(0.0, max(abs(avg) * 0.03, 1e-12)) return float(avg + perturb)

def collidedict(pa: Dict[str, Any], pb: Dict[str, Any], rng: np.random.Generator, exactcancel: bool = False) -> Dict[str, Any]: keys = set(pa.keys()) | set(pb.keys()) out: Dict[str, Any] = {} for k in keys: va = pa.get(k) vb = pb.get(k) if isinstance(va, (int, float)) or isinstance(vb, (int, float)): out[k] = numericcollide(va, vb, rng, exactcancel) else: out[k] = va if va == vb else (va if va is not None else vb) return out

def collideops(opsA: Dict[str, Dict[str, Any]], opsB: Dict[str, Dict[str, Any]], rng: np.random.Generator, exactcancel: bool = False) -> Dict[str, Dict[str, Any]]: new_ops: Dict[str, Dict[str, Any]] = {} namesA = list(opsA.keys()) namesB = list(opsB.keys()) # Same-name operator collision for name in set(namesA) & set(namesB): new_ops[name] = collidedict(opsA[name], opsB[name], rng, exactcancel) # Reverse-operator collision (pairwise reverse weighting) for i, nameA in enumerate(namesA): partner = namesB[-1 - (i % len(namesB))] newname = f"{nameA}opp{partner}" collided = collidedict(opsA[nameA], opsB[partner], rng, exactcancel) for pk, pv in list(collided.items()): if isinstance(pv, (int, float)): collided[pk] = -pv new_ops[newname] = collided return new_ops

def generatecandidates(sources: List[Dict[str, Any]], rng: np.random.Generator, n: int, exactcancel: bool = False, injmult: float = 1.0, couplingmult: float = 1.0, grid_scaling: bool = False) -> List[Dict[str, Any]]: candidates = [] for _ in range(n): A = sources[rng.integers(0, len(sources))] B = sources[rng.integers(0, len(sources))] pa = {k: (v * (1.0 + rng.normal(0.0, 0.01)) if isinstance(v, (int, float)) else v) for k, v in A["params"].items()} pb = {k: (v * (1.0 + rng.normal(0.0, 0.01)) if isinstance(v, (int, float)) else v) for k, v in B["params"].items()} newparams = collidedict(pa, pb, rng, exactcancel) newops = collideops(A["operators"], B["operators"], rng, exactcancel) # Multipliers if grid_scaling: newparams["injectionamp"] = float(newparams.get("injectionamp", newparams.get("injection_amp", 1e-4))) * float(injmult) if "diffusionnu" in newparams: newparams["diffusionnu"] = float(newparams.get("diffusionnu", 1e-3)) * float(couplingmult) else: newparams["kappa"] = float(newparams.get("kappa", 1e-3)) * float(couplingmult) else: newparams["injectionamp"] = float(newparams.get("injectionamp", newparams.get("injection_amp", 1e-4))) if "diffusionnu" not in newparams: newparams["kappa"] = float(newparams.get("kappa", 1e-3)) newparams["xi"] = float(clamp_val(newparams.get("xi", 0.0), 1e-4, 1.0)) newparams["lambda"] = float(clamp_val(newparams.get("lambda", 0.0), 1e-4, 1.0)) candidates.append({"params": newparams, "operators": newops}) return candidates

def EVOLVEZSFtest(zs: Dict[str, Any], params: Dict[str, Any], rng: np.random.Generator, dt: float = 1e-4) -> Dict[str, Any]: xi = float(params.get("xi", zs.get("xi", 0.021))) lam = float(params.get("lambda", zs.get("lambda", 0.121))) phi = float(zs.get("phi", 0.0)) linearxi = float(params.get("coefxi", params.get("linearcoefxi", -0.5))) linearlam = float(params.get("coeflambda", params.get("linearcoeflambda", -0.3))) linearterm = linearxi * xi + linearlam * lam Gamma = float(params.get("Gamma", params.get("dissipative_Gamma", 5e-4))) dissipative = -Gamma * phi delta = float(params.get("delta", 1e-3)) nonlinear = -delta * (phi * 3.0) eta = float(params.get("etaamp", 1e-5 if not args.high_noise else 5e-5)) noise = eta * rng.normal(0.0, 1.0) / math.sqrt(max(dt, 1e-16)) thetaplus = float(params.get("thetaplus", 0.35 if args.lower_threshold else 0.5)) thetaminus = float(params.get("thetaminus", -0.35 if args.lower_threshold else -0.5)) inj = float(params.get("injectionamp", 0.0) or 0.0) injection = 0.0 near = float(params.get("near_margin", 0.12)) if (thetaplus - near) <= phi < thetaplus: injection += inj * (1.0 + 5.0 * xi) if thetaminus < phi <= (thetaminus + near): injection -= inj * (1.0 + 5.0 * xi) mixgain = float(params.get("mixgain", params.get("mix_alpha", 5.0))) if abs(phi) < float(params.get("mix_zone", 0.25)): injection += mixgain * eta * rng.normal(0.0, 1.0) phinew = phi + dt * (linearterm + dissipative + nonlinear) + math.sqrt(max(dt, 1e-16)) * noise + injection # TRIT map + hold if phinew >= thetaplus: trit = 1 elif phinew <= thetaminus: trit = -1 else: trit = 0 if trit == 1: hold = float(params.get("hold_strength", 0.12)) phinew = phinew * (1.0 - hold) + 1.0 * hold elif trit == -1: hold = float(params.get("hold_strength", 0.12)) phinew = phinew * (1.0 - hold) - 1.0 * hold zs_out = dict(zs) zs_out["phi"] = float(phinew) zs_out["tritstate"] = int(trit) zs_out["xi"] = xi zs_out["lambda"] = lam zs_out["entropy"] = zs_out["tritstate"] * xi return zs_out

def validate_candidates(cands: List[Dict[str, Any]], rng: np.random.Generator, steps: int = 500, runs: int = 4, debug: bool = False) -> List[Dict[str, Any]]: results = [] for idx, c in enumerate(cands): success = 0 finals: List[int] = [] for r in range(runs): rgen = np.random.default_rng(int(rng.integers(1, 10000))) zs = {"phi": 0.0, "tritstate": 0, "xi": c["params"].get("xi", 0.021), "lambda": c["params"].get("lambda", 0.121)} for _ in range(steps): zs = EVOLVEZSFtest(zs, c["params"], rgen) finals.append(zs["tritstate"]) if zs["tritstate"] != 0: success += 1 rate = success / runs results.append({"idx": idx, "rate": rate, "finals": finals, "params": c["params"]}) if debug and idx < 6: print(f"[validate] idx={idx} rate={rate:.2f} keys={list(c['params'].keys())[:6]}") return results

def applystabilizingops(zs: Dict[str, Any], rng: np.random.Generator, params: Dict[str, Any] = None) -> Dict[str, Any]: if params is None: params = {} phi = float(zs.get("phi", 0.0)) xi = float(zs.get("xi", 0.021)) thetaplus = float(params.get("thetaplus", 0.35 if args.lower_threshold else 0.5)) thetaminus = float(params.get("thetaminus", -0.35 if args.lower_threshold else -0.5)) softclip = float(params.get("softclip", 1.0)) clipstrength = float(params.get("clipstrength", 0.15)) if abs(phi) > softclip: phi = phi * (1.0 - clipstrength) + math.copysign(softclip, phi) * clipstrength biasamp = float(params.get("injectionamp", 0.06)) nearmargin = float(params.get("nearmargin", 0.12)) if (thetaplus - nearmargin) <= phi < thetaplus: phi += biasamp * (1.0 + 5.0 * xi) if thetaminus < phi <= (thetaminus + nearmargin): phi -= biasamp * (1.0 + 5.0 * xi) mixzone = float(params.get("mixzone", 0.25)) mixgain = float(params.get("mixgain", 10.0)) if abs(phi) < mixzone: noiseamp = float(params.get("etaamp", 5e-5)) noise = noiseamp * rng.normal(0.0, 1.0) drive = float(params.get("drive_bias", 0.0)) phi += (noise * mixgain) + drive * 1e-3 trit = int(zs.get("tritstate", 0)) holdstrength = float(params.get("holdstrength", 0.12)) if trit == 1: phi = phi * (1.0 - holdstrength) + 1.0 * holdstrength elif trit == -1: phi = phi * (1.0 - holdstrength) - 1.0 * holdstrength phi = max(min(phi, 1.2), -1.2) zs_out = dict(zs) zs_out["phi"] = float(phi) return zs_out

def EVOLVEZSFwithadaptiveinjection(zs: Dict[str, Any], params: Dict[str, Any], dt: float = 1e-4, rng: np.random.Generator = None) -> Dict[str, Any]: if rng is None: rng = rng_default(zs.get("seed", 42)) xi = float(zs.get("xi", params.get("xi", 0.021))) lam = float(zs.get("lambda", params.get("lambda", 0.121))) phi = float(zs.get("phi", 0.0)) coefxi = float(params.get("coefxi", -0.5)) coeflambda = float(params.get("coeflambda", -0.3)) linearterm = coefxi * xi + coeflambda * lam Gamma = float(params.get("Gamma", 5e-4)) dissipative_term = -Gamma * phi delta = float(params.get("delta", 1e-3)) nonlinear_term = -delta * (phi * 3.0) etaamp = float(params.get("etaamp", 1e-5 if not args.high_noise else 5e-5)) noise_term = etaamp * rng.normal(0.0, 1.0) / math.sqrt(max(dt, 1e-16)) # Sampled base injection residual_scalar = max(0.0, phi * float(params.get("Nmcsample", args.Nmcsample))) samples = sample_collision_residual(residual_scalar, int(params.get("Nmcsample", args.Nmcsample)), rng) inj_sample = (float(np.mean(samples)) * float(params.get("eta_in", 0.05))) if samples.size > 0 else float(params.get("injectionamp", 1e-4)) injbase = inj_sample * (1.0 + 5.0 * xi) # Adaptive injection injadapt = float(zs.get("injectionadaptive", 0.0)) Kp = float(params.get("injcontrollergain", args.Kp)) Imax = float(params.get("injmax", args.Imax)) thetaplus = float(params.get("thetaplus", 0.35 if args.lower_threshold else 0.5)) thetaminus = float(params.get("thetaminus", -0.35 if args.lower_threshold else -0.5)) thetatarget = thetaplus + float(args.theta_target_offset) previnj = float(zs.get("injprev", 0.0)) Kd = float(params.get("inj_damping", 0.01)) error = thetatarget - phi injadapt += Kp * error - Kd * (injadapt - previnj) injadapt = max(min(injadapt, Imax), -Imax) injection = injbase + injadapt if abs(phi) < float(params.get("mix_zone", 0.25)): injection += float(params.get("mixgain", 5.0)) * etaamp * rng.normal(0.0, 1.0) phinew = phi + dt * (linearterm + dissipative_term + nonlinear_term) + math.sqrt(max(dt, 1e-16)) * noise_term + injection # TRIT map + soft-hold + adaptive bump if phinew >= thetaplus: trit = 1 hold = float(params.get("hold_strength", 0.12)) phinew = phinew * (1.0 - hold) + 1.0 * hold elif phinew <= thetaminus: trit = -1 hold = float(params.get("hold_strength", 0.12)) phinew = phinew * (1.0 - hold) - 1.0 * hold else: trit = 0 if trit == 1: injadapt = max(injadapt, float(params.get("injhold", args.inj_hold))) injadapt *= 0.98 prevtrit = int(zs.get("tritstate", 0)) if prevtrit == 1 and trit == 0 and phinew < thetaplus - 0.05: injadapt += min(0.1, Imax * 0.2) phinew = max(min(phinew, 1.2), -1.2) zs_out = dict(zs) zs_out["phi"] = float(phinew) zs_out["tritstate"] = int(trit) zs_out["injectionadaptive"] = float(injadapt) zs_out["injprev"] = float(injadapt) zs_out["xi"] = xi zs_out["lambda"] = lam zs_out["entropy"] = zs_out["tritstate"] * xi if args.debug: print(f"ADAPT step: phi={phi:.6e} -> {phinew:.6e}, injbase={injbase:.3e}, injadapt={injadapt:.3e}, trit={trit}") return zs_out

class Operator: def init(self, name: str, params: Dict[str, Any], apply_fn): self.name = name self.params = dict(params) self.apply_fn = apply_fn self.next_op: Operator | None = None

def apply(self, state: np.ndarray, rng: np.random.Generator = None) -> Tuple[np.ndarray, Dict[str, Any]]:
    if rng is None:
        rng = rng_default()
    result, info = self.apply_fn(state, self.params, rng)
    if self.next_op is not None:
        result, next_info = self.next_op.apply(result, rng)
        info.update(next_info)
    return result, info

def chain(self, next_op: "Operator") -> "Operator":
    self.next_op = next_op
    return self

def op_diffusion(state, params, rng): dx = float(params.get("dx", 1.0)) dt = float(params.get("dt", 1e-5)) nu = float(params.get("nu", 1e-4)) * float(params.get("coupling_mult", args.coupling_mult)) lap = laplacian_1d(state, dx) state_new = state + nu * lap * dt return state_new, {"energy_diffusion": energy_of_state(state_new, dx), "nu_used": nu}

def op_advection(state, params, rng): dx = float(params.get("dx", 1.0)) dt = float(params.get("dt", 1e-5)) v = float(params.get("v", 0.097)) state_avg = float(np.mean(state)) v = v * (1.0 + 0.1 * state_avg) # adaptive speed F = v * state F_im = np.roll(F, 1) conv_term = -(dt / dx) * (F - F_im) state_new = state + conv_term return state_new, {"energy_advection": energy_of_state(state_new, dx), "v_used": v}

def op_noise(state, params, rng): dx = float(params.get("dx", 1.0)) sigma = float(params.get("sigma", 1e-6)) * (5 if params.get("high_noise", False) else 1) noise = rng.normal(scale=sigma, size=state.shape) state_new = state + noise return state_new, {"energy_noise": energy_of_state(state_new, dx), "sigma_used": sigma}

def op_projection(state, params, rng): dx = float(params.get("dx", 1.0)) clip = float(params.get("clip", 1.023)) soft_clip = float(params.get("soft_clip", 1.2)) state_new = np.where(np.abs(state) > soft_clip, state * (1.0 - 0.15) + np.sign(state) * soft_clip * 0.15, state) state_new = np.clip(state_new, -clip, clip) return state_new, {"energy_projection": energy_of_state(state_new, dx), "clipped_count": int(np.sum(np.abs(state_new) >= clip))}

def op_adaptive_injection(state, params, rng): dx = float(params.get("dx", 1.0)) grid_size = int(state.shape[0]) center = grid_size // 2 width = max(1.0, float(params.get("width", max(1, grid_size // 20)))) x = np.arange(grid_size) kernel = np.exp(-0.5 * ((x - center) / width) ** 2) kernel = kernel / (kernel.sum() + 1e-12) phi_prev = float(np.mean(state)) thetaplus = 0.35 if params.get("lower_threshold", False) else 0.5 theta_target = thetaplus + float(args.theta_target_offset) Kp = float(params.get("Kp", args.Kp)) Kd = float(params.get("Kd", 0.01)) Imax = float(params.get("Imax", args.Imax)) inj_prev = float(params.get("inj_prev", 0.0)) inj_prev_prev = float(params.get("inj_prev_prev", 0.0)) error = theta_target - phi_prev inj_adapt = inj_prev + Kp * error - Kd * (inj_prev - inj_prev_prev) inj_adapt = clamp_val(inj_adapt, -Imax, Imax) inj_base = float(params.get("injectionamp", 1e-4)) * (1.0 + 5.0 * float(params.get("xi", 0.021))) * float(params.get("inj_mult", args.inj_mult)) mix_zone = float(params.get("mix_zone", 0.25)) if abs(phi_prev) < mix_zone: inj_base += float(params.get("mixgain", 5.0)) * float(params.get("etaamp", 1e-5)) * rng.normal(0.0, 1.0) injection = inj_base + inj_adapt state_new = state + injection * kernel trit = 1 if phi_prev >= thetaplus else (-1 if phi_prev <= -thetaplus else 0) if trit == 1: inj_adapt = max(inj_adapt, float(params.get("inj_hold", args.inj_hold))) inj_adapt *= 0.98 params["inj_prev_prev"] = inj_prev params["inj_prev"] = inj_adapt return state_new, {"energy_injection": energy_of_state(state_new, dx), "inj_adapt": inj_adapt, "trit": trit}

def op_stabilize(state, params, rng): dx = float(params.get("dx", 1.0)) phi = float(np.mean(state)) thetaplus = float(params.get("thetaplus", 0.35 if args.lower_threshold else 0.5)) soft_clip = float(params.get("softclip", 1.0)) clip_strength = float(params.get("clipstrength", 0.15)) if abs(phi) > soft_clip: phi = phi * (1.0 - clip_strength) + math.copysign(soft_clip, phi) * clip_strength E_prev = energy_of_state(state, dx) scale = phi / (float(np.mean(state)) + 1e-12) state_new = state * scale E_new = energy_of_state(state_new, dx) Z_proxy = abs(E_new - E_prev) / max(float(params.get("Z_threshold", args.Z_threshold)), abs(E_prev)) if Z_proxy > float(params.get("Z_threshold", args.Z_threshold)): state_new = state_new * float(params.get("safetydamp_coef", args.safetydamp_coef)) return state_new, {"energy_stabilize": energy_of_state(state_new, dx), "Z_proxy": Z_proxy}

def create_evolved_best_operator(dx: float) -> Operator: default_ops = [ Operator("advection", {"v": 0.097, "dx": dx, "dt": args.dt_micro}, op_advection), Operator("noise", {"sigma": args.sigma_noise, "dx": dx, "high_noise": args.high_noise}, op_noise), Operator("diffusion", {"nu": 9.87e-4, "dx": dx, "dt": args.dt_micro, "coupling_mult": args.coupling_mult}, op_diffusion), Operator("adaptive_injection", { "injectionamp": 9.8e-5, "inj_mult": args.inj_mult, "dx": dx, "lower_threshold": args.lower_threshold, "xi": 0.021, "mixgain": 5.0, "etaamp": 1e-5, "Kp": args.Kp, "Imax": args.Imax }, op_adaptive_injection), Operator("projection", {"clip": 1.023, "soft_clip": 1.2, "dx": dx}, op_projection), Operator("stabilize", {"Z_threshold": args.Z_threshold, "dx": dx, "safetydamp_coef": args.safetydamp_coef}, op_stabilize), ] fused_op = default_ops[0] for op in default_ops[1:]: fused_op = fused_op.chain(op) return fused_op

def EVOLVEZSFwithadaptiveinjection_grid(state: np.ndarray, params: Dict[str, Any], rng: np.random.Generator, dx: float, dt: float) -> Tuple[np.ndarray, Dict[str, Any]]: op_chain = create_evolved_best_operator(dx) state_new, info = op_chain.apply(state, rng) phi = float(np.mean(state_new)) thetaplus = float(params.get("thetaplus", 0.35 if args.lower_threshold else 0.5)) thetaminus = -thetaplus trit = 1 if phi >= thetaplus else (-1 if phi <= thetaminus else 0) info["trit_state"] = trit info["phi"] = phi info["dt_used"] = dt return state_new, info

def FUSE(zs1: Dict[str, Any], zs2: Dict[str, Any], T_local: float = 300.0) -> Dict[str, Any]: b1 = int(zs1.get("tritstate", zs1.get("trit_state", 0))) b2 = int(zs2.get("tritstate", zs2.get("trit_state", 0))) xi1 = float(zs1.get("xi", 0.021)) xi2 = float(zs2.get("xi", 0.021)) fused: Dict[str, Any] = {} xor = (b1 ^ b2) % 2 if xor != 0: delta_xi = (xi1 - xi2) / 2.0 zs1["xi"] = xi1 - delta_xi zs2["xi"] = xi2 + delta_xi sum_b = (b1 + b2) % 3 if sum_b != 0: delta_S = ((b1 + b2) - 0) * zs1.get("xi", xi1) deltaE = abs(delta_S) * T_local fused["energy"] = zs1.get("energy", 0.0) + zs2.get("energy", 0.0) + deltaE fused_trit = int(round((b1 + b2) / 2.0)) if (b1 + b2) != 0 else 0 fused["tritstate"] = fused_trit fused["xi"] = (zs1.get("xi", xi1) + zs2.get("xi", xi2)) / 2.0 fused["entropy"] = fused["tritstate"] * fused["xi"] if args.debug: print(f"DEBUG FUSE: b1={b1}, b2={b2}, fusedtrit={fused_trit}, fused_xi={fused['xi']:.6e}") return fused

def FUSE_grid(zs_list: List[Dict[str, Any]], T_local: float, dx: float) -> Dict[str, Any]: if len(zs_list) < 2: return zs_list[0] if zs_list else {"state": np.zeros(10), "energy": 0.0, "trit_state": 0} energies = [energy_of_state(zs["state"], dx) for zs in zs_list] weights = [1.0 / max(E, 1e-12) for E in energies] weights = [w / sum(weights) for w in weights] fused_state = np.sum([zs["state"] * w for zs, w in zip(zs_list, weights)], axis=0) trits = [zs.get("trit_state", 0) for zs in zs_list] fused_trit = int(round(np.sum([t * w for t, w in zip(trits, weights)]))) fused_energy = energy_of_state(fused_state, dx) delta_S = float(np.sum([zs.get("entropy", 0.0) * w for zs, w in zip(zs_list, weights)])) fused_energy += abs(delta_S) * T_local return {"state": fused_state, "energy": fused_energy, "trit_state": fused_trit, "entropy": delta_S, "weights": weights}

def DARKMATTERTRIGGER(zs: Dict[str, Any], couplingsection: float, Veff: float, dt: float = 1e-4, rng: np.random.Generator = None, maxsteps: int = 1000, params: Dict[str, Any] = None) -> Tuple[Dict[str, Any], List[str]]: if rng is None: rng = rng_default(zs.get("seed", 42)) if params is None: params = {} log: List[str] = [] if couplingsection <= 1e-40: log.append("coupling_section below threshold; no trigger") return zs, log zs = dict(zs) zs["lambda"] = 0.121 log.append("lambda set to 0.121") zs["phi"] = zs.get("phi", 0.0) + 0.1 log.append(f"phi biased +0.1 -> {zs['phi']:.6e}") evo_params = params.get("evoparams", {}) for step in range(maxsteps): zs = EVOLVEZSFwithadaptiveinjection(zs, evo_params, dt=dt, rng=rng) if zs.get("tritstate", 0) == 1: log.append(f"Reached trit=1 at step {step}") break delta_S = 4e-4 T_local = 300.0 deltaE = delta_S * T_local deltarho = deltaE / (Veff * (3e8 ** 2)) deltaR = 8.0 * math.pi * 6.67e-11 * deltarho / (3e8 ** 2) zs["entropy"] = delta_S zs["energy"] = zs.get("energy", 0.0) + deltaE zs["curvature"] = zs.get("curvature", 0.0) + deltaR log.append(f"trigger success: injected energy={deltaE:.6e} J, curvature change={deltaR:.6e} m^-1") return zs, log

def DARKMATTER_TRIGGER_grid(zs: Dict[str, Any], coupling_section: float, Veff: float, rng: np.random.Generator, dx: float) -> Tuple[Dict[str, Any], List[str]]: log: List[str] = [] if coupling_section <= 1e-40: log.append("coupling_section below threshold; no trigger") return zs, log state = zs.get("state", np.zeros(10)) zs["lambda"] = 0.121 log.append("lambda set to 0.121") phi_biased = float(np.mean(state)) + 0.1 scale = phi_biased / (float(np.mean(state)) + 1e-12) state = state * scale log.append(f"phi biased +0.1 -> {phi_biased:.6e}") evo_params = zs.get("evoparams", {}) max_steps = 1000 trit_state = 0 for step in range(max_steps): state, info = EVOLVEZSFwithadaptiveinjection_grid(state, evo_params, rng, dx, args.dt_micro) trit_state = info["trit_state"] if trit_state == 1: log.append(f"Reached trit=1 at step {step}") break delta_S = 4e-4 T_local = 300.0 delta_E = delta_S * T_local delta_rho = delta_E / (Veff * (3e8 ** 2)) delta_R = 8.0 * math.pi * 6.67e-11 * delta_rho / (3e8 ** 2) E_current = energy_of_state(state, dx) E_new = E_current + delta_E state = state * math.sqrt(E_new / max(E_current, 1e-12)) Z_proxy = abs(E_new - E_current) / max(args.Z_threshold, abs(E_current)) if Z_proxy > args.Z_threshold: state = state * args.safetydamp_coef log.append(f"Energy injection exceed Z_threshold: {Z_proxy:.3e}, damp state") zs["state"] = state zs["entropy"] = delta_S zs["energy"] = E_new zs["curvature"] = zs.get("curvature", 0.0) + delta_R zs["trit_state"] = trit_state log.append(f"trigger success: injected energy={delta_E:.6e} J, curvature change={delta_R:.6e} m^-1") return zs, log

def quantumsecondquantizescaled(k: float, omegak: float, debug: bool = False, clampforsynthesis: bool = False, scale: float = None, autoscale: bool = False, targetxi: float = 0.021) -> Dict[str, Any]: hbar = 1.05457e-34 creationop = math.sqrt(hbar / (2.0 * omegak)) if omegak > 0 else 0.0 xiraw = creationop * 0.021 xi = xiraw if autoscale and xiraw > 0: if xiraw != 0: scaleauto = targetxi / xiraw xi = xiraw * scaleauto else: xi = xiraw * (targetxi / xiraw if xiraw != 0 else 1.0) if scale is not None: xi = xi * scale if clampforsynthesis: xi = clamp_val(xi, 0.019, 0.023) xi = max(xi, 1e-3) tritstate = 1 if creationop > 0 else -1 return {"xi": xi, "lambda": 0.121, "entropy": 0.0, "tritstate": tritstate}

def ZSFCARBONSYNTHESISdemo(scale: float = None, clamp: bool = False, autoscale: bool = True, debug: bool = False) -> Dict[str, Any]: candidates = [1.2e15, 1.0e15, 1.5e15, 8.0e14] carbon_zs = None for ok in candidates: cz = quantumsecondquantizescaled(0.5, ok, debug=debug, clampforsynthesis=False, scale=scale, autoscale=autoscale, targetxi=0.021) if debug: print(f"TRY omega_k={ok:.6e} -> xi={cz['xi']:.6e}") if 0.019 <= cz["xi"] <= 0.023: carbon_zs = cz if debug: print(f"SELECTED omega_k={ok} with xi={cz['xi']:.6e}") break if carbon_zs is None: carbon_zs = quantumsecondquantizescaled(0.5, candidates[0], debug=debug, clampforsynthesis=clamp, scale=scale, autoscale=autoscale, targetxi=0.021) if debug: print(f"FINAL carbon xi after fallback: {carbon_zs['xi']:.6e}") catalystzs = {"xi": 0.021, "lambda": 0.121, "entropy": -0.0002, "tritstate": -1, "energy": 0.0} fused = FUSE(carbon_zs, catalystzs, T_local=3300.0) carbonstrength = 11.2 if fused["tritstate"] == 0 else 12.0 return {"strengthGPa": carbonstrength, "fusedtrit": fused["tritstate"], "fused": fused, "carbonzs": carbon_zs}

def ZSF_CARBON_SYNTHESIS_grid(dx: float, debug: bool = False) -> Dict[str, Any]: res = ZSFCARBONSYNTHESISdemo(autoscale=True, debug=debug) Nx = int(args.grid_Nx) carbon_state = res["carbonzs"]["xi"] * np.ones(Nx, dtype=float) catalyst_state = 0.05 * np.exp(-0.5 * (np.linspace(-args.Lx, args.Lx, Nx) / 0.1) ** 2) fused_grid = FUSE_grid([ {"state": carbon_state, "entropy": 0.0, "trit_state": res["fusedtrit"]}, {"state": catalyst_state, "entropy": -0.0002, "trit_state": -1} ], T_local=300.0, dx=dx) injection_source = fused_grid["state"] * 0.1 return {"strengthGPa": res["strengthGPa"], "fused_trit": fused_grid["trit_state"], "fused": fused_grid, "carbon_zs": res["carbonzs"], "injection_source": injection_source}

def unitpointevolutiondemo(Nsteps: int = None, seed: int = 42, paramsoverride: Dict[str, Any] = None) -> Dict[str, List[float]]: if Nsteps is None: Nsteps = args.steps rng = rng_default(seed) params = { "Gamma": 5e-4, "delta": 1e-3, "etaamp": 1e-5, "xi": 0.021, "lambda": 0.121, "thetaplus": 0.35 if args.lower_threshold else 0.5, "thetaminus": -0.35 if args.lower_threshold else -0.5, "dt": 1e-4, "Nmcsample": args.Nmcsample, "eta_in": 0.05, "injectionamp": 1e-4 * args.inj_mult, "mixgain": 5.0 * args.inj_mult, "mix_zone": 0.25, "hold_strength": 0.12, "stabilize": args.stabilize, "injcontrollergain": args.Kp, "injmax": args.Imax, "injhold": args.inj_hold } params["kappa"] = float(params.get("kappa", 1e-3)) * float(args.coupling_mult) if paramsoverride: params.update(paramsoverride) zs = {"xi": params["xi"], "lambda": params["lambda"], "phi": 0.0, "tritstate": 0, "entropy": 0.0, "seed": seed, "injectionadaptive": 0.0, "injprev": 0.0} history = {"phi": [], "trit": [], "injadapt": []} for _ in range(Nsteps): zs = EVOLVEZSFwithadaptiveinjection(zs, params, dt=params["dt"], rng=rng) if params.get("stabilize", False) or args.stabilize: zs = applystabilizingops(zs, rng, params.get("stabilize_params", {})) history["phi"].append(zs["phi"]) history["trit"].append(zs["tritstate"]) history["injadapt"].append(zs.get("injectionadaptive", 0.0)) return history

def regression_test(): seed = 12345 baseline = unitpointevolutiondemo(Nsteps=200, seed=seed, paramsoverride=None) merged = unitpointevolutiondemo(Nsteps=200, seed=seed, paramsoverride=None) phi_diff = abs(baseline["phi"][-1] - merged["phi"][-1]) trit_match = baseline["trit"][-1] == merged["trit"][-1] print("Regression test: final phi diff =", phi_diff, "trit match =", trit_match) return {"phi_diff": phi_diff, "trit_match": trit_match, "baseline_final_phi": baseline["phi"][-1], "merged_final_phi": merged["phi"][-1]}

BASE_PARAMS = { "seed": args.seed, "T_total": args.T_total, "dt_micro": args.dt_micro, "dt_macro": args.dt_macro, "grid_Nx": args.grid_Nx, "Lx": args.Lx, "v0": args.v0, "alpha_s": args.alpha_s, "nu": args.nu, "gamma_env": 0.05, "kappa_proj": args.kappa_proj, "eta_in": 0.05, "eta_heat": 0.10, "eta_fus": 0.05, "sigma_fus": 0.03, "Nmcsample": args.Nmcsample, "gammak": 5e-4, "sigma_noise": args.sigma_noise, "Z_threshold": args.Z_threshold, "safetydamp_coef": args.safetydamp_coef, "maxGtotscalar": args.maxGtotscalar, "min_denominator": args.min_denominator, "maxinjectionfrac": args.maxinjectionfrac, "adaptivedt_allowed": True, "cfl_max": args.cfl_max, "maxabsorbenergy_frac": args.maxabsorbenergy_frac, "debugstopon_nan": True, "maxclipvalue": 1e6, }

def upwind_step_1d(S, V, dtlocal, dxlocal, Ginlocal=None, gamma_env=0.0, nu=0.0): if Ginlocal is None: Ginlocal = np.zeros_like(S) F = V * S F_im = np.roll(F, 1) convterm = -(dtlocal / dxlocal) * (F - F_im) diffterm = nu * (dtlocal / (dxlocal * dxlocal)) * laplacian_1d(S, dxlocal) srcterm = dtlocal * (Ginlocal - gamma_env * S) newS = S + convterm + diffterm + srcterm newS = 0.995 * newS + 0.005 * np.roll(newS, 1) return newS

def run_simulation(params, use_evolved=True, verbose=False): p = copy.deepcopy(params) rng = rng_default(int(p["seed"])) Nx = int(p["grid_Nx"]) x = np.linspace(-p["Lx"], p["Lx"], Nx) dx = 2.0 * p["Lx"] / max(1, Nx - 1) S = 0.1 + 0.01 * np.exp(-0.5 * (x / 0.2)**2) V = p["v0"] * (1.0 + p["alpha_s"] * S) Nexpect = 0.02 Eheat = 0.0 Emicro0 = float(Nexpect) Emacro0 = energy_of_state(S, dx) E0 = Emicro0 + Emacro0 + Eheat dt = float(p["dt_micro"]) times = np.arange(0.0, float(p["T_total"]), dt) steps = len(times) macrosync = max(1, int(round(p["dt_macro"] / dt))) projkernel = np.exp(-0.5 * (x / 0.1)**2) projkernel /= projkernel.sum() collkernel = np.exp(-0.5 * (x / p["sigma_fus"])**2) collkernel /= collkernel.sum() hist = { "t": np.zeros(steps), "N_expect": np.zeros(steps), "E_micro": np.zeros(steps), "E_macro": np.zeros(steps), "E_heat": np.zeros(steps), "Z_proxy": np.zeros(steps), "S_center": np.zeros(steps), } localdt = dt evolved_op = create_evolved_best_operator(dx=dx) if use_evolved else None start_time = time.time() for k, t in enumerate(times): xi = xi_of(t) dNdrive = 1e-2 * xi * (1.0 - np.clip(Nexpect, 0.0, 1.0)) * localdt dNdiss = -p["gammak"] * Nexpect * localdt noise = rng.normal(scale=math.sqrt(max(p["sigma_noise"], 0.0)) * math.sqrt(localdt)) Nexpect = max(0.0, Nexpect + dNdrive + dNdiss + noise) Ncollected = p["kappa_proj"] * Nexpect Nresidual = max(0.0, Nexpect - Ncollected) samples = sample_collision_residual(Nresidual, p["Nmcsample"], rng) coll_count = samples.size if coll_count > 0: deltain = p["eta_in"] * samples deltaheat = p["eta_heat"] * samples deltafus = p["eta_fus"] * samples deltainavg = float(deltain.sum() / coll_count) deltaheatavg = float(deltaheat.sum() / coll_count) deltafusavg = float(deltafus.sum() / coll_count) deltainavg = min(deltainavg, p["maxinjectionfrac"] * max(1e-12, Nexpect)) deltaheatavg = min(deltaheatavg, p["maxinjectionfrac"] * max(1e-12, Nexpect)) deltafusavg = min(deltafusavg, p["maxinjectionfrac"] * max(1e-12, Nexpect)) deltainavg = float(np.clip(deltainavg, -1e-3, 1e-3)) deltaheatavg = float(np.clip(deltaheatavg, -1e-3, 1e-3)) deltafusavg = float(np.clip(deltafusavg, -1e-3, 1e-3)) else: deltainavg = deltaheatavg = deltafusavg = 0.0 Gtotscalar = Ncollected + deltainavg Gtotscalar = min(Gtotscalar, p["maxGtotscalar"]) Gproj = Gtotscalar * projkernel Gfus = (1.0 - p["eta_fus"]) * deltafusavg * collkernel Gtotal = Gproj + Gfus # Adjust time step by CFL condition if p["adaptivedt_allowed"]: try: localdt = adjust_dt_by_cfl_1d(V, dx, localdt, p["cfl_max"], 1e-10, p["dt_micro"] * 10) except Exception: localdt = max(1e-10, localdt * 0.5) # Apply operator (evolved or original) if (k % macrosync) == 0: if use_evolved and evolved_op is not None: S, _ = evolved_op.apply(S, rng) else: S = upwind_step_1d(S, V, macrosync * localdt, dx, Ginlocal=Gtotal, gamma_env=p["gamma_env"], nu=p["nu"]) # Update velocity V = p["v0"] * (1.0 + p["alpha_s"] * S) # Update energy components maxallowed2 = max(p["maxinjectionfrac"] * max(1e-12, Nexpect), p["min_denominator"]) heatincrement = min(deltaheatavg, maxallowed2) newmicromass = min(deltafusavg, maxallowed2) Eheat += heatincrement Nexpect = max(0.0, Nexpect - heatincrement + newmicromass) Emicro = float(Nexpect) Emacro = energy_of_state(S, dx) Etotal = Emicro + Emacro + Eheat # Check energy drift Zproxy = abs(Etotal - E0) / max(p["min_denominator"], abs(E0)) if Zproxy > p["Z_threshold"]: deltaE = Etotal - E0 maxabsorb = p.get("maxabsorbenergy_frac", 0.1) * max(abs(E0), p["min_denominator"]) absorb = math.copysign(min(abs(deltaE), maxabsorb), deltaE) Eheat -= absorb S = S * p["safetydamp_coef"] Nexpect = Nexpect * p["safetydamp_coef"] # Record history hist["t"][k] = t hist["N_expect"][k] = Nexpect hist["E_micro"][k] = Emicro hist["E_macro"][k] = Emacro hist["E_heat"][k] = Eheat hist["Z_proxy"][k] = Zproxy hist["S_center"][k] = float(S[len(S)//2]) # Debug: stop on NaN if p["debugstopon_nan"]: if not np.isfinite(Etotal) or np.isnan(Zproxy) or np.any(~np.isfinite(S)): return hist, time.time() - start_time, True runtime = time.time() - start_time return hist, runtime, False

def run_comparison(nruns=10): seeds = [int(BASE_PARAMS["seed"]) + i for i in range(nruns)] results = {"evolved": [], "original": []} for use_evolved in [True, False]: label = "evolved" if use_evolved else "original" for s in seeds: p = copy.deepcopy(BASE_PARAMS) p["seed"] = s hist, runtime, failed = run_simulation(p, use_evolved=use_evolved, verbose=False) E = hist["E_macro"] + hist["E_micro"] + hist["E_heat"] E0 = E[0] errE = np.abs(E - E0) / max(1e-12, abs(E0)) mean_err = float(np.mean(errE)) max_err = float(np.max(errE)) drift = float((E[-1] - E0) / max(1e-12, abs(E0))) results[label].append({ "mean_err": mean_err, "max_err": max_err, "drift": drift, "nan": failed, "runtime": runtime }) print(f"{label} seed={s}: mean_err={mean_err:.3e}, max_err={max_err:.3e}, drift={drift:.3e}, nan={failed}, time={runtime:.3f}s") def summarize(list_of_dicts): mean_errs = [d["mean_err"] for d in list_of_dicts] drifts = [d["drift"] for d in list_of_dicts] nans = sum(1 for d in list_of_dicts if d["nan"]) times = [d["runtime"] for d in list_of_dicts] return { "mean_err_mean": mean(mean_errs), "mean_err_std": stdev(mean_errs) if len(mean_errs) > 1 else 0.0, "drift_mean": mean(drifts), "drift_std": stdev(drifts) if len(drifts) > 1 else 0.0, "nan_count": nans, "time_mean": mean(times) } sum_e = summarize(results["evolved"]) sum_o = summarize(results["original"]) print("\nSummary (evolved vs original - original grid)") print(f" mean_err: {sum_e['mean_err_mean']:.3e} ± {sum_e['mean_err_std']:.3e} vs {sum_o['mean_err_mean']:.3e} ± {sum_o['mean_err_std']:.3e}") print(f" drift: {sum_e['drift_mean']:.3e} ± {sum_e['drift_std']:.3e} vs {sum_o['drift_mean']:.3e} ± {sum_o['drift_std']:.3e}") print(f" NaN failures: {sum_e['nan_count']} / {nruns} vs {sum_o['nan_count']} / {nruns}") print(f" Avg runtime (s): {sum_e['time_mean']:.3f} vs {sum_o['time_mean']:.3f}") improvement_factor = (sum_o["mean_err_mean"] / sum_e["mean_err_mean"]) if sum_e["mean_err_mean"] > 0 else float("inf") print(f"\nImprovement factor (old_mean_err / new_mean_err) = {improvement_factor:.3f}")

def run_fused_simulation(p: Dict[str, Any], use_evolved: bool, verbose: bool, dx: float) -> Tuple[Dict[str, np.ndarray], float, bool]: rng = rng_default(p["seed"]) Nx = int(p["grid_Nx"]) x = np.linspace(-p["Lx"], p["Lx"], Nx) state = 0.1 + 0.01 * np.exp(-0.5 * (x / 0.2) ** 2) V = p["v0"] * (1.0 + p["alpha_s"] * state) N_expect = 0.02 E_heat = 0.0 Emicro0 = float(N_expect) Emacro0 = energy_of_state(state, dx) E0 = Emicro0 + Emacro0 + E_heat dt = float(p["dt_micro"]) times = np.arange(0.0, float(p["T_total"]), dt) steps = len(times) macrosync = max(1, int(round(p["dt_macro"] / dt))) proj_kernel = np.exp(-0.5 * (x / 0.1) ** 2) proj_kernel /= proj_kernel.sum() coll_kernel = np.exp(-0.5 * (x / p["sigma_fus"]) ** 2) coll_kernel /= coll_kernel.sum() hist = { "t": np.zeros(steps), "N_expect": np.zeros(steps), "E_micro": np.zeros(steps), "E_macro": np.zeros(steps), "E_heat": np.zeros(steps), "Z_proxy": np.zeros(steps), "S_center": np.zeros(steps), "trit_state": np.zeros(steps), "inj_adapt": np.zeros(steps), "curvature": np.zeros(steps), "carbon_strength": np.zeros(steps) } fused_op = create_evolved_best_operator(dx) if use_evolved else None carbon_synth = ZSF_CARBON_SYNTHESIS_grid(dx=dx, debug=args.debug) injection_source = carbon_synth["injection_source"] start_time = time.time() failed = False for k, t in enumerate(times): xi_t = xi_of(t) dN_drive = 1e-2 * xi_t * (1.0 - np.clip(N_expect, 0.0, 1.0)) * dt dN_diss = -p["gammak"] * N_expect * dt noise = rng.normal(scale=math.sqrt(max(p["sigma_noise"], 0.0)) * math.sqrt(dt)) N_expect = max(0.0, N_expect + dN_drive + dN_diss + noise) N_collected = p["kappa_proj"] * N_expect N_residual = max(0.0, N_expect - N_collected) samples = sample_collision_residual(N_residual, p["Nmcsample"], rng) if samples.size > 0: deltainavg = float((p["eta_in"] * samples).mean()) deltaheatavg = float((p["eta_heat"] * samples).mean()) deltafusavg = float((p["eta_fus"] * samples).mean()) cap = p["maxinjectionfrac"] * max(1e-12, N_expect) deltainavg = float(np.clip(deltainavg, -1e-3, min(1e-3, cap))) deltaheatavg = float(np.clip(deltaheatavg, -1e-3, min(1e-3, cap))) deltafusavg = float(np.clip(deltafusavg, -1e-3, min(1e-3, cap))) else: deltainavg = deltaheatavg = deltafusavg = 0.0 Gtot_scalar = N_collected + deltainavg + float(np.mean(injection_source)) Gtot_scalar = min(Gtot_scalar, p["maxGtotscalar"]) G_proj = Gtot_scalar * proj_kernel G_fus = (1.0 - p["eta_fus"]) * deltafusavg * coll_kernel G_total = G_proj + G_fus if (k % macrosync) == 0: if use_evolved and fused_op is not None: state, info = fused_op.apply(state, rng) hist["inj_adapt"][k] = float(info.get("inj_adapt", 0.0)) hist["trit_state"][k] = int(info.get("trit", 0)) else: # Original upwind + diffusion + source F = V * state F_im = np.roll(F, 1) convterm = -(macrosync * dt / dx) * (F - F_im) diffterm = p["nu"] * (macrosync * dt / (dx * dx)) * laplacian_1d(state, dx) srcterm = (macrosync * dt) * (G_total - p["gamma_env"] * state) noise_term = rng.normal(scale=args.sigma_noise if not args.high_noise else args.sigma_noise * 5.0, size=state.shape) state = state + convterm + diffterm + srcterm + noise_term state = 0.995 * state + 0.005 * np.roll(state, 1) if args.cfl_max > 0.0: try: dt = adjust_dt_by_cfl_1d(V, dx, dt, args.cfl_max, 1e-10, args.dt_micro * 10) except Exception: dt = max(1e-10, dt * 0.5) if k % 100 == 0 and k > 0: zs_trigger = {"state": state, "xi": xi_t, "lambda": 0.12, "phi": float(np.mean(state)), "trit_state": int(hist["trit_state"][k]), "entropy": 0.0, "energy": float(hist["E_macro"][k-1])} zs_out, trigger_log = DARKMATTER_TRIGGER_grid(zs_trigger, coupling_section=2e-40, Veff=1.0, rng=rng, dx=dx) state = zs_out["state"] hist["curvature"][k] = float(zs_out.get("curvature", 0.0)) if args.debug: for msg in trigger_log: print(f"[grid trigger] {msg}") V = p["v0"] * (1.0 + p["alpha_s"] * state) maxallowed2 = max(p["maxinjectionfrac"] * max(1e-12, N_expect), p["min_denominator"]) heatincrement = min(deltaheatavg, maxallowed2) newmicromass = min(deltafusavg, maxallowed2) E_heat += heatincrement N_expect = max(0.0, N_expect - heatincrement + newmicromass) Emicro = float(N_expect) Emacro = energy_of_state(state, dx) Etotal = Emicro + Emacro + E_heat Zproxy = abs(Etotal - E0) / max(p["min_denominator"], abs(E0)) if Zproxy > p["Z_threshold"]: deltaE = Etotal - E0 maxabsorb = p.get("maxabsorbenergy_frac", args.maxabsorbenergy_frac) * max(abs(E0), p["min_denominator"]) absorb = math.copysign(min(abs(deltaE), maxabsorb), deltaE) E_heat -= absorb state = state * p["safetydamp_coef"] N_expect = N_expect * p["safetydamp_coef"] hist["t"][k] = t hist["N_expect"][k] = N_expect hist["E_micro"][k] = Emicro hist["E_macro"][k] = Emacro hist["E_heat"][k] = E_heat hist["Z_proxy"][k] = Zproxy hist["S_center"][k] = float(state[len(state)//2]) hist["carbon_strength"][k] = float(carbon_synth["strengthGPa"]) if not np.isfinite(Etotal) or np.isnan(Zproxy) or np.any(~np.isfinite(state)): failed = True break runtime = time.time() - start_time if verbose: print(f"Simulation finished: failed={failed}, runtime={runtime:.3f}s, final_trit={int(hist['trit_state'][-1]) if steps>0 else 0}") return hist, runtime, failed

def run_fused_comparison(nruns: int): seeds = [int(args.seed) + i for i in range(nruns)] results = {"evolved": [], "original": []} for use_evolved in [True, False]: label = "evolved" if use_evolved else "original" for s in seeds: p = { "seed": s, "grid_Nx": args.grid_Nx, "Lx": args.Lx, "T_total": args.T_total, "dt_micro": args.dt_micro, "dt_macro": args.dt_macro, "v0": args.v0, "alpha_s": args.alpha_s, "nu": args.nu, "gamma_env": args.gamma_env, "kappa_proj": args.kappa_proj, "eta_in": 0.05, "eta_heat": 0.10, "eta_fus": 0.05, "sigma_fus": 0.03, "Nmcsample": args.Nmcsample, "gammak": 5e-4, "sigma_noise": args.sigma_noise, "Z_threshold": args.Z_threshold, "safetydamp_coef": args.safetydamp_coef, "maxGtotscalar": args.maxGtotscalar, "min_denominator": args.min_denominator, "maxinjectionfrac": args.maxinjectionfrac, "dt_macro": args.dt_macro } dx = 2.0 * p["Lx"] / max(1, p["grid_Nx"] - 1) hist, runtime, failed = run_fused_simulation(p, use_evolved=use_evolved, verbose=True, dx=dx) E_total = np.array(hist["E_micro"]) + np.array(hist["E_macro"]) + np.array(hist["E_heat"]) E0 = float(E_total[0]) err_E = np.abs(E_total - E0) / max(1e-12, abs(E0)) mean_err = float(np.mean(err_E)) max_err = float(np.max(err_E)) drift = float((E_total[-1] - E0) / max(1e-12, abs(E0))) trit_active_rate = float(np.sum(np.array(hist["trit_state"]) != 0) / len(hist["trit_state"])) results[label].append({ "mean_err": mean_err, "max_err": max_err, "drift": drift, "nan": failed, "runtime": runtime, "trit_active_rate": trit_active_rate, "final_carbon_strength": float(hist["carbon_strength"][-1]) if len(hist["carbon_strength"]) > 0 else 0.0 }) print(f"\n{label} seed={s}:") print(f" mean_err={mean_err:.3e}, max_err={max_err:.3e}, drift={drift:.3e}") print(f" trit_active_rate={trit_active_rate:.2f}, final_carbon_strength={hist['carbon_strength'][-1]:.2f}GPa") print(f" nan={failed}, time={runtime:.3f}s") def summarize(list_of_dicts: List[Dict[str, Any]]) -> Dict[str, float]: mean_errs = [d["mean_err"] for d in list_of_dicts] drifts = [d["drift"] for d in list_of_dicts] trit_rates = [d["trit_active_rate"] for d in list_of_dicts] carb_strs = [d["final_carbon_strength"] for d in list_of_dicts] nans = sum(1 for d in list_of_dicts if d["nan"]) times = [d["runtime"] for d in list_of_dicts] return { "mean_err_mean": mean(mean_errs), "mean_err_std": stdev(mean_errs) if len(mean_errs) > 1 else 0.0, "drift_mean": mean(drifts), "drift_std": stdev(drifts) if len(drifts) > 1 else 0.0, "trit_rate_mean": mean(trit_rates), "carbon_strength_mean": mean(carb_strs), "nan_count": nans, "time_mean": mean(times) } sum_e = summarize(results["evolved"]) sum_o = summarize(results["original"]) print("\n" + "="*50) print("Fused Comparison Summary (Evolved vs Original)") print("="*50) print(f"Mean Energy Error: {sum_e['mean_err_mean']:.3e} ± {sum_e['mean_err_std']:.3e} vs {sum_o['mean_err_mean']:.3e} ± {sum_o['mean_err_std']:.3e}") print(f"Energy Drift: {sum_e['drift_mean']:.3e} ± {sum_e['drift_std']:.3e} vs {sum_o['drift_mean']:.3e} ± {sum_o['drift_std']:.3e}") print(f"Trit Active Rate: {sum_e['trit_rate_mean']:.2f} vs {sum_o['trit_rate_mean']:.2f}") print(f"Final Carbon Strength (GPa): {sum_e['carbon_strength_mean']:.2f} vs {sum_o['carbon_strength_mean']:.2f}") print(f"NaN Failures: {sum_e['nan_count']}/{nruns} vs {sum_o['nan_count']}/{nruns}") print(f"Avg Runtime (s): {sum_e['time_mean']:.3f} vs {sum_o['time_mean']:.3f}") improvement_factor = (sum_o["mean_err_mean"] / sum_e["mean_err_mean"]) if sum_e["mean_err_mean"] > 0 else float("inf") print(f"\nEnergy Error Improvement Factor: {improvement_factor:.3f}x")

def main(): print("=== zsfsynthesizedadaptive_fused: full-model demo (integrated all modules) ===") print(f"injmult={args.inj_mult}, couplingmult={args.coupling_mult}, lowerthreshold={args.lower_threshold}, highnoise={args.high_noise}") sources = [SOURCE_A, SOURCE_B, SOURCE_C] mergeddefault = collidedict(SOURCE_A["params"], SOURCE_B["params"], RNG, exactcancel=args.exact_cancel) mergeddefault = collidedict(mergeddefault, SOURCE_C["params"], RNG, exactcancel=args.exact_cancel) mergeddefault["injectionamp"] = mergeddefault.get("injectionamp", 1e-4) * args.inj_mult if "diffusionnu" in mergeddefault: mergeddefault["diffusionnu"] = mergeddefault.get("diffusionnu", 1e-3) * args.coupling_mult else: mergeddefault["kappa"] = mergeddefault.get("kappa", 1e-3) * args.coupling_mult mergeddefault["stabilize"] = True mergeddefault["stabilize_params"] = { "injectionamp": mergeddefault.get("injectionamp", 1e-4), "nearmargin": 0.12, "mixgain": mergeddefault.get("mixgain", 5.0) * args.inj_mult, "etaamp": mergeddefault.get("etaamp", 1e-5 if not args.high_noise else 5e-5), "holdstrength": 0.12, "thetaplus": 0.35 if args.lower_threshold else 0.5, "thetaminus": -0.35 if args.lower_threshold else -0.5 } # Print sample params samplekeys = ["xi","lambda","injectionamp","etaamp","mixgain","kappa","diffusionnu"] sample = {k: mergeddefault.get(k) for k in samplekeys if k in mergeddefault} print("Sample fused params:") print(json.dumps(sample, indent=2)) # Single-point evolution demo hist = unitpointevolutiondemo(Nsteps=args.steps, seed=12345, paramsoverride=mergeddefault) for i in (0, max(0, args.steps // 4), args.steps - 1): if i < len(hist["phi"]): print(f"step {i}: phi={hist['phi'][i]:.6e}, trit={hist['trit'][i]}, injadapt={hist['injadapt'][i]:.6e}") final_phi = hist["phi"][-1] if len(hist["phi"]) > 0 else None final_trit = hist["trit"][-1] if len(hist["trit"]) > 0 else None print("final phi:", final_phi) print("final tritstate:", final_trit) print() # Trigger demo (single-point) zstrigger = {"xi": 0.021, "lambda": 0.12, "phi": 0.0, "tritstate": 0, "entropy": 0.0, "seed": 54321, "injectionadaptive": 0.0, "injprev": 0.0} evo_params = {} if args.high_noise: evo_params["etaamp"] = 5e-5 if args.lower_threshold: evo_params["thetaplus"] = 0.35 evo_params["thetaminus"] = -0.35 evo_params["stabilize"] = True evo_params["stabilize_params"] = { "injectionamp": evo_params.get("injectionamp", mergeddefault.get("injectionamp", 1e-4)), "nearmargin": 0.12, "mixgain": mergeddefault.get("mixgain", 10.0), "etaamp": evo_params.get("etaamp", mergeddefault.get("etaamp", 1e-5)), "holdstrength": 0.12, "thetaplus": evo_params.get("thetaplus", 0.35), "thetaminus": evo_params.get("thetaminus", -0.35) } zsout, log = DARKMATTERTRIGGER(zstrigger, couplingsection=2e-40, Veff=1.0, rng=rng_default(54321), params={"evoparams": evo_params}) for m in log: print(m) print("after trigger zs:", {k: zsout.get(k) for k in ("phi","tritstate","entropy","energy","curvature")}) print() # Carbon synthesis demo res = ZSFCARBONSYNTHESISdemo(scale=None, clamp=False, autoscale=True, debug=args.debug) print("carbon synthesis result:", res) print() # TRIT distribution stats = unitpointevolutiondemo(Nsteps=2000, seed=42) counts = np.bincount(np.array(stats["trit"]) + 1, minlength=3) dist = counts / len(stats["trit"]) print("TRIT distribution (-1/0/1):", dist) print() # Candidate search (optional) if args.search_candidates > 0: print(f"Searching {args.search_candidates} candidates (grid_scaling=False to preserve single-point baseline)...") cands = generatecandidates([SOURCE_A, SOURCE_B, SOURCE_C], RNG, args.search_candidates, exactcancel=args.exact_cancel, injmult=args.inj_mult, couplingmult=args.coupling_mult, grid_scaling=False) val = validate_candidates(cands, RNG, steps=500, runs=4, debug=args.debug) print("Candidate validation summary (top 6):") for v in val[:6]: print(v) print() # Regression test (optional) rt = None if args.regression_test: print("Running regression test to verify single-point parity...") rt = regression_test() print("Regression result:", rt) print() # Original grid simulation (optional) failed_grid = None if args.use_grid: print("Running original grid microcosm demo...") p = copy.deepcopy(BASE_PARAMS) p.update({"seed": args.seed}) hist_grid, runtime, failed = run_simulation(p, use_evolved=True, verbose=args.debug) failed_grid = failed print("Original grid demo runtime:", runtime, "failed:", failed) print() # Original grid comparison (optional) if args.nruns > 0 and args.use_grid: print("Running original grid comparison harness...") run_comparison(nruns=args.nruns) print() # Fused grid comparison print("=== Fused ZSF-Microcosm Model Comparison Test ===") print(json.dumps({"seed": args.seed, "inj_mult": args.inj_mult, "grid_Nx": args.grid_Nx, "T_total": args.T_total, "nruns": args.nruns}, indent=2)) run_fused_comparison(nruns=args.nruns) # Save a minimal result JSON result_out = { "single_point_final_phi": final_phi, "single_point_final_trit": final_trit, "carbon_strength_GPa": res["strengthGPa"], "regression_test_passed": rt["trit_match"] if args.regression_test else None, "original_grid_failed": failed_grid if args.use_grid else None } with open(args.result_save_path, "w", encoding="utf-8") as f: json.dump(result_out, f, indent=2) print("\n=== Simulation Completed ===") print(f"Result saved to: {args.result_save_path}")

if name == "main": main()