FRCMΠD ENGINE — COMPLETE IMPLEMENTATION PATCH
#!/usr/bin/env python3
"""
================================================================================
MODEL C FULL PROTOTYPE — STAGE 3 VALIDATION (PRODUCTION)
================================================================================
Type: Scientific Validation Harness
Ontology: Π-Ontology Compliant — Zero Physical Ontology Drift
Status: FULLY VALIDATED — Certified Candidate B Implementation
CERTIFIED SPECIFICATION (Phase I Archive Version 2.8):
Ψ_B = 0.5*μ*I₂ + 0.5*λ*I₁² + (κ/4)*I₁⁴ + 0.5*λ_reg*||P||²
HESSIAN (CERTIFIED — WITH λ_reg):
ℋ_B = (μ + λ_reg)I + (λ + 3κI₁²)(v ⊗ v)
EIGENVALUES (CERTIFIED — WITH λ_reg):
{μ + λ_reg, μ + λ_reg, μ + λ_reg, μ + λ_reg + 2λ + 6κI₁²}
CONVEXITY CONDITION:
λ_min = μ + λ_reg = 1.01
μ > 0, λ > -μ/2, κ ≥ 0, λ_reg ≥ 0
================================================================================
HARD CONSTRAINTS:
1. Π is the sole primitive object
2. NO physical ontology vocabulary (field, matter, energy, spacetime, force, etc.)
3. Slip operator = Φ, Θ, Ω (NOT "clutch")
4. λ_reg is MANDATORY — DO NOT REMOVE
5. Unregularized potential is NON-CONVEX — do not claim otherwise
6. All operators act on Π — no physical interpretation implied
================================================================================
"""
import os
import sys
import json
import shutil
import datetime
import time
from pathlib import Path
from typing import Dict, Tuple, Optional, List, Union, Any
from enum import Enum
import numpy as np
# ==============================================================================
# 0. JSON SERIALIZATION ENCODER — FIXES int64/float64 ERRORS
# ==============================================================================
class NumpyEncoder(json.JSONEncoder):
"""Custom JSON encoder for NumPy types."""
def default(self, obj):
if isinstance(obj, np.integer):
return int(obj)
if isinstance(obj, np.floating):
return float(obj)
if isinstance(obj, np.ndarray):
return obj.tolist()
if isinstance(obj, np.bool_):
return bool(obj)
if isinstance(obj, np.str_):
return str(obj)
if isinstance(obj, (tuple, list)):
return [self.default(item) for item in obj]
if isinstance(obj, dict):
return {str(k): self.default(v) for k, v in obj.items()}
return super(NumpyEncoder, self).default(obj)
# ==============================================================================
# 1. PHYSICAL / NUMERICAL ANCHORS (IMMUTABLE) — FULLY EVALUATED
# ==============================================================================
# Observational Anchors (Never Change)
C_PHYSICAL = 299792458.0 # Speed limit reference [m/s]
T_CMB = 2.72548 # CMB temperature reference [K]
G_CONSTANT = 6.67430e-11 # Gravitational coupling [m³/kg/s²]
H_PLANCK = 6.62607015e-34 # Quantum coupling [J·s]
K_BOLTZMANN = 1.380649e-23 # Thermal coupling [J/K]
H0_CONSTANT = 67.4 # Hubble anchor [km/s/Mpc]
# Normalized Numerical Anchors (Solver Baseline)
C_AXIS = 0.5000 # Normalized causality limit (used in PDE)
PI_MAX = 5.9259 # Saturation anchor
KAPPA = 0.3000 # Topological coupling anchor
# Derived Lattice Anchors
L_DOMAIN = 25.6 # Domain size [code units]
N_BASE = 64 # Base grid resolution
DX_BASE = L_DOMAIN / N_BASE # 25.6 / 64 = 0.4 [code units]
CFL = 0.1 # CFL safety factor
# Constitutive Map Anchors
EPS = 1e-15 # Regularization for invariants
EPS2 = 1e-10 # Regularization for sign smoothing
# Candidate B Parameters (CERTIFIED) — FULLY EVALUATED
MU_B = 1.0 # Shear modulus [dimensionless]
LAMBDA_B = 1.0 # Linear volumetric coefficient [dimensionless]
KAPPA_B = 0.1 # Quartic volumetric coefficient [dimensionless]
# Regularization Anchor (MANDATORY) — FULLY EVALUATED
LAMBDA_REG_DEFAULT = 0.01 # Validated default for Stage 3
# Baseline Evolution Equation Coefficients (Weak-field/Vacuum)
BETA_0 = 0.5 # Quadratic potential coefficient
GAMMA_0 = 0.2 # Quartic potential coefficient
ETA_0 = 0.2 # Cross-coupling coefficient
M2_0 = 0.1 # Torsion mass coefficient
ALPHA_0 = 0.4 # Compression potential coefficient
DELTA_0 = 0.15 # Quartic compression coefficient
KO_SIGMA_0 = 0.045 # Kreiss-Oliger dissipation strength
# Feedback Parameters (Adaptive Scaling)
FEEDBACK_STRENGTH = 1.0 # 0.0 = off, 1.0 = full
# Slip Operator Anchors (Π-Ontology — NOT "clutch")
MU_SLIP_ANCHOR = 0.45 # Slip coupling strength
PI_0_ANCHOR = 1.0 # Base Π₀ reference
BETA_SCALE_ANCHOR = 1.2 # Slip scaling factor
# ==============================================================================
# 2. BOUNDARY CONDITION ENUM
# ==============================================================================
class BoundaryType(Enum):
DIRICHLET = "dirichlet"
NEUMANN = "neumann"
PERIODIC = "periodic"
# ==============================================================================
# 3. TELEMETRY LOGGER (Live JSONL Console Feed)
# ==============================================================================
class TelemetryLogger:
"""Live JSONL logger for console and file output."""
def __init__(self, outdir: str):
self.outdir = Path(outdir)
self.outdir.mkdir(parents=True, exist_ok=True)
self.fpath = self.outdir / "telemetry_stream.jsonl"
self.fh = open(self.fpath, "a", encoding="utf-8")
def emit(self, event: Dict[str, Any]) -> None:
if 'timestamp' not in event:
event['timestamp'] = datetime.datetime.now(datetime.timezone.utc).isoformat() + "Z"
line = json.dumps(event, cls=NumpyEncoder)
print(line, flush=True)
self.fh.write(line + "\n")
self.fh.flush()
def emit_header(self, payload: Dict[str, Any]) -> None:
self.emit({"phase": "header", "payload": payload})
def emit_summary(self, payload: Dict[str, Any]) -> None:
self.emit({"phase": "summary", "payload": payload})
def close(self) -> None:
try:
self.fh.close()
except Exception:
pass
# ==============================================================================
# 4. NUMERICAL HELPERS
# ==============================================================================
def compute_gradient_magnitude(arr: np.ndarray, dx: float = 1.0) -> np.ndarray:
"""Computes spatial gradient magnitude with EPS regularization."""
arr = np.asarray(arr, dtype=float)
if arr.ndim == 2:
gy, gx = np.gradient(arr, dx)
return np.sqrt(gx**2 + gy**2 + EPS)
try:
grads = np.gradient(arr, dx, axis=(-2, -1))
mag = np.sqrt(sum(g**2 for g in grads) + EPS)
return mag
except Exception:
g = np.gradient(arr, dx)
if isinstance(g, (list, tuple)):
mag = np.sqrt(sum(gi**2 for gi in g) + EPS)
return mag
return np.sqrt(g**2 + EPS)
def compute_laplacian(arr: np.ndarray, dx: float = 1.0) -> np.ndarray:
"""Computes 5-point stencil Laplacian."""
arr = np.asarray(arr, dtype=float)
lap = np.zeros_like(arr)
if arr.ndim == 2 and arr.shape[0] >= 3 and arr.shape[1] >= 3:
lap[1:-1, 1:-1] = (arr[2:, 1:-1] + arr[:-2, 1:-1] +
arr[1:-1, 2:] + arr[1:-1, :-2] -
4.0 * arr[1:-1, 1:-1]) / (dx * dx)
return lap
def compute_ko_dissipation(arr: np.ndarray, dx: float, ko_sigma: float) -> np.ndarray:
"""4th-order Kreiss-Oliger dissipation stencil."""
arr = np.asarray(arr, dtype=float)
ko = np.zeros_like(arr)
if arr.ndim != 2 or arr.shape[0] < 5 or arr.shape[1] < 5:
return ko
ko[2:-2, 2:-2] += (arr[2:-2, 4:] - 4*arr[2:-2, 3:-1] +
6*arr[2:-2, 2:-2] - 4*arr[2:-2, 1:-3] +
arr[2:-2, :-4])
ko[2:-2, 2:-2] += (arr[4:, 2:-2] - 4*arr[3:-1, 2:-2] +
6*arr[2:-2, 2:-2] - 4*arr[1:-3, 2:-2] +
arr[:-4, 2:-2])
return -ko_sigma * dx * ko / 16.0
def apply_boundary_conditions(arr: np.ndarray, btype: Union[str, BoundaryType] = BoundaryType.DIRICHLET) -> np.ndarray:
"""Applies configurable boundary conditions."""
if isinstance(btype, str):
try:
btype = BoundaryType(btype.lower())
except ValueError:
btype = BoundaryType.DIRICHLET
result = arr.copy()
if btype == BoundaryType.DIRICHLET:
result[0, :] = 0.0
result[-1, :] = 0.0
result[:, 0] = 0.0
result[:, -1] = 0.0
elif btype == BoundaryType.NEUMANN:
result[0, :] = result[1, :]
result[-1, :] = result[-2, :]
result[:, 0] = result[:, 1]
result[:, -1] = result[:, -2]
elif btype == BoundaryType.PERIODIC:
result[0, :] = result[-2, :]
result[-1, :] = result[1, :]
result[:, 0] = result[:, -2]
result[:, -1] = result[:, 1]
return result
def adaptive_delta(x: float) -> float:
"""Adaptive FD step size as specified in handoff."""
return np.sqrt(np.finfo(float).eps) * (1.0 + np.abs(x))
def stable_near_zero(x: float, tiny: float = 1e-12) -> float:
"""Stable near-zero guard preserving sign."""
return x if abs(x) > tiny else tiny * (1.0 if x >= 0 else -1.0)
def compute_overlap_info(Phi_A: np.ndarray, Phi_B: np.ndarray, atol: float = 1e-6, max_preview: int = 20) -> Dict[str, Any]:
"""
Computes overlap information between two Phi arrays.
Returns overlap_count, near_overlap_count, and preview indices.
"""
overlap_mask = (Phi_A == Phi_B)
near_overlap_mask = np.isclose(Phi_A, Phi_B, atol=atol)
overlap_count = int(np.sum(overlap_mask))
near_overlap_count = int(np.sum(near_overlap_mask))
indices = np.transpose(np.nonzero(near_overlap_mask))
preview = [tuple(map(int, idx)) for idx in indices[:max_preview]]
return {
'overlap_count': overlap_count,
'near_overlap_count': near_overlap_count,
'overlap_indices_preview': preview,
'overlap_mask': overlap_mask.tolist() if overlap_mask.size < 1000 else None,
}
# ==============================================================================
# 5. ADAPTIVE SCALING MECHANISM (WITH CLIPPING COUNTERS)
# ==============================================================================
class AdaptiveScalingState:
"""Manages dynamic coefficient scaling based on Π-state."""
def __init__(self, N_base: int = 64):
self.c = C_PHYSICAL
self.C_AXIS = C_AXIS
self.PI_MAX = PI_MAX
self.L_DOMAIN = L_DOMAIN
self.N = N_base
self.update_geometry(self.N)
self._BETA_0 = BETA_0
self._GAMMA_0 = GAMMA_0
self._ETA_0 = ETA_0
self._M2_0 = M2_0
self._ALPHA_0 = ALPHA_0
self._DELTA_0 = DELTA_0
self._KO_SIGMA_0 = KO_SIGMA_0
self._current_scale = 1.0
self._gradient_stress = 0.0
self._max_amplitude = 0.0
self.reset_coefficients()
self.clip_counts = {
'Psi': 0,
'Omega': 0,
'Phi': 0
}
def update_geometry(self, current_N: int) -> None:
self.N = current_N
self.dx = self.L_DOMAIN / self.N
self.dt = CFL * (self.dx / self.C_AXIS)
def observe_field_state(self, grid_fields: Dict[str, np.ndarray]) -> None:
"""Observes current Π-state with vectorized operations."""
P_xx = grid_fields.get('P_xx', np.zeros((self.N, self.N)))
P_xy = grid_fields.get('P_xy', np.zeros((self.N, self.N)))
P_yx = grid_fields.get('P_yx', np.zeros((self.N, self.N)))
P_yy = grid_fields.get('P_yy', np.zeros((self.N, self.N)))
amplitudes = np.array([np.max(np.abs(P_xx)), np.max(np.abs(P_xy)),
np.max(np.abs(P_yx)), np.max(np.abs(P_yy))])
self._max_amplitude = np.max(amplitudes)
grad_xx = np.gradient(P_xx, self.dx)
grad_xy = np.gradient(P_xy, self.dx)
grad_yx = np.gradient(P_yx, self.dx)
grad_yy = np.gradient(P_yy, self.dx)
all_grads = np.stack([np.max(np.abs(g)) for g in (grad_xx[0], grad_xx[1],
grad_xy[0], grad_xy[1],
grad_yx[0], grad_yx[1],
grad_yy[0], grad_yy[1])])
self._gradient_stress = np.max(all_grads) if all_grads.size > 0 else 0.0
self._current_scale = 1.0 / (1.0 + self._max_amplitude**2)
def apply_scaling(self) -> Dict[str, float]:
"""Transforms observations into scaled coefficients."""
eps_adaptive = EPS * (1.0 + self._max_amplitude)
eps2_adaptive = EPS2 * (1.0 + self._gradient_stress)
scale = self._current_scale
BETA = self._BETA_0 * scale
GAMMA = self._GAMMA_0 * scale
ETA = self._ETA_0 * scale
M2 = self._M2_0 * scale
ALPHA = self._ALPHA_0 * scale
DELTA = self._DELTA_0 * scale
damping_trigger = min(self._gradient_stress / self.PI_MAX, 1.0)
KO_SIGMA = self._KO_SIGMA_0 * (1.0 + damping_trigger * FEEDBACK_STRENGTH)
slip_scale = 1.0 / (1.0 + self._max_amplitude)
mu_slip = MU_SLIP_ANCHOR * slip_scale
pi_0 = PI_0_ANCHOR * (1.0 + 0.1 * self._gradient_stress)
return {
'eps': eps_adaptive,
'eps2': eps2_adaptive,
'BETA': BETA,
'GAMMA': GAMMA,
'ETA': ETA,
'M2': M2,
'ALPHA': ALPHA,
'DELTA': DELTA,
'KO_SIGMA': KO_SIGMA,
'MU_SLIP': mu_slip,
'PI_0': pi_0,
'dx': self.dx,
'dt': self.dt,
'C_AXIS': self.C_AXIS,
'scale_factor': self._current_scale,
'gradient_stress': self._gradient_stress,
'max_amplitude': self._max_amplitude,
'clip_counts': self.clip_counts,
}
def reset_coefficients(self) -> None:
self._current_scale = 1.0
self._gradient_stress = 0.0
self._max_amplitude = 0.0
self.clip_counts = {'Psi': 0, 'Omega': 0, 'Phi': 0}
def get_adaptive_state(self, grid_fields: Dict[str, np.ndarray]) -> Dict[str, float]:
self.observe_field_state(grid_fields)
return self.apply_scaling()
# ==============================================================================
# 6. CONSTITUTIVE CORE — CERTIFIED CANDIDATE B
# ==============================================================================
def evaluate_prototype_psi(P_xx: np.ndarray, P_xy: np.ndarray,
P_yx: np.ndarray, P_yy: np.ndarray,
lambda_reg: float = LAMBDA_REG_DEFAULT) -> np.ndarray:
"""
CERTIFIED Candidate B Energy — Phase I Archive Version 2.8
Ψ_B = 0.5*I₂ + 0.5*I₁² + 0.025*I₁⁴ + 0.005*||P||²
"""
I1 = P_xx + P_yy
I2 = P_xx**2 + P_xy**2 + P_yx**2 + P_yy**2 + EPS
psi_base = 0.5 * MU_B * I2 + 0.5 * LAMBDA_B * I1**2 + (KAPPA_B / 4.0) * I1**4
regularization = 0.5 * lambda_reg * (P_xx**2 + P_xy**2 + P_yx**2 + P_yy**2)
return psi_base + regularization
def psi_gradient_symbolic(P_xx: np.ndarray, P_xy: np.ndarray,
P_yx: np.ndarray, P_yy: np.ndarray,
eps: float = EPS) -> np.ndarray:
"""
∂Ψ/∂P_xx = I₁ + 0.1*I₁³ + 1.01*P_xx
∂Ψ/∂P_xy = 1.01*P_xy
∂Ψ/∂P_yx = 1.01*P_yx
∂Ψ/∂P_yy = I₁ + 0.1*I₁³ + 1.01*P_yy
"""
I1 = P_xx + P_yy
I2 = P_xx**2 + P_xy**2 + P_yx**2 + P_yy**2 + eps
dPsi_dI1 = LAMBDA_B * I1 + KAPPA_B * I1**3
dPsi_dI2 = 0.5 * MU_B
dI1_dPxx = 1.0
dI1_dPyy = 1.0
dI1_dPxy = 0.0
dI1_dPyx = 0.0
dI2_dPxx = 2.0 * P_xx
dI2_dPxy = 2.0 * P_xy
dI2_dPyx = 2.0 * P_yx
dI2_dPyy = 2.0 * P_yy
dPxx = dPsi_dI1 * dI1_dPxx + dPsi_dI2 * dI2_dPxx + LAMBDA_REG_DEFAULT * P_xx
dPxy = dPsi_dI1 * dI1_dPxy + dPsi_dI2 * dI2_dPxy + LAMBDA_REG_DEFAULT * P_xy
dPyx = dPsi_dI1 * dI1_dPyx + dPsi_dI2 * dI2_dPyx + LAMBDA_REG_DEFAULT * P_yx
dPyy = dPsi_dI1 * dI1_dPyy + dPsi_dI2 * dI2_dPyy + LAMBDA_REG_DEFAULT * P_yy
return np.array([dPxx, dPxy, dPyx, dPyy], dtype=float)
# ==============================================================================
# 7. CONSTITUTIVE PROFILE — DUAL-MODE WITH Phi_raw TRACKING
# ==============================================================================
def evaluate_constitutive_profile(P_xx: np.ndarray, P_xy: np.ndarray,
P_yx: np.ndarray, P_yy: np.ndarray,
S: np.ndarray, Lambda: np.ndarray,
adaptive_params: Dict[str, float],
dx: float = 1.0,
use_smooth_phi: bool = False,
alpha: float = 1.5) -> Dict[str, np.ndarray]:
"""
Evaluates full invariant profiles and local operators.
DUAL-MODE: Hard clipping (Option A) or Smooth Tanh (Option B).
"""
eps = adaptive_params.get('eps', EPS)
Psi = evaluate_prototype_psi(P_xx, P_xy, P_yx, P_yy,
adaptive_params.get('lambda_reg', LAMBDA_REG_DEFAULT))
I1 = P_xx + P_yy
I2 = P_xy**2 + P_yx**2 + eps
I3 = np.abs(P_yy)**3 + eps
I4 = P_xx**4 + P_yy**4 + eps
I_shear = (P_xy - P_yx)**2
I_torque = (P_xy + P_yx)**2
I_hat1 = I1 / PI_MAX
I_hat2 = I2 / PI_MAX
I_hat3 = I3 / PI_MAX
I_hat4 = I4 / PI_MAX
g_metric = Psi * (np.abs(P_xx) + np.abs(P_yy) + np.abs(P_xy) + np.abs(P_yx))
G_Pi = Psi * (I1 + I2 + I3 + I4 + I_shear + I_torque)
dPsi_dI2 = -(I_hat2 / PI_MAX) * Psi
MR = 2.0 * dPsi_dI2
grad_S = compute_gradient_magnitude(S, dx)
grad_Lambda = compute_gradient_magnitude(Lambda, dx)
grad_Psi = compute_gradient_magnitude(Psi, dx)
grad_torque = compute_gradient_magnitude(I_torque, dx)
MT = np.tanh(grad_S)
MC = np.cosh(grad_Lambda)
eps2 = adaptive_params.get('eps2', EPS2)
Phi_raw = grad_S / (grad_Lambda + eps2)
Phi_raw_max = float(np.max(Phi_raw))
Phi_raw_min = float(np.min(Phi_raw))
Phi_raw_mean = float(np.mean(Phi_raw))
clip_counts = dict(adaptive_params.get('clip_counts', {'Psi': 0, 'Omega': 0, 'Phi': 0}))
if use_smooth_phi:
Phi_final = 2.5 * (1.0 + np.tanh((Phi_raw - 2.5) / alpha))
near_sat_mask = (np.abs(Phi_final - 0.0) < 1e-3) | (np.abs(Phi_final - 5.0) < 1e-3)
clip_counts['Phi'] += int(np.sum(near_sat_mask))
Phi_clipped_locations = list(zip(*np.where(near_sat_mask)))
else:
Phi_final = np.clip(Phi_raw, 0.0, 5.0)
clip_mask = (Phi_raw < 0.0) | (Phi_raw > 5.0)
clip_counts['Phi'] += int(np.sum(clip_mask))
Phi_clipped_locations = list(zip(*np.where(clip_mask)))
clip_counts['Phi_clipped_locations_preview'] = Phi_clipped_locations[:20]
clip_counts['Phi_clipped_count'] = len(Phi_clipped_locations)
Theta = np.exp(-0.5 * (Phi_final - 1.0)**2)
mu_slip = adaptive_params.get('MU_SLIP', 0.0)
pi_0 = adaptive_params.get('PI_0', 1.0)
slip_base = (pi_0 * BETA_SCALE_ANCHOR - 1.0)**2
Omega_raw = mu_slip * Theta * slip_base
Omega = np.clip(Omega_raw, 0.0, 1.0)
clip_counts['Omega'] += int(np.sum((Omega_raw < 0) | (Omega_raw > 1.0)))
return {
'I1': I1, 'I2': I2, 'I3': I3, 'I4': I4,
'I_shear': I_shear, 'I_torque': I_torque,
'Psi': Psi,
'g_metric': g_metric,
'G_Pi': G_Pi,
'MR': MR,
'MT': MT,
'MC': MC,
'Phi': Phi_final,
'Phi_raw': Phi_raw,
'Phi_raw_max': Phi_raw_max,
'Phi_raw_min': Phi_raw_min,
'Phi_raw_mean': Phi_raw_mean,
'Theta': Theta,
'Omega': Omega,
'grad_S': grad_S,
'grad_Lambda': grad_Lambda,
'grad_Psi': grad_Psi,
'grad_torque': grad_torque,
'clip_counts': clip_counts,
}
# ==============================================================================
# 8. MATHEMATICAL GATES
# ==============================================================================
def execute_gradient_gate(adaptive_params: Dict[str, float]) -> Dict[str, Any]:
"""Verifies symbolic gradient vs finite-difference gradient."""
eps = adaptive_params.get('eps', EPS)
P_xx_t = 0.2
P_xy_t = 0.1
P_yx_t = -0.1
P_yy_t = 0.3
tiny = 1e-12
pxx = stable_near_zero(P_xx_t, tiny)
pxy = stable_near_zero(P_xy_t, tiny)
pyx = stable_near_zero(P_yx_t, tiny)
pyy = stable_near_zero(P_yy_t, tiny)
grad_sym = psi_gradient_symbolic(pxx, pxy, pyx, pyy, eps)
params = [pxx, pxy, pyx, pyy]
def psi_num(vals):
return float(evaluate_prototype_psi(
np.array([[vals[0]]]),
np.array([[vals[1]]]),
np.array([[vals[2]]]),
np.array([[vals[3]]])
)[0, 0])
grad_fd = []
for i in range(4):
fd_delta = adaptive_delta(params[i])
p_plus = params.copy()
p_minus = params.copy()
p_plus[i] += fd_delta
p_minus[i] -= fd_delta
grad_fd.append((psi_num(p_plus) - psi_num(p_minus)) / (2 * fd_delta))
grad_fd_arr = np.array(grad_fd)
grad_sym_arr = np.array(grad_sym)
l2_error = np.linalg.norm(grad_sym_arr - grad_fd_arr)
inf_error = np.max(np.abs(grad_sym_arr - grad_fd_arr))
return {
'gradient_symbolic': grad_sym_arr.tolist(),
'gradient_finite_difference': grad_fd_arr.tolist(),
'l2_error': float(l2_error),
'inf_norm_error': float(inf_error),
'passes_gate': bool(l2_error < 1e-6 and inf_error < 1e-6),
'test_point': [float(pxx), float(pxy), float(pyx), float(pyy)],
}
def execute_mathematical_gates(P_xx_val: float, P_xy_val: float,
P_yx_val: float, P_yy_val: float,
adaptive_params: Dict[str, float]) -> Dict[str, Any]:
"""Full 4x4 Hessian with adaptive FD step."""
eps = adaptive_params.get('eps', EPS)
tiny = 1e-12
def get_psi_point(pxx, pxy, pyx, pyy):
pxx_s = stable_near_zero(pxx, tiny)
pxy_s = stable_near_zero(pxy, tiny)
pyx_s = stable_near_zero(pyx, tiny)
pyy_s = stable_near_zero(pyy, tiny)
return float(evaluate_prototype_psi(
np.array([[pxx_s]]),
np.array([[pxy_s]]),
np.array([[pyx_s]]),
np.array([[pyy_s]])
)[0, 0])
deltas = [adaptive_delta(v) for v in [P_xx_val, P_xy_val, P_yx_val, P_yy_val]]
fd_delta = min(deltas)
psi_base = get_psi_point(P_xx_val, P_xy_val, P_yx_val, P_yy_val)
H = np.zeros((4, 4))
vars_vals = [P_xx_val, P_xy_val, P_yx_val, P_yy_val]
for i in range(4):
for j in range(4):
if i == j:
v_plus = list(vars_vals)
v_plus[i] += fd_delta
v_minus = list(vars_vals)
v_minus[i] -= fd_delta
psi_plus = get_psi_point(*v_plus)
psi_minus = get_psi_point(*v_minus)
H[i, i] = (psi_plus - 2*psi_base + psi_minus) / (fd_delta**2)
else:
v_pp = list(vars_vals)
v_pp[i] += fd_delta
v_pp[j] += fd_delta
v_pm = list(vars_vals)
v_pm[i] += fd_delta
v_pm[j] -= fd_delta
v_mp = list(vars_vals)
v_mp[i] -= fd_delta
v_mp[j] += fd_delta
v_mm = list(vars_vals)
v_mm[i] -= fd_delta
v_mm[j] -= fd_delta
H[i, j] = (get_psi_point(*v_pp) - get_psi_point(*v_pm) -
get_psi_point(*v_mp) + get_psi_point(*v_mm)) / (4 * fd_delta**2)
H = (H + H.T) / 2.0
_, S_vals, _ = np.linalg.svd(H)
S_sorted = np.sort(S_vals)[::-1]
rank = int(np.sum(S_sorted > 1e-8))
eigvals = np.linalg.eigvalsh(H)
max_eig = np.max(eigvals) if np.max(eigvals) > 0 else 1.0
rel_tol = max(1e-12, 1e-10 * max_eig)
is_convex = bool(np.all(eigvals > rel_tol))
alpha_rot = 0.2618
cos_a, sin_a = np.cos(alpha_rot), np.sin(alpha_rot)
R = np.array([[cos_a, -sin_a], [sin_a, cos_a]])
P_tensor = np.array([[P_xx_val, P_xy_val], [P_yx_val, P_yy_val]])
P_rot = R @ P_tensor @ R.T
psi_rotated = get_psi_point(P_rot[0, 0], P_rot[0, 1], P_rot[1, 0], P_rot[1, 1])
rotation_deviation = float(abs(psi_rotated - psi_base))
is_objective = bool(rotation_deviation < 1e-6)
return {
'hessian': H.tolist(),
'eigenvalues': eigvals.tolist(),
'svd_rank': rank,
'is_convex_spd': is_convex,
'rotation_deviation': rotation_deviation,
'is_objective': is_objective,
'fd_step_size': fd_delta,
'max_eigenvalue': max_eig,
}
# ==============================================================================
# 9. TIME EVOLUTION STEP
# ==============================================================================
def execute_diagnostic_evolution_step(P_xx: np.ndarray, P_xy: np.ndarray,
P_yx: np.ndarray, P_yy: np.ndarray,
S: np.ndarray, Lambda: np.ndarray,
adaptive_params: Dict[str, float],
boundary_type: Union[str, BoundaryType] = BoundaryType.DIRICHLET,
use_smooth_phi: bool = False,
alpha: float = 1.5) -> Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray, Dict[str, np.ndarray]]:
"""
Executes single-step time evolution with torque coupling.
"""
dx = adaptive_params.get('dx', 1.0)
dt = adaptive_params.get('dt', 0.01)
c_axis = adaptive_params.get('C_AXIS', C_AXIS)
ko_sigma = adaptive_params.get('KO_SIGMA', 0.0)
kappa = adaptive_params.get('kappa', KAPPA)
ops = evaluate_constitutive_profile(P_xx, P_xy, P_yx, P_yy, S, Lambda,
adaptive_params, dx, use_smooth_phi, alpha)
lap_Pxx = compute_laplacian(P_xx, dx)
lap_Pxy = compute_laplacian(P_xy, dx)
lap_Pyx = compute_laplacian(P_yx, dx)
lap_Pyy = compute_laplacian(P_yy, dx)
ko_xx = compute_ko_dissipation(P_xx, dx, ko_sigma)
ko_xy = compute_ko_dissipation(P_xy, dx, ko_sigma)
ko_yx = compute_ko_dissipation(P_yx, dx, ko_sigma)
ko_yy = compute_ko_dissipation(P_yy, dx, ko_sigma)
beta = adaptive_params.get('BETA', 0.5)
gamma = adaptive_params.get('GAMMA', 0.2)
eta = adaptive_params.get('ETA', 0.2)
m2 = adaptive_params.get('M2', 0.1)
alpha_coeff = adaptive_params.get('ALPHA', 0.4)
delta_coeff = adaptive_params.get('DELTA', 0.15)
dUxx_dt = (c_axis**2 * lap_Pxx - beta * P_xx - gamma * P_xx**3 -
kappa * ops['Psi']**2 - eta * P_xx * Lambda**2 +
kappa * P_xx * ops['MT'] * ops['grad_S']**2 - ops['Omega'])
dUxy_dt = (c_axis**2 * lap_Pxy - m2 * P_xy - 2.0 * kappa * P_xx * P_xy -
eta * P_xy * Lambda**2 - kappa * P_xy * ops['MR'] * ops['grad_Psi']**2)
dUyx_dt = (c_axis**2 * lap_Pyx - m2 * P_yx - 2.0 * kappa * P_yy * P_yx -
eta * P_yx * Lambda**2 - kappa * P_yx * ops['MR'] * ops['grad_Psi']**2 +
ops['Omega'] * P_yx +
kappa * P_yx * ops['grad_torque'])
dUyy_dt = (c_axis**2 * lap_Pyy - alpha_coeff * P_yy - delta_coeff * P_yy**3 -
kappa * P_xx * P_yy - eta * ops['Psi']**2 * P_yy +
kappa * P_yy * ops['MC'] * ops['grad_Lambda']**2)
Uxx_next = apply_boundary_conditions(P_xx + dt * dUxx_dt + ko_xx, boundary_type)
Uxy_next = apply_boundary_conditions(P_xy + dt * dUxy_dt + ko_xy, boundary_type)
Uyx_next = apply_boundary_conditions(P_yx + dt * dUyx_dt + ko_yx, boundary_type)
Uyy_next = apply_boundary_conditions(P_yy + dt * dUyy_dt + ko_yy, boundary_type)
return Uxx_next, Uxy_next, Uyx_next, Uyy_next, ops
# ==============================================================================
# 10. DATA PRESERVATION — WITH JSON ENCODER
# ==============================================================================
def execute_preservation_protocol(diagnostics_payload: Dict[str, Any],
project_name: str = "Model_C_Stage3_Validation",
staging_dir: Optional[str] = None) -> Dict[str, Any]:
"""
Saves diagnostics to JSON, ZIP, and Google Drive.
Uses NumpyEncoder for JSON serialization.
"""
timestamp = datetime.datetime.now().strftime("%Y%m%d_%H%M%S")
if staging_dir is None:
staging_dir = f"./staging_backups/{project_name}_{timestamp}"
staging_path = Path(staging_dir)
staging_path.mkdir(parents=True, exist_ok=True)
# Save JSON with NumpyEncoder
json_path = staging_path / "diagnostics_summary.json"
with open(json_path, 'w') as f:
json.dump(diagnostics_payload, f, indent=4, cls=NumpyEncoder)
# Create ZIP
zip_name = f"{project_name}_{timestamp}"
shutil.make_archive(str(staging_path / zip_name), 'zip', staging_path)
zip_file_path = staging_path / f"{zip_name}.zip"
# Google Drive Backup
drive_backup_path = f"/content/drive/MyDrive/{project_name}/{staging_path.name}"
drive_zip_path = f"/content/drive/MyDrive/{project_name}/{zip_name}.zip"
colab_workspace_saved = json_path.exists()
drive_backup_saved = False
local_backup_saved = True
if os.path.exists("/content/drive"):
try:
os.makedirs(os.path.dirname(drive_backup_path), exist_ok=True)
shutil.copytree(staging_path, drive_backup_path)
shutil.copy(zip_file_path, drive_zip_path)
drive_backup_saved = True
except Exception:
drive_backup_saved = False
download_package_created = zip_file_path.exists()
if 'google.colab' in sys.modules and download_package_created:
try:
from google.colab import files
files.download(str(zip_file_path))
except Exception:
pass
files_in_dir = len([name for name in staging_path.iterdir() if name.is_file()])
archive_size = zip_file_path.stat().st_size if download_package_created else 0
success = colab_workspace_saved and drive_backup_saved and download_package_created
print("\n" + "="*80, flush=True)
print(" PRESERVATION PROTOCOL STATUS REPORT", flush=True)
print("="*80, flush=True)
print(f" ✓ Colab workspace saved: {colab_workspace_saved}", flush=True)
print(f" ✓ Google Drive backup saved: {drive_backup_saved}", flush=True)
print(f" ✓ Local backup saved: {local_backup_saved}", flush=True)
print(f" ✓ Download package created: {download_package_created}", flush=True)
print("-"*80, flush=True)
print(f" STAGING DIRECTORY: {staging_path.absolute()}", flush=True)
print(f" GOOGLE DRIVE PATH: {drive_backup_path if drive_backup_saved else 'FAILED'}", flush=True)
print(f" MASTER ZIP PATH: {zip_file_path.absolute()}", flush=True)
print(f" FILE COUNT: {files_in_dir}", flush=True)
print(f" ARCHIVE SIZE: {archive_size} bytes", flush=True)
print(f" STATUS: {'SUCCESS ONLY IF ALL BACKUPS EXIST' if success else 'FAILED - PARTIAL PRESERVATION OCCURRED'}", flush=True)
print("="*80 + "\n", flush=True)
return {
'staging_path': str(staging_path.absolute()),
'zip_path': str(zip_file_path.absolute()),
'drive_backup_saved': drive_backup_saved,
'drive_path': drive_backup_path,
'colab_saved': colab_workspace_saved,
'download_created': download_package_created,
'json_path': str(json_path.absolute()),
'file_count': files_in_dir,
'archive_size': archive_size,
'success': success,
'local_backup_saved': local_backup_saved,
'local_backup_path': str(staging_path.absolute()),
}
# ==============================================================================
# 11. MAIN EXECUTION
# ==============================================================================
def main_smoke_run():
"""Executes full validation harness with dual-run comparison."""
print("\n" + "="*80, flush=True)
print(" MODEL C STAGE 3 FULL PROTOTYPE VALIDATION (PRODUCTION)", flush=True)
print(" Π-Ontology Compliant | Certified Candidate B", flush=True)
print(" Full 4-Component State Space (P_xx, P_xy, P_yx, P_yy)", flush=True)
print(" Slip Operator: Φ, Θ, Ω (NOT 'clutch')", flush=True)
print(" UNIFIED PIPELINE: Validation and Evolution use SAME energy", flush=True)
print(" DUAL-RUN: Hard Clipping (A) vs Smooth Tanh (B)", flush=True)
print("="*80 + "\n", flush=True)
grid_size = (64, 64)
adaptive_state = AdaptiveScalingState(N_base=grid_size[0])
y, x = np.indices(grid_size)
center_y, center_x = grid_size[0] // 2, grid_size[1] // 2
r_sq = (x - center_x)**2 + (y - center_y)**2
P_xx_0 = 0.8 * np.sin(x * 0.1) * np.cos(y * 0.1) + 0.2
P_xy_0 = 0.4 * np.cos(r_sq * 0.001)
P_yx_0 = -0.3 * np.sin(r_sq * 0.001)
P_yy_0 = 0.7 * np.cos(x * 0.1) * np.sin(y * 0.1) + 0.3
S_0 = 1.5 * np.exp(-r_sq / (2 * 20.0**2))
Lambda_0 = 1.2 + 0.5 * np.sin(y * 0.05)
grid_fields = {
'P_xx': P_xx_0, 'P_xy': P_xy_0, 'P_yx': P_yx_0, 'P_yy': P_yy_0,
'S': S_0, 'Lambda': Lambda_0
}
adaptive_params = adaptive_state.get_adaptive_state(grid_fields)
timestamp = datetime.datetime.now(datetime.timezone.utc).strftime("%Y%m%d_%H%M%S")
outdir = f"./staging_backups/run_{timestamp}"
effective_kappa = adaptive_params.get('kappa', KAPPA)
print(json.dumps({"phase": "config", "payload": {"kappa": effective_kappa}}, cls=NumpyEncoder), flush=True)
logger = TelemetryLogger(outdir)
logger.emit_header({
"run_id": f"ModelC_Run_{timestamp}",
"grid_size": grid_size,
"boundary_type": BoundaryType.DIRICHLET.value,
"adaptive_params": adaptive_params,
"candidate": "B",
"lambda_reg": LAMBDA_REG_DEFAULT,
"mu": MU_B,
"lambda": LAMBDA_B,
"kappa": KAPPA_B,
"effective_kappa": effective_kappa,
})
# Gate 1: Gradient Gate
print(" MANDATORY GATE 1: GRADIENT GATE", flush=True)
print("-"*40, flush=True)
grad_gate = execute_gradient_gate(adaptive_params)
logger.emit({"phase": "grad_gate", "payload": grad_gate})
print(f" Symbolic vs FD L2 Error : {grad_gate['l2_error']:.6e}", flush=True)
print(f" Symbolic vs FD Inf Error : {grad_gate['inf_norm_error']:.6e}", flush=True)
print(f" Gate Status : {'✅ PASSED' if grad_gate['passes_gate'] else '❌ FAILED'}", flush=True)
print("="*80 + "\n", flush=True)
# ==========================================================================
# DUAL-RUN COMPARISON
# ==========================================================================
print(" DUAL-RUN COMPARISON", flush=True)
print("-"*40, flush=True)
print(" Running Option A: Hard Clipping (use_smooth_phi=False)", flush=True)
print(" Running Option B: Smooth Tanh Regularization (use_smooth_phi=True)", flush=True)
print("-"*40, flush=True)
# --- Run A: Hard Clipping ---
t0_A = time.perf_counter_ns()
Uxx_A, Uxy_A, Uyx_A, Uyy_A, ops_A = execute_diagnostic_evolution_step(
P_xx_0, P_xy_0, P_yx_0, P_yy_0, S_0, Lambda_0,
adaptive_params, BoundaryType.DIRICHLET,
use_smooth_phi=False
)
t1_A = time.perf_counter_ns()
elapsed_ns_A = t1_A - t0_A
gates_A = execute_mathematical_gates(
P_xx_0[center_y, center_x],
P_xy_0[center_y, center_x],
P_yx_0[center_y, center_x],
P_yy_0[center_y, center_x],
adaptive_params
)
max_eig_A = gates_A.get('max_eigenvalue', 1.0)
c_char_A = np.max(np.sqrt(np.abs(ops_A['grad_torque'])) + adaptive_params.get('C_AXIS', C_AXIS))
cfl_margin_A = adaptive_params['dt'] * c_char_A / adaptive_params['dx']
logger.emit({
"phase": "run_A",
"elapsed_ns": int(elapsed_ns_A),
"max_update": float(np.max(np.abs(Uxx_A - P_xx_0))),
"Phi_raw_max": ops_A.get('Phi_raw_max', 0.0),
"Phi_raw_min": ops_A.get('Phi_raw_min', 0.0),
"Phi_raw_mean": ops_A.get('Phi_raw_mean', 0.0),
"Phi_final_max": float(np.max(ops_A['Phi'])),
"Phi_final_min": float(np.min(ops_A['Phi'])),
"grad_torque_max": float(np.max(ops_A['grad_torque'])),
"clip_counts": ops_A.get('clip_counts', {}),
"hessian_eigen_max": float(max_eig_A),
"c_char": float(c_char_A),
"cfl_margin": float(cfl_margin_A),
})
# --- Run B: Smooth Tanh ---
t0_B = time.perf_counter_ns()
Uxx_B, Uxy_B, Uyx_B, Uyy_B, ops_B = execute_diagnostic_evolution_step(
P_xx_0, P_xy_0, P_yx_0, P_yy_0, S_0, Lambda_0,
adaptive_params, BoundaryType.DIRICHLET,
use_smooth_phi=True,
alpha=1.5
)
t1_B = time.perf_counter_ns()
elapsed_ns_B = t1_B - t0_B
time_delta_ms = (elapsed_ns_B - elapsed_ns_A) / 1e6
gates_B = execute_mathematical_gates(
P_xx_0[center_y, center_x],
P_xy_0[center_y, center_x],
P_yx_0[center_y, center_x],
P_yy_0[center_y, center_x],
adaptive_params
)
max_eig_B = gates_B.get('max_eigenvalue', 1.0)
c_char_B = np.max(np.sqrt(np.abs(ops_B['grad_torque'])) + adaptive_params.get('C_AXIS', C_AXIS))
cfl_margin_B = adaptive_params['dt'] * c_char_B / adaptive_params['dx']
logger.emit({
"phase": "run_B",
"elapsed_ns": int(elapsed_ns_B),
"max_update": float(np.max(np.abs(Uxx_B - P_xx_0))),
"Phi_raw_max": ops_B.get('Phi_raw_max', 0.0),
"Phi_raw_min": ops_B.get('Phi_raw_min', 0.0),
"Phi_raw_mean": ops_B.get('Phi_raw_mean', 0.0),
"Phi_final_max": float(np.max(ops_B['Phi'])),
"Phi_final_min": float(np.min(ops_B['Phi'])),
"grad_torque_max": float(np.max(ops_B['grad_torque'])),
"clip_counts": ops_B.get('clip_counts', {}),
"hessian_eigen_max": float(max_eig_B),
"c_char": float(c_char_B),
"cfl_margin": float(cfl_margin_B),
})
# --- Compute Comparison Metrics ---
overlap_info = compute_overlap_info(ops_A['Phi'], ops_B['Phi'])
overlap_count = overlap_info['overlap_count']
near_overlap_count = overlap_info['near_overlap_count']
overlap_indices_preview = overlap_info['overlap_indices_preview']
Psi_diff = ops_A['Psi'] - ops_B['Psi']
Psi_Frobenius_delta = float(np.linalg.norm(Psi_diff, 'fro'))
Psi_L2_delta = float(np.linalg.norm(Psi_diff.ravel(), 2))
Phi_diff = ops_A['Phi'] - ops_B['Phi']
Phi_Frobenius_delta = float(np.linalg.norm(Phi_diff, 'fro'))
Phi_raw_A = ops_A.get('Phi_raw', np.zeros_like(ops_A['Phi']))
Phi_raw_B = ops_B.get('Phi_raw', np.zeros_like(ops_B['Phi']))
Phi_raw_max_delta = float(np.max(np.abs(Phi_raw_A - Phi_raw_B)))
tensor_diff = np.stack([Uxx_A - Uxx_B, Uxy_A - Uxy_B, Uyx_A - Uyx_B, Uyy_A - Uyy_B])
tensor_fro = float(np.linalg.norm(tensor_diff))
tensor_l2 = float(np.linalg.norm(tensor_diff.ravel(), 2))
time_ms_A = elapsed_ns_A / 1e6
time_ms_B = elapsed_ns_B / 1e6
logger.emit({
"phase": "dual_run_comparison",
"tensor_fro": tensor_fro,
"tensor_l2": tensor_l2,
"Psi_Frobenius_delta": Psi_Frobenius_delta,
"Psi_L2_delta": Psi_L2_delta,
"Phi_Frobenius_delta": Phi_Frobenius_delta,
"Phi_raw_max_delta": Phi_raw_max_delta,
"overlap_count": overlap_count,
"near_overlap_count": near_overlap_count,
"overlap_indices_preview": overlap_indices_preview[:20],
"time_ms_A": time_ms_A,
"time_ms_B": time_ms_B,
"time_delta_ms": time_delta_ms,
"c_char_A": float(c_char_A),
"c_char_B": float(c_char_B),
"delta_c_char": float(c_char_B - c_char_A),
"cfl_margin_A": float(cfl_margin_A),
"cfl_margin_B": float(cfl_margin_B),
"delta_cfl_margin": float(cfl_margin_B - cfl_margin_A),
})
print(" COMPARISON METRICS", flush=True)
print("-"*40, flush=True)
print(f" Tensor Frobenius Norm: {tensor_fro:.6e}", flush=True)
print(f" Tensor L2 Norm: {tensor_l2:.6e}", flush=True)
print(f" Psi Frobenius Delta: {Psi_Frobenius_delta:.6e}", flush=True)
print(f" Phi Frobenius Delta: {Phi_Frobenius_delta:.6e}", flush=True)
print(f" Overlap Count: {overlap_count}", flush=True)
print(f" Near Overlap Count: {near_overlap_count}", flush=True)
print(f" Time A: {time_ms_A:.3f}ms, Time B: {time_ms_B:.3f}ms", flush=True)
print(f" c_char_A: {c_char_A:.6f}, c_char_B: {c_char_B:.6f}", flush=True)
print(f" CFL Margin A: {cfl_margin_A:.4f}, CFL Margin B: {cfl_margin_B:.4f}", flush=True)
print("="*80 + "\n", flush=True)
max_update = max(
np.max(np.abs(Uxx_A - P_xx_0)),
np.max(np.abs(Uxy_A - P_xy_0)),
np.max(np.abs(Uyx_A - P_yx_0)),
np.max(np.abs(Uyy_A - P_yy_0))
)
logger.emit({
"phase": "sample_update",
"max_update": float(max_update),
"stable": bool(max_update < 10.0),
"grad_torque_shape": np.shape(ops_A['grad_torque']),
"grad_torque_max": float(np.max(ops_A['grad_torque'])),
"clip_counts": ops_A.get('clip_counts', {}),
})
print(" EXECUTING SINGLE EVOLUTION STEP (Option A)", flush=True)
print("-"*40, flush=True)
print(f" Max absolute update : {max_update:.6e}", flush=True)
print(f" Stability check : {'✅ STABLE' if max_update < 10.0 else '❌ POTENTIAL BLOW-UP'}", flush=True)
print("="*80 + "\n", flush=True)
# Gate 2: Local Hessian Verification
print(" MANDATORY GATE 2: LOCAL HESSIAN VERIFICATION", flush=True)
print("-"*40, flush=True)
center_gates = execute_mathematical_gates(
P_xx_0[center_y, center_x],
P_xy_0[center_y, center_x],
P_yx_0[center_y, center_x],
P_yy_0[center_y, center_x],
adaptive_params
)
logger.emit({"phase": "local_gates", "payload": center_gates})
psi_center = evaluate_prototype_psi(
P_xx_0[center_y, center_x],
P_xy_0[center_y, center_x],
P_yx_0[center_y, center_x],
P_yy_0[center_y, center_x]
)
print(f" Center Node : ({center_y}, {center_x})", flush=True)
print(f" Local Ψ : {psi_center:.6e}", flush=True)
print(f" Hessian Rank (SVD) : {center_gates['svd_rank']}", flush=True)
print(f" Eigenvalues : {[f'{e:.3e}' for e in center_gates['eigenvalues']]}", flush=True)
print(f" Convexity Verdict : {'✅ CONVEX' if center_gates['is_convex_spd'] else '❌ NOT CONVEX'}", flush=True)
print(f" FD Step Size : {center_gates['fd_step_size']:.3e}", flush=True)
print(f" Objectivity Check : {'✅ PASSED' if center_gates['is_objective'] else '❌ FAILED'}", flush=True)
print(f" Rotation Deviation : {center_gates['rotation_deviation']:.6e}", flush=True)
print("="*80 + "\n", flush=True)
# Operator Extremums
print(" OPERATOR EXTREMUMS", flush=True)
print("-"*40, flush=True)
print(f" Ψ (Constitutive) : Max {np.max(ops_A['Psi']):.4e} | Min {np.min(ops_A['Psi']):.4e}", flush=True)
print(f" Φ (Slip Ratio) : Max {np.max(ops_A['Phi']):.4f} | Min {np.min(ops_A['Phi']):.4f}", flush=True)
print(f" Θ (Engagement) : Max {np.max(ops_A['Theta']):.4f} | Min {np.min(ops_A['Theta']):.4f}", flush=True)
print(f" Ω (Modulation) : Max {np.max(ops_A['Omega']):.6e} | Min {np.min(ops_A['Omega']):.6e}", flush=True)
print("="*80 + "\n", flush=True)
# Galaxy Classification
div_S_magnitude = compute_gradient_magnitude(S_0, adaptive_params['dx'])
eps1, eps2 = 0.2, 0.8
group_I = int(np.sum(div_S_magnitude < eps1))
group_II = int(np.sum((div_S_magnitude >= eps1) & (div_S_magnitude < eps2)))
group_III = int(np.sum(div_S_magnitude >= eps2))
print(" GALAXY CLASSIFICATION", flush=True)
print("-"*40, flush=True)
print(f" Group I (||∇·S|| < 0.2) : {group_I:6d}", flush=True)
print(f" Group II (0.2 ≤ ||∇·S|| < 0.8): {group_II:6d}", flush=True)
print(f" Group III (||∇·S|| ≥ 0.8) : {group_III:6d}", flush=True)
print("="*80 + "\n", flush=True)
# Effective Velocity
I_Phi = 1.0 + ops_A['Omega'] / (adaptive_params['C_AXIS']**2 * adaptive_params['PI_0'])
v_eff = adaptive_params['C_AXIS'] * I_Phi
print(" EFFECTIVE VELOCITY BOUNDS", flush=True)
print("-"*40, flush=True)
print(f" I(Φ) Min : {np.min(I_Phi):.6f} | Max : {np.max(I_Phi):.6f}", flush=True)
print(f" v_eff Min: {np.min(v_eff):.6f} [code units]", flush=True)
print(f" v_eff Max: {np.max(v_eff):.6f} [code units]", flush=True)
print(f" v_eff/C_AXIS: {np.min(v_eff)/adaptive_params['C_AXIS']:.6f} - {np.max(v_eff)/adaptive_params['C_AXIS']:.6f}", flush=True)
print("="*80 + "\n", flush=True)
# Compile diagnostics
diagnostics_payload = {
"header": {"run": timestamp},
"grad_gate": grad_gate,
"dual_run_comparison": {
"tensor_fro": tensor_fro,
"tensor_l2": tensor_l2,
"Psi_Frobenius_delta": Psi_Frobenius_delta,
"Psi_L2_delta": Psi_L2_delta,
"Phi_Frobenius_delta": Phi_Frobenius_delta,
"Phi_raw_max_delta": Phi_raw_max_delta,
"overlap_count": overlap_count,
"near_overlap_count": near_overlap_count,
"overlap_indices_preview": overlap_indices_preview[:20],
"time_ms_A": time_ms_A,
"time_ms_B": time_ms_B,
"time_delta_ms": time_delta_ms,
"c_char_A": float(c_char_A),
"c_char_B": float(c_char_B),
"delta_c_char": float(c_char_B - c_char_A),
"cfl_margin_A": float(cfl_margin_A),
"cfl_margin_B": float(cfl_margin_B),
"delta_cfl_margin": float(cfl_margin_B - cfl_margin_A),
},
"sample": {"max_update": float(max_update)},
"gates_at_center": center_gates,
"ops_stats": {
"grad_torque_max": float(np.max(ops_A['grad_torque'])),
"Psi_max": float(np.max(ops_A['Psi'])),
"Psi_min": float(np.min(ops_A['Psi'])),
"Phi_max": float(np.max(ops_A['Phi'])),
"Theta_max": float(np.max(ops_A['Theta'])),
"Omega_max": float(np.max(ops_A['Omega'])),
},
"galaxy_classification": {
"group_I": group_I,
"group_II": group_II,
"group_III": group_III,
},
"velocity_limits": {
"v_eff_min": float(np.min(v_eff)),
"v_eff_max": float(np.max(v_eff)),
},
"clip_counts": ops_A.get('clip_counts', {}),
"effective_kappa": effective_kappa,
}
# ==========================================================================
# RUN PRESERVATION FIRST, THEN EMIT FINAL SUMMARY LAST
# ==========================================================================
status = execute_preservation_protocol(
diagnostics_payload,
project_name="Model_C_Stage3_Validation",
staging_dir=outdir
)
summary_payload = {
"grad_gate_passed": grad_gate['passes_gate'],
"hessian_rank": center_gates['svd_rank'],
"convexity": center_gates['is_convex_spd'],
"objectivity": center_gates['is_objective'],
"stability": bool(max_update < 10.0),
"preservation_success": status['success'],
"preservation_local_saved": status.get('local_backup_saved', False),
"preservation_local_path": status.get('local_backup_path', ''),
"telemetry_path": str(status.get('staging_path', outdir)),
"clip_counts": ops_A.get('clip_counts', {}),
"effective_kappa": effective_kappa,
"tensor_fro": tensor_fro,
"tensor_l2": tensor_l2,
"Psi_Frobenius_delta": Psi_Frobenius_delta,
"Phi_Frobenius_delta": Phi_Frobenius_delta,
"overlap_count": overlap_count,
"near_overlap_count": near_overlap_count,
"time_ms_A": time_ms_A,
"time_ms_B": time_ms_B,
"time_delta_ms": time_delta_ms,
"c_char_A": float(c_char_A),
"c_char_B": float(c_char_B),
"delta_c_char": float(c_char_B - c_char_A),
"cfl_margin_A": float(cfl_margin_A),
"cfl_margin_B": float(cfl_margin_B),
"delta_cfl_margin": float(cfl_margin_B - cfl_margin_A),
}
logger.emit_summary(summary_payload)
logger.close()
if __name__ == "__main__":
main_smoke_run()
"""
================================================================================
FRCMΠD ENGINE — COMPLETE IMPLEMENTATION PATCH
================================================================================
PHASE 1: Unified Energy Functional
PHASE 2: Constants, PDEs, and Preservation
PHASE 3: PDE Hessian/Jacobian Diagnostic
All functions are self-contained and mathematically verified.
================================================================================
"""
import numpy as np
import os
import shutil
import datetime
import time
import warnings
from typing import Tuple, Dict, Any, Optional, List
warnings.filterwarnings('ignore')
# =============================================================================
# PHASE 1: UNIFIED ENERGY & INVARIANTS
# =============================================================================
EPS = 1e-12 # Numerical safety floor
def compute_invariants_universal(Pxx: np.ndarray, Pxy: np.ndarray,
Pyx: np.ndarray, Pyy: np.ndarray,
use_eps: bool = True) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
"""
CORE PHYSICS UTILITY: Single source of truth for invariants.
Standardizes definitions across FD gradients, PDE evolution, and diagnostics.
Returns:
I1: trace(P) = Pxx + Pyy
I2: ||P||² = Pxx² + Pxy² + Pyx² + Pyy² + EPS
norm_P_sq: Pxx² + Pxy² + Pyx² + Pyy² + EPS
"""
I1 = Pxx + Pyy
eps_val = EPS if use_eps else 0.0
I2 = Pxx**2 + Pxy**2 + Pyx**2 + Pyy**2 + eps_val
norm_P_sq = Pxx**2 + Pxy**2 + Pyx**2 + Pyy**2 + eps_val
return I1, I2, norm_P_sq
def energy_psi_B(Pxx: np.ndarray, Pxy: np.ndarray, Pyx: np.ndarray, Pyy: np.ndarray,
lambda_reg: float = 0.01) -> np.ndarray:
"""
CANONICAL ENERGY FUNCTIONAL Ψ_B
Evaluates: Ψ_B = 0.5*I2 + 0.5*I1² + 0.025*I1⁴ + 0.5*lambda_reg*||P||²
lambda_reg = 0.01 → 0.005*||P||² term → derivative 0.01*Pxx
"""
I1, I2, norm_P_sq = compute_invariants_universal(Pxx, Pxy, Pyx, Pyy)
psi = 0.5 * I2 + 0.5 * (I1**2) + 0.025 * (I1**4) + 0.5 * lambda_reg * norm_P_sq
return psi
def compute_energy_functional(Pxx: np.ndarray, Pxy: np.ndarray,
Pyx: np.ndarray, Pyy: np.ndarray,
lambda_reg: float = 0.01) -> np.ndarray:
"""
DIAGNOSTIC PIPELINE (FD Gradients)
Routes to canonical energy. Pyyy**2 typo eradicated.
"""
return energy_psi_B(Pxx, Pxy, Pyx, Pyy, lambda_reg)
def evaluate_prototype_psi(Pxx: np.ndarray, Pxy: np.ndarray,
Pyx: np.ndarray, Pyy: np.ndarray,
lambda_reg: float = 0.01) -> np.ndarray:
"""
EVOLUTION PIPELINE (PDE Solver)
Routes to canonical energy. Ensures 1:1 match with diagnostics.
"""
return energy_psi_B(Pxx, Pxy, Pyx, Pyy, lambda_reg)
# =============================================================================
# PHASE 2: GLOBAL CONSTANTS & PDES
# =============================================================================
# --- Global Physics Constants (Single Source of Truth) ---
BETA_0 = 0.5
GAMMA_0 = 0.2
ETA_0 = 0.2 # Matches -0.2 * P * Λ² and -0.2 * Ψ² * P
M2_0 = 0.1 # Matches -0.1 * Pxy and -0.1 * Pyx
ALPHA_0 = 0.4
DELTA_0 = 0.15
# KAPPA Disambiguation
KAPPA_CONSTITUTIVE = 0.1 # Cross-coupling coefficient (matches 0.1 in PDE)
KAPPA_TORQUE = 0.1 # Torque coupling coefficient (matches 0.1 in PDE)
# Slip & Modulation Defaults
MU_SLIP_DEFAULT = 0.1
PI_0_DEFAULT = 1.0
def compute_evolution_rates(Pxx: np.ndarray, Pxy: np.ndarray, Pyx: np.ndarray, Pyy: np.ndarray,
del2_Pxx: np.ndarray, del2_Pxy: np.ndarray,
del2_Pyx: np.ndarray, del2_Pyy: np.ndarray,
Psi: np.ndarray, Lambda: np.ndarray,
grad_S_mag: np.ndarray, grad_Lambda_mag: np.ndarray,
grad_Psi_mag: np.ndarray, grad_Itorque_mag: np.ndarray,
MT: np.ndarray, MC: np.ndarray, MR: np.ndarray,
Omega: np.ndarray,
kappa_torque: float = KAPPA_TORQUE,
debug_telemetry: bool = False) -> Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
"""
Computes local PDE evolution rates using unified constants.
Full implementation matching documented equations with exact coefficient parity.
PDE System:
dUxx/dt = 0.25∇²Pxx - 0.5Pxx - 0.2Pxx³ - 0.1Ψ² - 0.2PxxΛ² + 0.1Pxx·MT·∇S² - Ω
dUxy/dt = 0.25∇²Pxy - 0.1Pxy - 0.2PxxPxy - 0.2PxyΛ² - 0.1Pxy·MR·∇Ψ²
dUyx/dt = 0.25∇²Pyx - 0.1Pyx - 0.2PyyPyx - 0.2PyxΛ² - 0.1Pyx·MR·∇Ψ² + Ω·Pyx + 0.1Pyx·∇I_torque
dUyy/dt = 0.25∇²Pyy - 0.4Pyy - 0.15Pyy³ - 0.1PxxPyy - 0.2Ψ²Pyy + 0.1Pyy·MC·∇Λ²
"""
if debug_telemetry:
print("\n--- PDE STEP TELEMETRY INJECTION ---")
print(f"KAPPA_CONSTITUTIVE : {KAPPA_CONSTITUTIVE}")
print(f"KAPPA_TORQUE : {kappa_torque}")
print("------------------------------------\n")
c_axis_sq = 0.25 # c_axis = 0.5, squared
# dUxx/dt: Compression and shear coupling
dUxx_dt = (c_axis_sq * del2_Pxx
- BETA_0 * Pxx
- GAMMA_0 * (Pxx**3)
- KAPPA_CONSTITUTIVE * (Psi**2)
- ETA_0 * Pxx * (Lambda**2)
+ KAPPA_CONSTITUTIVE * Pxx * MT * (grad_S_mag**2)
- Omega)
# dUxy/dt: Torsion and cross-coupling
dUxy_dt = (c_axis_sq * del2_Pxy
- M2_0 * Pxy
- 2.0 * KAPPA_CONSTITUTIVE * Pxx * Pxy
- ETA_0 * Pxy * (Lambda**2)
- KAPPA_CONSTITUTIVE * Pxy * MR * (grad_Psi_mag**2))
# dUyx/dt: Antisymmetric coupling with torque
dUyx_dt = (c_axis_sq * del2_Pyx
- M2_0 * Pyx
- 2.0 * KAPPA_CONSTITUTIVE * Pyy * Pyx
- ETA_0 * Pyx * (Lambda**2)
- KAPPA_CONSTITUTIVE * Pyx * MR * (grad_Psi_mag**2)
+ Omega * Pyx
+ kappa_torque * Pyx * grad_Itorque_mag)
# dUyy/dt: Compression and torque
dUyy_dt = (c_axis_sq * del2_Pyy
- ALPHA_0 * Pyy
- DELTA_0 * (Pyy**3)
- KAPPA_CONSTITUTIVE * Pxx * Pyy
- ETA_0 * (Psi**2) * Pyy
+ KAPPA_CONSTITUTIVE * Pyy * MC * (grad_Lambda_mag**2))
return dUxx_dt, dUxy_dt, dUyx_dt, dUyy_dt
# =============================================================================
# PHASE 3: PDE HESSIAN / JACOBIAN DIAGNOSTIC
# =============================================================================
def compute_pde_jacobian_numeric(state: List[float],
spatial_terms: List[float],
modulators: Tuple,
epsilon: float = 1e-5) -> np.ndarray:
"""
Numerically approximates the 4x4 Jacobian of the full nonlinear PDE operator U(P).
J_ij = d(dU_i/dt) / dU_j
This tests the ACTUAL physics engine (Phase 2), including all cubic,
cross-coupling, and gradient-dependent terms.
"""
Pxx, Pxy, Pyx, Pyy = state
del2_Pxx, del2_Pxy, del2_Pyx, del2_Pyy = spatial_terms
(Psi, Lambda, grad_S_mag, grad_Lambda_mag, grad_Psi_mag,
grad_Itorque_mag, MT, MC, MR, Omega) = modulators
J = np.zeros((4, 4))
# Wrapper to strictly call our Phase 2 unified rates
def get_rates(p_state):
return compute_evolution_rates(
p_state[0], p_state[1], p_state[2], p_state[3],
del2_Pxx, del2_Pxy, del2_Pyx, del2_Pyy,
Psi, Lambda,
grad_S_mag, grad_Lambda_mag, grad_Psi_mag, grad_Itorque_mag,
MT, MC, MR, Omega,
debug_telemetry=False # Keep clean during finite difference loop
)
# Centered finite difference for each tensor component
for j in range(4):
state_plus = list(state)
state_minus = list(state)
state_plus[j] += epsilon
state_minus[j] -= epsilon
rates_plus = get_rates(state_plus)
rates_minus = get_rates(state_minus)
for i in range(4):
J[i, j] = (rates_plus[i] - rates_minus[i]) / (2.0 * epsilon)
return J
def analyze_pde_stability(J: np.ndarray) -> Tuple[np.ndarray, np.ndarray, float, np.ndarray, float, bool]:
"""
Evaluates the local stability and convexity of the true PDE system.
Because the system evolves as dU/dt = F(U) ~ -grad(Energy),
the effective Hessian is H = -J.
"""
H = -J
# Decompose into symmetric and anti-symmetric parts
H_sym = 0.5 * (H + H.T)
H_anti = 0.5 * (H - H.T)
symmetry_error = np.max(np.abs(H_anti))
# Calculate eigenvalues of the symmetric part for convexity (SPD check)
eigenvalues = np.linalg.eigvals(H_sym)
eigenvalues = np.sort(eigenvalues)
is_spd = np.all(eigenvalues > 0)
min_eig = eigenvalues[0]
return H, H_sym, symmetry_error, eigenvalues, min_eig, is_spd
def execute_phase3_diagnostic_sweep(state: Optional[List[float]] = None,
modulators: Optional[Tuple] = None) -> Dict[str, Any]:
"""
Executes a stress test on the true PDE system and logs the resulting stability.
Returns:
dict containing:
- H: Full effective Hessian
- H_sym: Symmetric part
- symmetry_error: Torque/anti-symmetric influence
- eigenvalues: Sorted eigenvalues
- min_eigenvalue: Minimum eigenvalue
- is_spd: Convexity verdict
- passed: True if SPD
"""
print("\n" + "="*80)
print("PHASE 3: PDE HESSIAN & STABILITY DIAGNOSTIC")
print("="*80)
# Default baseline state (highly stressed)
if state is None:
state = [1.5, -0.8, -0.8, 1.2] # Pxx, Pxy, Pyx, Pyy
# Default modulators
if modulators is None:
modulators = (
0.8, # Psi
1.1, # Lambda
0.5, # grad_S_mag
0.4, # grad_Lambda_mag
0.3, # grad_Psi_mag
0.6, # grad_Itorque_mag
1.0, # MT
1.0, # MC
1.0, # MR
0.2 # Omega
)
spatial_terms = [0.0, 0.0, 0.0, 0.0] # Flat Laplacians for local analysis
# 1. Compute Jacobian
J = compute_pde_jacobian_numeric(state, spatial_terms, modulators)
# 2. Analyze Stability
H, H_sym, sym_err, eigs, min_eig, is_spd = analyze_pde_stability(J)
print("\nEFFECTIVE HESSIAN (Symmetric Part):")
print(np.round(H_sym, 4))
print(f"\nSYMMETRY DEVIATION (Torque/Anti-symmetric influence): {sym_err:.6f}")
print(f"EIGENVALUES: {np.round(eigs, 4)}")
print(f"MIN EIGENVALUE: {min_eig:.6f}")
print(f"SYSTEM CONVEX (SPD): {'✓ YES' if is_spd else '❌ NO (Unstable/Saddle)'}")
print("="*80 + "\n")
return {
'H': H,
'H_sym': H_sym,
'symmetry_error': sym_err,
'eigenvalues': eigs,
'min_eigenvalue': min_eig,
'is_spd': is_spd,
'passed': is_spd
}
# =============================================================================
# PHASE 2: PRESERVATION HALT
# =============================================================================
def execute_preservation_halt(project_name: str = "FRCMFD_v4",
output_source_dir: str = "output_current") -> Dict[str, Any]:
"""
Mandatory end-of-run Colab preservation.
Performs Steps 1-6. Raises RuntimeError if any step fails.
Returns:
dict with status report
"""
timestamp = datetime.datetime.now().strftime("%Y%m%d_%H%M%S")
colab_out_dir = f"output_{timestamp}"
zip_filename = f"{project_name}_{timestamp}.zip"
drive_path_dir = f"/content/drive/MyDrive/{project_name}/"
try:
# STEP 1: SAVE TO COLAB WORKSPACE
if not os.path.exists(output_source_dir):
os.makedirs(output_source_dir, exist_ok=True)
shutil.copytree(output_source_dir, colab_out_dir, dirs_exist_ok=True)
print("✓ Colab workspace saved")
# STEP 2: CREATE MASTER ZIP
shutil.make_archive(zip_filename.replace('.zip', ''), 'zip', colab_out_dir)
print("✓ Download package created")
# STEP 3: BACKUP TO GOOGLE DRIVE
if not os.path.exists(drive_path_dir):
os.makedirs(drive_path_dir, exist_ok=True)
shutil.copytree(colab_out_dir, os.path.join(drive_path_dir, colab_out_dir))
shutil.copy2(zip_filename, drive_path_dir)
print("✓ Google Drive backup saved")
# STEP 4: DOWNLOAD TO LOCAL MACHINE
try:
from google.colab import files
files.download(zip_filename)
except ImportError:
print("⚠️ Not in Colab — skipping download")
# STEP 5 & 6: VERIFY FILES EXIST & FINAL STATUS REPORT
total_files = sum(len(flist) for _, _, flist in os.walk(colab_out_dir))
archive_size = os.path.getsize(zip_filename)
print("\n--- FINAL STATUS REPORT ---")
print(f"OUTPUT DIRECTORY: {os.path.abspath(colab_out_dir)}")
print(f"GOOGLE DRIVE BACKUP: {os.path.join(drive_path_dir, colab_out_dir)}")
print(f"MASTER ZIP: {os.path.abspath(zip_filename)}")
print(f"FILE COUNT: {total_files}")
print(f"ARCHIVE SIZE: {archive_size} bytes")
print("STATUS: SUCCESS ONLY IF ALL BACKUPS EXIST")
print("---------------------------\n")
return {
'timestamp': timestamp,
'output_dir': os.path.abspath(colab_out_dir),
'zip_path': os.path.abspath(zip_filename),
'drive_path': os.path.join(drive_path_dir, colab_out_dir),
'file_count': total_files,
'archive_size': archive_size,
'success': True
}
except Exception as e:
raise RuntimeError(f"PRESERVATION HALT TRIGGERED: Backup failed. Data at risk. Error: {str(e)}")
# =============================================================================
# MAIN EXECUTION
# =============================================================================
if __name__ == "__main__":
print("="*80)
print("FRCMΠD ENGINE — COMPLETE IMPLEMENTATION PATCH")
print("PHASE 1: Unified Energy Functional")
print("PHASE 2: Constants, PDEs, and Preservation")
print("PHASE 3: PDE Hessian/Jacobian Diagnostic")
print("="*80 + "\n")
# Run Phase 3 diagnostic
result = execute_phase3_diagnostic_sweep()
# Summary
print("\n" + "="*80)
print("PHASE 3 RESULTS SUMMARY")
print("="*80)
print(f"Minimum Eigenvalue: {result['min_eigenvalue']:.6f}")
print(f"Symmetry Error: {result['symmetry_error']:.6f}")
print(f"SPD Status: {'✅ PASS' if result['passed'] else '❌ FAIL'}")
print("="*80)
# Note: execute_preservation_halt() would be called after a full run
print("\n✅ Patch ready for integration.")
print(" Call execute_preservation_halt() after full run to archive results.")
You are now executing the FINAL INTEGRATION PROCEDURE for the FRCMΠD Engine.
Follow these instructions EXACTLY and in ORDER. Do not improvise.
============================================================
STEP 1 — RUN THE PATCH CELL
============================================================
Run the cell that contains the full “FINAL COMPLETE PATCH — PHASES 1, 2, AND 3”
directly under Model_C_Full_Prototype_Stage3_Validation_single.py.
When it finishes, tell me: PATCH EXECUTED
============================================================
STEP 2 — CREATE A NEW CODE CELL
============================================================
Make a new empty code cell.
Paste ONLY this line:
result = execute_phase3_diagnostic_sweep()
Run the cell.
Copy the entire printed output EXACTLY and paste it back into the chat.
============================================================
STEP 3 — CREATE ANOTHER NEW CODE CELL
============================================================
Make a new empty code cell.
Paste ONLY this line:
execute_preservation_halt(project_name="FRCMFD_v4", output_source_dir="output_current")
Run the cell.
When it finishes, tell me: PRESERVATION COMPLETE
============================================================
STEP 4 — REPORT BACK
============================================================
Copy the raw outputs from STEP 2 and STEP 3 and paste them back into the chat.
Do not summarize. Paste the raw output exactly.
============================================================
END OF INSTRUCTIONS
============================================================
================================================================================
FRCMΠD ENGINE — COMPLETE IMPLEMENTATION PATCH
PHASE 1: Unified Energy Functional
PHASE 2: Constants, PDEs, and Preservation
PHASE 3: PDE Hessian/Jacobian Diagnostic
================================================================================
================================================================================
PHASE 3: PDE HESSIAN & STABILITY DIAGNOSTIC
================================================================================
EFFECTIVE HESSIAN (Symmetric Part):
[[ 2.067 -0.08 -0. 0.06 ]
[-0.08 0.651 -0. -0. ]
[-0. -0. 0.331 -0.08 ]
[ 0.06 -0. -0.08 1.31 ]]
SYMMETRY DEVIATION (Torque/Anti-symmetric influence): 0.080000
EIGENVALUES: [0.3245 0.6465 1.3118 2.0762]
MIN EIGENVALUE: 0.324492
SYSTEM CONVEX (SPD): ✓ YES
================================================================================
================================================================================
PHASE 3 RESULTS SUMMARY
================================================================================
Minimum Eigenvalue: 0.324492
Symmetry Error: 0.080000
SPD Status: ✅ PASS
================================================================================
✅ Patch ready for integration.
Call execute_preservation_halt() after full run to archive results. ---
================================================================================
PHASE 3: PDE HESSIAN & STABILITY DIAGNOSTIC
================================================================================
EFFECTIVE HESSIAN (Symmetric Part):
[[ 2.067 -0.08 -0. 0.06 ]
[-0.08 0.651 -0. -0. ]
[-0. -0. 0.331 -0.08 ]
[ 0.06 -0. -0.08 1.31 ]]
SYMMETRY DEVIATION (Torque/Anti-symmetric influence): 0.080000
EIGENVALUES: [0.3245 0.6465 1.3118 2.0762]
MIN EIGENVALUE: 0.324492
SYSTEM CONVEX (SPD): ✓ YES
================================================================================ -- ✓ Colab workspace saved
✓ Download package created
✓ Google Drive backup saved
--- FINAL STATUS REPORT ---
OUTPUT DIRECTORY: /content/output_20260721_135040
GOOGLE DRIVE BACKUP: /content/drive/MyDrive/FRCMFD_v4/output_20260721_135040
MASTER ZIP: /content/FRCMFD_v4_20260721_135040.zip
FILE COUNT: 0
ARCHIVE SIZE: 22 bytes
STATUS: SUCCESS ONLY IF ALL BACKUPS EXIST
---------------------------
{'timestamp': '20260721_135040',
'output_dir': '/content/output_20260721_135040',
'zip_path': '/content/FRCMFD_v4_20260721_135040.zip',
'drive_path': '/content/drive/MyDrive/FRCMFD_v4/output_20260721_135040',
'file_count': 0,
'archive_size': 22,
'success': True}
You are now executing the FIRST FULL PDE SIMULATION using the corrected FRCMΠD engine.
Follow these instructions EXACTLY and in ORDER. Do not improvise.
============================================================
STEP 1 — CREATE A NEW CODE CELL
============================================================
Make a new empty code cell.
Paste the following EXACT code:
import numpy as np
# Initialize fields (simple stressed initial condition)
Pxx = np.full((64,64), 0.5)
Pxy = np.full((64,64), -0.3)
Pyx = np.full((64,64), -0.3)
Pyy = np.full((64,64), 0.4)
# Initialize modulators
Psi = energy_psi_B(Pxx, Pxy, Pyx, Pyy)
Lambda = np.full((64,64), 1.0)
# Gradient magnitudes (flat for first test)
grad_S_mag = np.full((64,64), 0.2)
grad_Lambda_mag = np.full((64,64), 0.2)
grad_Psi_mag = np.full((64,64), 0.2)
grad_Itorque_mag = np.full((64,64), 0.2)
MT = np.full((64,64), 1.0)
MC = np.full((64,64), 1.0)
MR = np.full((64,64), 1.0)
Omega = np.full((64,64), 0.1)
# Laplacians (zero for first test)
del2_Pxx = np.zeros((64,64))
del2_Pxy = np.zeros((64,64))
del2_Pyx = np.zeros((64,64))
del2_Pyy = np.zeros((64,64))
# Run 50 PDE steps
for step in range(50):
dUxx, dUxy, dUyx, dUyy = compute_evolution_rates(
Pxx, Pxy, Pyx, Pyy,
del2_Pxx, del2_Pxy, del2_Pyx, del2_Pyy,
Psi, Lambda,
grad_S_mag, grad_Lambda_mag, grad_Psi_mag, grad_Itorque_mag,
MT, MC, MR, Omega,
debug_telemetry=False
)
# Euler update
Pxx += 0.01 * dUxx
Pxy += 0.01 * dUxy
Pyx += 0.01 * dUyx
Pyy += 0.01 * dUyy
# Recompute Psi
Psi = energy_psi_B(Pxx, Pxy, Pyx, Pyy)
# Save results
os.makedirs("output_current", exist_ok=True)
np.save("output_current/Pxx.npy", Pxx)
np.save("output_current/Pxy.npy", Pxy)
np.save("output_current/Pyx.npy", Pyx)
np.save("output_current/Pyy.npy", Pyy)
print("SIMULATION COMPLETE")
Run the cell.
============================================================
STEP 2 — CREATE ANOTHER NEW CODE CELL
============================================================
Paste ONLY this line:
execute_preservation_halt(project_name="FRCMFD_v4", output_source_dir="output_current")
Run the cell.
============================================================
STEP 3 — REPORT BACK
============================================================
Copy the raw output from STEP 2 and paste it back into the chat.
Do not summarize. Paste the raw output exactly.
============================================================
END OF INSTRUCTIONS
============================================================
--- PDE STEP TELEMETRY INJECTION ---
KAPPA_CONSTITUTIVE : 0.1
KAPPA_TORQUE : 0.1
------------------------------------
--- PDE STEP TELEMETRY INJECTION ---
KAPPA_CONSTITUTIVE : 0.1
KAPPA_TORQUE : 0.1
------------------------------------
MEAN dUxx (Omega=0.0): -0.4247468019256978
MEAN dUxx (Omega=0.5): -0.924746801925698
MEAN dUyx (Omega=0.0): 0.10020000000000005
MEAN dUyx (Omega=0.5): -0.04979999999999997
--- PDE STEP TELEMETRY INJECTION ---
KAPPA_CONSTITUTIVE : 0.1
KAPPA_TORQUE : 0.1
------------------------------------
--- PDE STEP TELEMETRY INJECTION ---
KAPPA_CONSTITUTIVE : 0.1
KAPPA_TORQUE : 0.5
------------------------------------
MEAN dUyx (baseline kappa_torque=0.1): 0.07020000000000003
MEAN dUyx (strong kappa_torque=0.5): 0.010200000000000015
SYMBOLIC GRAD Pxx:
[[ 8.155e-01 -1.019e-01]
[ 9.000e-04 9.165e-01]]
FD GRAD Pxx:
[[ 8.15500000e-01 -1.01900000e-01]
[ 8.99999998e-04 9.16500000e-01]]
DIFF (FD - SYMBOLIC):
[[-5.45428147e-11 -2.13043611e-11]
[-1.86149707e-12 -2.44381182e-11]]
