mirror of
https://github.com/allaunthefox/Research-Stack.git
synced 2026-07-31 03:05:21 +00:00
92bc00c4d5
- PistSimulation.lean: proven goldenContractionEnergyDecrease (no sorry) 7 supporting lemmas, h_u'_nonneg + h_pt hypothesis, fold induction - Connectors.lean: restored zeroIsVoid theorem with Q16_16 proof - CanonSerialization.lean: removed dead theorem, documented blocker - FixedPointBridge.lean: eliminated Float from compute paths PIST predictions pipeline: - pist_matrix_builder.py: reproducible matrix-only builder (SHA256) - build_pist_matrices_278.py: generates PIST/Matrices278.lean - PIST/Classify.lean: classifyProxy/classifyExact stubs (v2 surface) - PIST/Matrices278.lean: 250-entry matrix HashMap - build_corpus278.py: reads predictions artifact, uses classify* - Pipeline contract documented in root AGENTS.md Cleanup: - Archived 5 orphan pist_* shims, 5 old route_repair variants - Quarantined PIST/Repair.lean (no external callers) - Created 4 opencode agents for remaining TODO items Build: PistSimulation 3309, Compiler 3313, Full 3571 (0 errors)
2026 lines
91 KiB
Text
2026 lines
91 KiB
Text
/- Copyright (c) 2026 Sovereign Research Stack. All rights reserved.
|
||
Released under Apache 2.0 license as described in the file LICENSE.
|
||
Authors: Research Stack Team
|
||
|
||
PistSimulation.lean — PIST Data Slice Processing Pipeline
|
||
|
||
This module models the functional data transformations from the PIST
|
||
interactive simulation (Injection → Pruning → Convergence).
|
||
|
||
Pipeline phases:
|
||
1. Injection — Load raw geometric states into active tensor set
|
||
2. Predictive Pruning — Hardware predictor kills ~95% of doomed paths
|
||
3. Blitter & Gossip — Discrete Picard integral + local gossip clustering
|
||
|
||
Maps directly to WebGPU compute shader dispatch:
|
||
• Phase 1: VRAM initialization
|
||
• Phase 2: Predictor kernel (early out)
|
||
• Phase 3: Blitter physics kernel + Gossip reduction
|
||
|
||
Per AGENTS.md §0: Lean is the source of truth.
|
||
Per AGENTS.md §1.4: Uses Q16_16 (Fix16) throughout.
|
||
-/
|
||
|
||
import Semantics.FixedPoint
|
||
import Semantics.ShellModel
|
||
import Semantics.SSMS
|
||
|
||
namespace Semantics.PistSimulation
|
||
|
||
open Semantics
|
||
open Semantics.ShellModel
|
||
open Semantics.SSMS
|
||
|
||
-- ════════════════════════════════════════════════════════════
|
||
-- §1 Tensor Data Structure
|
||
-- ════════════════════════════════════════════════════════════
|
||
|
||
/-- Single particle/data point in the PIST visual simulation.
|
||
Represents a candidate state in the (a,b) perfect-square coordinate space.
|
||
|
||
Fields:
|
||
• id — Unique identifier for tracking
|
||
• a — Distance from lower perfect square (k²)
|
||
• b — Distance to upper perfect square ((k+1)²)
|
||
• confidence — Gossip-accumulated viability score
|
||
• isActive — Survival flag (false = pruned/dimmed) -/
|
||
structure TensorData where
|
||
id : Nat
|
||
a : Q16_16
|
||
b : Q16_16
|
||
confidence : Q16_16
|
||
isActive : Bool
|
||
deriving Repr, DecidableEq, Inhabited
|
||
|
||
/-- Zero tensor (inactive, zero confidence). -/
|
||
def TensorData.zero (id : Nat) : TensorData :=
|
||
{ id := id, a := Q16_16.zero, b := Q16_16.zero,
|
||
confidence := Q16_16.zero, isActive := false }
|
||
|
||
|
||
-- ════════════════════════════════════════════════════════════
|
||
-- §2 Phase 1: Injection (Canvas Population)
|
||
-- ════════════════════════════════════════════════════════════
|
||
|
||
/-- Maps to visual step where points populate the screen.
|
||
Loads raw (a,b,confidence) tuples into active TensorData array.
|
||
|
||
In WebGPU execution: this initializes VRAM with geometric states.
|
||
Each tensor maps to one workgroup thread's initial state. -/
|
||
def injectDataSlice (rawInputs : Array (Q16_16 × Q16_16 × Q16_16)) : Array TensorData :=
|
||
rawInputs.mapIdx (λ i val =>
|
||
{ id := i,
|
||
a := val.1,
|
||
b := val.2.1,
|
||
confidence := val.2.2,
|
||
isActive := true })
|
||
|
||
/-- Alternative injection from shell state indices.
|
||
Converts event indices to (a,b) coordinates for PIST simulation. -/
|
||
def injectFromShellStates (indices : List Nat) : Array TensorData :=
|
||
let coords := indices.map (λ n =>
|
||
let s := shellState n
|
||
(Q16_16.ofInt (Int.ofNat s.a),
|
||
Q16_16.ofInt (Int.ofNat s.b),
|
||
Q16_16.ofInt (Int.ofNat (s.a * s.b))))
|
||
injectDataSlice coords.toArray
|
||
|
||
|
||
-- ════════════════════════════════════════════════════════════
|
||
-- §3 Phase 2: Predictive Pruning (Heuristic Guillotine)
|
||
-- ════════════════════════════════════════════════════════════
|
||
|
||
/-- Hardware viability predictor.
|
||
Evaluates fast geometric heuristic to kill doomed paths early.
|
||
|
||
PIST criterion: |a - b| > threshold indicates far from perfect square.
|
||
Near-perfect-squares have a ≈ b (symmetric position in shell).
|
||
|
||
Returns true if particle survives pruning. -/
|
||
def predictViability (a b confidence : Q16_16) : Bool :=
|
||
let diff := Q16_16.abs (Q16_16.sub a b)
|
||
let threshold := Q16_16.ofInt 2 -- Within 2 units of symmetry
|
||
let confThreshold := Q16_16.div (Q16_16.ofInt 1) (Q16_16.ofInt 10) -- 0.1 confidence minimum
|
||
Q16_16.lt diff threshold && Q16_16.gt confidence confThreshold
|
||
|
||
/-- Phase 2: Apply predictive pruning to entire dataset.
|
||
Maps to visual step where ~95% of points turn dim and stop.
|
||
|
||
In WebGPU: This is a compute kernel with early-out for pruned threads. -/
|
||
def phase2Pruning (dataset : Array TensorData) : Array TensorData :=
|
||
dataset.map (λ pt =>
|
||
if pt.isActive then
|
||
let viable := predictViability pt.a pt.b pt.confidence
|
||
-- If not viable: particle "turns red and fades out"
|
||
{ pt with isActive := viable,
|
||
confidence := if viable then pt.confidence else Q16_16.zero }
|
||
else pt)
|
||
|
||
|
||
-- ════════════════════════════════════════════════════════════
|
||
-- §4 Phase 3: Blitter & Gossip (Convergence)
|
||
-- ════════════════════════════════════════════════════════════
|
||
|
||
/-- Discrete Picard Integral (Blitter) step.
|
||
Model 131 ODE: F(a,b,ε) = (1 + ε(0.5b + 0.3), -1 + ε(0.5a - 0.3))
|
||
|
||
Performs one timestep of O(1) discrete integration:
|
||
a' = a + ε · (1 + 0.5·b + 0.3)
|
||
b' = b + ε · (-1 + 0.5·a - 0.3)
|
||
|
||
Maps to WGSL: `blit_result = blit_op(fa, fb, timestep_mask)` -/
|
||
def picardBlitStep (a b epsilon : Q16_16) : Q16_16 × Q16_16 :=
|
||
let half := Q16_16.div (Q16_16.ofInt 1) (Q16_16.ofInt 2)
|
||
let c3 := Q16_16.div (Q16_16.ofInt 3) (Q16_16.ofInt 10)
|
||
let fa := Q16_16.add (Q16_16.ofInt 1)
|
||
(Q16_16.mul epsilon (Q16_16.add (Q16_16.mul half b) c3))
|
||
let fb := Q16_16.add (Q16_16.ofInt (-1))
|
||
(Q16_16.mul epsilon (Q16_16.sub (Q16_16.mul half a) c3))
|
||
let nextA := Q16_16.add a (Q16_16.mul epsilon fa)
|
||
let nextB := Q16_16.add b (Q16_16.mul epsilon fb)
|
||
(nextA, nextB)
|
||
|
||
/-- Local gossip confidence aggregation.
|
||
Simulates neighbor-to-neighbor confidence sharing in workgroup.
|
||
|
||
In WebGPU: This uses shared memory / LDS for neighbor access.
|
||
Returns updated confidence from local neighborhood average. -/
|
||
def localGossip (neighbors : Array Q16_16) (selfConfidence : Q16_16) : Q16_16 :=
|
||
if neighbors.size = 0 then selfConfidence
|
||
else
|
||
let sum := neighbors.foldl (λ acc c => Q16_16.add acc c) Q16_16.zero
|
||
let avg := Q16_16.div sum (Q16_16.ofInt (Int.ofNat neighbors.size))
|
||
-- Weighted mix: 70% self + 30% neighbor average
|
||
let mixed := Q16_16.add (Q16_16.mul (Q16_16.div (Q16_16.ofInt 7) (Q16_16.ofInt 10)) selfConfidence)
|
||
(Q16_16.mul (Q16_16.div (Q16_16.ofInt 3) (Q16_16.ofInt 10)) avg)
|
||
mixed
|
||
|
||
/-- Phase 3: Single simulation tick.
|
||
Maps to visual step where surviving points cluster together.
|
||
|
||
One "frame" of physics simulation:
|
||
1. Blitter update (move toward perfect square)
|
||
2. Local gossip (pull toward neighbor confidence)
|
||
|
||
In WebGPU: Dispatch compute shader with barrier between steps. -/
|
||
def phase3Tick (dataset : Array TensorData) : Array TensorData :=
|
||
dataset.map (λ pt =>
|
||
if pt.isActive then
|
||
-- Step 1: Discrete Picard Integral (particle moves toward resonance)
|
||
let epsilon := Q16_16.div (Q16_16.ofInt 1) (Q16_16.ofInt 10) -- ε = 0.1
|
||
let (nextA, nextB) := picardBlitStep pt.a pt.b epsilon
|
||
|
||
-- Step 2: Local Gossip (pull toward neighbor confidence)
|
||
-- Neighbors are mod 8 in workgroup for L1 cache efficiency
|
||
let neighborIds := List.range 8 |>.map (λ i => (pt.id + i) % dataset.size)
|
||
let neighbors := neighborIds.filterMap (λ i =>
|
||
if i < dataset.size then some (dataset[i]!.confidence) else none)
|
||
let gossipConf := localGossip neighbors.toArray pt.confidence
|
||
|
||
{ pt with a := nextA, b := nextB, confidence := gossipConf }
|
||
else pt)
|
||
|
||
|
||
-- ════════════════════════════════════════════════════════════
|
||
-- §5 Full Pipeline Execution
|
||
-- ════════════════════════════════════════════════════════════
|
||
|
||
/-- Execute complete PIST simulation pipeline.
|
||
|
||
Steps:
|
||
1. Inject raw (a,b,confidence) states
|
||
2. Apply predictive pruning (kill doomed paths)
|
||
3. Run Blitter+Gossip for N frames
|
||
|
||
Returns final clustered states (the Perfect Square solutions).
|
||
|
||
Maps to WebGPU sequence:
|
||
• vkCmdDispatch(Phase1_Init)
|
||
• vkCmdDispatch(Phase2_Prune)
|
||
• for i in 0..frames: vkCmdDispatch(Phase3_BlitGossip) -/
|
||
def executePipeline (rawInputs : Array (Q16_16 × Q16_16 × Q16_16)) (frames : Nat) : Array TensorData :=
|
||
-- Step 1: Populate canvas
|
||
let injected := injectDataSlice rawInputs
|
||
|
||
-- Step 2: Apply heuristic (kill doomed paths instantly)
|
||
let pruned := phase2Pruning injected
|
||
|
||
-- Step 3: Run physics for 'frames' iterations
|
||
let rec loop (data : Array TensorData) (f : Nat) : Array TensorData :=
|
||
match f with
|
||
| 0 => data
|
||
| f' + 1 => loop (phase3Tick data) f'
|
||
loop pruned frames
|
||
|
||
/-- Execute pipeline from shell event indices.
|
||
Convenience wrapper for AVMR/SSMS integration. -/
|
||
def executeFromShellIndices (indices : List Nat) (frames : Nat) : Array TensorData :=
|
||
let rawInputs := indices.map (λ n =>
|
||
let s := shellState n
|
||
(Q16_16.ofInt (Int.ofNat s.a),
|
||
Q16_16.ofInt (Int.ofNat s.b),
|
||
Q16_16.ofInt (Int.ofNat (s.a * s.b))))
|
||
executePipeline rawInputs.toArray frames
|
||
|
||
|
||
-- ════════════════════════════════════════════════════════════
|
||
-- §4b Spectral Refinement — Gauss-Jordan Least-Squares (Rietveld-PIST)
|
||
-- ════════════════════════════════════════════════════════════
|
||
|
||
-- All types below are inlined from the former Astrophysics modules
|
||
-- to keep PistSimulation.lean self-contained and policy-clean.
|
||
|
||
/-- Bragg peak descriptor (pure Q16_16, no strings, no Float). -/
|
||
structure BraggPeak where
|
||
position : Q16_16
|
||
height : Q16_16
|
||
width : Q16_16
|
||
deriving Repr
|
||
|
||
/-- Simplified pseudo-Voigt profile in Q16_16.
|
||
Gaussian approximation: height / (1 + dist²/width²). -/
|
||
def pseudoVoigtQ16 (peak : BraggPeak) (x : Q16_16) : Q16_16 :=
|
||
let dist := if Q16_16.lt x peak.position then Q16_16.sub peak.position x else Q16_16.sub x peak.position
|
||
let distSq := Q16_16.mul dist dist
|
||
let widthSq := Q16_16.mul peak.width peak.width
|
||
let denom := Q16_16.add Q16_16.one (Q16_16.div distSq widthSq)
|
||
Q16_16.div peak.height denom
|
||
|
||
/-- χ² between observed and model value lists. -/
|
||
def chiSqWindow (observed model : List Q16_16) : Q16_16 :=
|
||
let pairs := observed.zip model
|
||
let sqDiffs := pairs.map (λ (o, m) =>
|
||
let diff := Q16_16.sub o m
|
||
Q16_16.mul diff diff)
|
||
sqDiffs.foldl Q16_16.add Q16_16.zero
|
||
|
||
/-- Magnetic domain descriptor (pure Q16_16, no string fields). -/
|
||
structure MagneticDomain where
|
||
domainId : Nat
|
||
centerFreq : Q16_16
|
||
sizeQ16 : Q16_16
|
||
totalIntensity : Q16_16
|
||
numPeaks : Nat
|
||
deltaNegative : Bool -- true = NÉEL (Δ<0), false = BLOCH (Δ≥0)
|
||
deriving Repr, Inhabited
|
||
|
||
/-- Scale width to domain size proxy. -/
|
||
def domainSizeFromWidth (w : Q16_16) : Q16_16 := Q16_16.mul w (Q16_16.ofNat 2)
|
||
|
||
/-- Semantic regime as an inductive type (compute-path, no strings). -/
|
||
inductive MagneticRegime
|
||
| bloch
|
||
| neel
|
||
| uglyAsymmetricPruning
|
||
| horribleManifoldTearing
|
||
deriving Repr, BEq
|
||
|
||
/-- Display function: maps inductive regime to human string at the boundary. -/
|
||
def regimeToString (r : MagneticRegime) : String :=
|
||
match r with
|
||
| .bloch => "BLOCH"
|
||
| .neel => "NÉEL"
|
||
| .uglyAsymmetricPruning => "uglyAsymmetricPruning"
|
||
| .horribleManifoldTearing => "horribleManifoldTearing"
|
||
|
||
/-- Classify a list of domains into a regime. -/
|
||
def domainRegime (domains : List MagneticDomain) : MagneticRegime :=
|
||
if domains.isEmpty then MagneticRegime.uglyAsymmetricPruning
|
||
else
|
||
let d := domains.head!
|
||
if d.numPeaks == 0 then MagneticRegime.uglyAsymmetricPruning
|
||
else if d.numPeaks > 1 then MagneticRegime.horribleManifoldTearing
|
||
else if d.deltaNegative then MagneticRegime.neel
|
||
else MagneticRegime.bloch
|
||
|
||
/-- Small-matrix row: up to 4 elements for 3×3 + RHS systems. -/
|
||
abbrev MatRow := Array Q16_16
|
||
|
||
/-- Augmented matrix for Gauss-Jordan: n rows × (n+1) columns. -/
|
||
abbrev AugMat := Array MatRow
|
||
|
||
def swapRows (m : AugMat) (i j : Nat) : AugMat :=
|
||
if i < m.size && j < m.size then
|
||
let ri := m[i]!
|
||
let rj := m[j]!
|
||
m.set! i rj |>.set! j ri
|
||
else m
|
||
|
||
def scaleRow (row : MatRow) (factor : Q16_16) : MatRow :=
|
||
row.map (λ v => Q16_16.mul v factor)
|
||
|
||
def addScaledRow (target : MatRow) (source : MatRow) (scale : Q16_16) : MatRow :=
|
||
target.zip source |>.map (λ (t, s) => Q16_16.add t (Q16_16.mul scale s))
|
||
|
||
private def findPivot (m : AugMat) (col startRow : Nat) : Nat :=
|
||
let rec search (r : Nat) (bestRow : Nat) (bestVal : Q16_16) : Nat :=
|
||
if r >= m.size then bestRow
|
||
else
|
||
let row := m[r]!
|
||
let v := row[col]!
|
||
let absV := Q16_16.abs v
|
||
if Q16_16.gt absV bestVal then search (r + 1) r absV
|
||
else search (r + 1) bestRow bestVal
|
||
search startRow startRow Q16_16.zero
|
||
|
||
private def eliminateColumn (m : AugMat) (col : Nat) (nRows : Nat) : AugMat :=
|
||
let pivotRow := m[col]!
|
||
let rec elim (m' : AugMat) (r : Nat) : AugMat :=
|
||
if r >= nRows then m'
|
||
else if r = col then elim m' (r + 1)
|
||
else
|
||
let row := m'[r]!
|
||
let factor := Q16_16.neg (row[col]!)
|
||
let newRow := addScaledRow row pivotRow factor
|
||
elim (m'.set! r newRow) (r + 1)
|
||
elim m 0
|
||
|
||
/-- Gauss-Jordan elimination for small linear systems (Ax = b).
|
||
Returns solution vector x (last column of reduced matrix).
|
||
Singular matrices yield a zero vector. -/
|
||
def gaussJordanSolve (mat : AugMat) : Array Q16_16 :=
|
||
let n := mat.size
|
||
let rec forward (m : AugMat) (col : Nat) : AugMat :=
|
||
match col with
|
||
| 0 => m
|
||
| c' + 1 =>
|
||
let pivotRow := findPivot m c' c'
|
||
let m1 := swapRows m c' pivotRow
|
||
let pivotRowData := m1[c']!
|
||
let pivotVal := pivotRowData[c']!
|
||
if pivotVal.val = 0 then m1
|
||
else
|
||
let invPivot := Q16_16.recip pivotVal
|
||
let rowNormed := scaleRow pivotRowData invPivot
|
||
let m2 := m1.set! c' rowNormed
|
||
let m3 := eliminateColumn m2 c' n
|
||
forward m3 c'
|
||
let reduced := forward mat n
|
||
reduced.map (λ row => if row.size > n then row[n]! else Q16_16.zero)
|
||
|
||
-- ════════════════════════════════════════════════════════════
|
||
-- §4c Matrix Packet (G = AᵀA, det, rank, spectral witness)
|
||
-- ════════════════════════════════════════════════════════════
|
||
|
||
def matGet (m : Array (Array Q16_16)) (i j : Nat) : Q16_16 :=
|
||
if h₁ : i < m.size then
|
||
let row := m[i]
|
||
if h₂ : j < row.size then row[j] else Q16_16.zero
|
||
else Q16_16.zero
|
||
|
||
def det3 (m : Array (Array Q16_16)) : Q16_16 :=
|
||
let a := matGet m 0 0; let b := matGet m 0 1; let c := matGet m 0 2
|
||
let d := matGet m 1 0; let e := matGet m 1 1; let f := matGet m 1 2
|
||
let g := matGet m 2 0; let h := matGet m 2 1; let i := matGet m 2 2
|
||
let term1 := Q16_16.mul a (Q16_16.sub (Q16_16.mul e i) (Q16_16.mul f h))
|
||
let term2 := Q16_16.mul b (Q16_16.sub (Q16_16.mul d i) (Q16_16.mul f g))
|
||
let term3 := Q16_16.mul c (Q16_16.sub (Q16_16.mul d h) (Q16_16.mul e g))
|
||
Q16_16.add term1 (Q16_16.sub term3 term2)
|
||
|
||
def trace3 (m : Array (Array Q16_16)) : Q16_16 :=
|
||
Q16_16.add (matGet m 0 0) (Q16_16.add (matGet m 1 1) (matGet m 2 2))
|
||
|
||
def rank3 (m : Array (Array Q16_16)) : Nat :=
|
||
if (det3 m).val != 0 then 3
|
||
else
|
||
let minors : List Q16_16 := [
|
||
Q16_16.sub (Q16_16.mul (matGet m 0 0) (matGet m 1 1)) (Q16_16.mul (matGet m 0 1) (matGet m 1 0)),
|
||
Q16_16.sub (Q16_16.mul (matGet m 0 0) (matGet m 1 2)) (Q16_16.mul (matGet m 0 2) (matGet m 1 0)),
|
||
Q16_16.sub (Q16_16.mul (matGet m 0 1) (matGet m 1 2)) (Q16_16.mul (matGet m 0 2) (matGet m 1 1)),
|
||
Q16_16.sub (Q16_16.mul (matGet m 0 0) (matGet m 2 1)) (Q16_16.mul (matGet m 0 1) (matGet m 2 0)),
|
||
Q16_16.sub (Q16_16.mul (matGet m 0 0) (matGet m 2 2)) (Q16_16.mul (matGet m 0 2) (matGet m 2 0)),
|
||
Q16_16.sub (Q16_16.mul (matGet m 0 1) (matGet m 2 2)) (Q16_16.mul (matGet m 0 2) (matGet m 2 1)),
|
||
Q16_16.sub (Q16_16.mul (matGet m 1 0) (matGet m 2 1)) (Q16_16.mul (matGet m 1 1) (matGet m 2 0)),
|
||
Q16_16.sub (Q16_16.mul (matGet m 1 0) (matGet m 2 2)) (Q16_16.mul (matGet m 1 2) (matGet m 2 0)),
|
||
Q16_16.sub (Q16_16.mul (matGet m 1 1) (matGet m 2 2)) (Q16_16.mul (matGet m 1 2) (matGet m 2 1))
|
||
]
|
||
if minors.any (λ v => v.val != 0) then 2
|
||
else
|
||
let elems := [matGet m 0 0, matGet m 0 1, matGet m 0 2,
|
||
matGet m 1 0, matGet m 1 1, matGet m 1 2,
|
||
matGet m 2 0, matGet m 2 1, matGet m 2 2]
|
||
if elems.any (λ v => v.val != 0) then 1
|
||
else 0
|
||
|
||
structure MatrixPacket where
|
||
gram : Array (Array Q16_16)
|
||
det : Q16_16
|
||
trace : Q16_16
|
||
rank : Nat
|
||
nullity : Nat
|
||
spectral : Array (Array Q16_16)
|
||
chiSq : Q16_16
|
||
deriving Repr
|
||
|
||
def spectralPlaceholder : Array (Array Q16_16) := #[
|
||
#[Q16_16.one, Q16_16.zero, Q16_16.zero],
|
||
#[Q16_16.zero, Q16_16.one, Q16_16.zero],
|
||
#[Q16_16.zero, Q16_16.zero, Q16_16.one]
|
||
]
|
||
|
||
def buildMatrixPacket (_window : List Q16_16) : MatrixPacket :=
|
||
let s0 := Q16_16.ofNat 8
|
||
let s1 := Q16_16.ofNat 28
|
||
let s2 := Q16_16.ofNat 140
|
||
let s3 := Q16_16.ofNat 784
|
||
let s4 := Q16_16.ofNat 4676
|
||
let gram := #[
|
||
#[s0, s1, s2],
|
||
#[s1, s2, s3],
|
||
#[s2, s3, s4]
|
||
]
|
||
let d := det3 gram
|
||
let r := rank3 gram
|
||
let t := trace3 gram
|
||
{ gram := gram, det := d, trace := t, rank := r, nullity := 3 - r,
|
||
spectral := spectralPlaceholder, chiSq := Q16_16.zero }
|
||
|
||
-- ════════════════════════════════════════════════════════════
|
||
-- §4d Quadratic LSQ Fit and Packet
|
||
-- ════════════════════════════════════════════════════════════
|
||
|
||
def quadraticFitCoeffs (window : List Q16_16) : Array Q16_16 :=
|
||
let n := 8
|
||
let s0 := Q16_16.ofNat n
|
||
let s1 := Q16_16.ofNat 28
|
||
let s2 := Q16_16.ofNat 140
|
||
let s3 := Q16_16.ofNat 784
|
||
let s4 := Q16_16.ofNat 4676
|
||
let rec computeSums (idx : Nat) (ys : List Q16_16)
|
||
(sy : Q16_16) (sxy : Q16_16) (sx2y : Q16_16) : Q16_16 × Q16_16 × Q16_16 :=
|
||
match ys with
|
||
| [] => (sy, sxy, sx2y)
|
||
| y :: rest =>
|
||
let x := Q16_16.ofNat idx
|
||
let x2 := Q16_16.mul x x
|
||
computeSums (idx + 1) rest
|
||
(Q16_16.add sy y)
|
||
(Q16_16.add sxy (Q16_16.mul x y))
|
||
(Q16_16.add sx2y (Q16_16.mul x2 y))
|
||
let (rhs0, rhs1, rhs2) := computeSums 0 window Q16_16.zero Q16_16.zero Q16_16.zero
|
||
let mat : AugMat := #[
|
||
#[s0, s1, s2, rhs0],
|
||
#[s1, s2, s3, rhs1],
|
||
#[s2, s3, s4, rhs2]
|
||
]
|
||
gaussJordanSolve mat
|
||
|
||
def quadraticDiscriminant (c0 c1 c2 : Q16_16) : Q16_16 :=
|
||
Q16_16.sub (Q16_16.mul c1 c1) (Q16_16.mul (Q16_16.mul (Q16_16.ofNat 4) c2) c0)
|
||
|
||
structure QuadraticPacket where
|
||
a : Q16_16
|
||
h : Q16_16
|
||
k : Q16_16
|
||
delta : Q16_16
|
||
width : Q16_16
|
||
chiSq : Q16_16
|
||
deriving Repr
|
||
|
||
def spectralEnergy (peak : BraggPeak) (window : List Q16_16) : Q16_16 :=
|
||
let modelVals := window.mapIdx (fun idx _obs => pseudoVoigtQ16 peak (Q16_16.ofNat idx))
|
||
chiSqWindow window modelVals
|
||
|
||
def coeffsToPacket (coeffs : Array Q16_16) (window : List Q16_16) : QuadraticPacket :=
|
||
if coeffs.size < 3 then
|
||
{ a := Q16_16.zero, h := Q16_16.zero, k := Q16_16.zero,
|
||
delta := Q16_16.zero, width := Q16_16.one, chiSq := Q16_16.zero }
|
||
else
|
||
let c0 := coeffs[0]!
|
||
let c1 := coeffs[1]!
|
||
let c2 := coeffs[2]!
|
||
let delta := quadraticDiscriminant c0 c1 c2
|
||
if c2.val = 0 then
|
||
let avg := Q16_16.div (window.foldl Q16_16.add Q16_16.zero) (Q16_16.ofNat window.length)
|
||
{ a := Q16_16.zero, h := Q16_16.zero, k := avg,
|
||
delta := delta, width := Q16_16.one, chiSq := Q16_16.zero }
|
||
else
|
||
let twoC2 := Q16_16.mul (Q16_16.ofNat 2) c2
|
||
let h := Q16_16.div (Q16_16.neg c1) twoC2
|
||
let c1sq := Q16_16.mul c1 c1
|
||
let fourC2 := Q16_16.mul (Q16_16.ofNat 4) c2
|
||
let k := Q16_16.sub c0 (Q16_16.div c1sq fourC2)
|
||
let negK := Q16_16.neg k
|
||
let twoA := Q16_16.mul (Q16_16.ofNat 2) c2
|
||
let ratio := if twoA.val = 0 then Q16_16.zero else Q16_16.div negK twoA
|
||
let widthApprox := if Q16_16.lt ratio Q16_16.zero then Q16_16.one
|
||
else Q16_16.add Q16_16.one (Q16_16.div ratio (Q16_16.ofNat 2))
|
||
let peak : BraggPeak := { position := h, height := k, width := widthApprox }
|
||
let chi := spectralEnergy peak window
|
||
{ a := c2, h := h, k := k, delta := delta,
|
||
width := widthApprox, chiSq := chi }
|
||
|
||
def packetToPeak (pkt : QuadraticPacket) : BraggPeak :=
|
||
{ position := pkt.h, height := pkt.k, width := pkt.width }
|
||
|
||
def coeffsToPeak (coeffs : Array Q16_16) : BraggPeak :=
|
||
packetToPeak (coeffsToPacket coeffs [])
|
||
|
||
def spectralRefineLsq (window : List Q16_16) : QuadraticPacket :=
|
||
coeffsToPacket (quadraticFitCoeffs window) window
|
||
|
||
def spectralPacketToDomain (pkt : QuadraticPacket) (domainId : Nat) : MagneticDomain :=
|
||
{ domainId := domainId
|
||
, centerFreq := pkt.h
|
||
, sizeQ16 := domainSizeFromWidth pkt.width
|
||
, totalIntensity := pkt.k
|
||
, numPeaks := 1
|
||
, deltaNegative := Q16_16.lt pkt.delta Q16_16.zero
|
||
}
|
||
|
||
def spectralWindowToRegime (window : List Q16_16) : MagneticRegime :=
|
||
let pkt := spectralRefineLsq window
|
||
let domain := spectralPacketToDomain pkt 0
|
||
domainRegime [domain]
|
||
|
||
def packetToRegime (pkt : QuadraticPacket) : MagneticRegime :=
|
||
let domain := spectralPacketToDomain pkt 0
|
||
domainRegime [domain]
|
||
|
||
-- ════════════════════════════════════════════════════════════
|
||
-- §4e Chaos Game Affine IFS Iterative Refinement
|
||
-- ════════════════════════════════════════════════════════════
|
||
|
||
structure ChaosState where
|
||
position : Q16_16
|
||
height : Q16_16
|
||
width : Q16_16
|
||
deriving Repr
|
||
|
||
def chaosContractionStep (s : ChaosState) (anchor : ChaosState) (r : Q16_16) : ChaosState :=
|
||
let oneMinusR := Q16_16.sub Q16_16.one r
|
||
let newPos := Q16_16.add (Q16_16.mul oneMinusR anchor.position) (Q16_16.mul r s.position)
|
||
let newHeight := Q16_16.add (Q16_16.mul oneMinusR anchor.height) (Q16_16.mul r s.height)
|
||
let newWidth := Q16_16.add (Q16_16.mul oneMinusR anchor.width) (Q16_16.mul r s.width)
|
||
{ position := newPos, height := newHeight, width := newWidth }
|
||
|
||
def chaosStateDiff (s1 s2 : ChaosState) : Q16_16 :=
|
||
let dPos := Q16_16.abs (Q16_16.sub s1.position s2.position)
|
||
let dHeight := Q16_16.abs (Q16_16.sub s1.height s2.height)
|
||
let dWidth := Q16_16.abs (Q16_16.sub s1.width s2.width)
|
||
Q16_16.add dPos (Q16_16.add dHeight dWidth)
|
||
|
||
def packetToChaosState (pkt : QuadraticPacket) : ChaosState :=
|
||
{ position := pkt.h, height := pkt.k, width := pkt.width }
|
||
|
||
def chaosStateToPacket (s : ChaosState) (window : List Q16_16) : QuadraticPacket :=
|
||
let peak : BraggPeak := { position := s.position, height := s.height, width := s.width }
|
||
let chi := spectralEnergy peak window
|
||
let a := Q16_16.recip (Q16_16.mul s.width s.width)
|
||
let delta := Q16_16.neg (Q16_16.mul (Q16_16.mul (Q16_16.ofNat 4) a) s.height)
|
||
{ a := a, h := s.position, k := s.height, delta := delta,
|
||
width := s.width, chiSq := chi }
|
||
|
||
def chaosConverge
|
||
(initial : ChaosState)
|
||
(anchor : ChaosState)
|
||
(_window : List Q16_16)
|
||
(r : Q16_16)
|
||
(threshold : Q16_16)
|
||
(maxSteps : Nat)
|
||
: ChaosState × Nat :=
|
||
let rec loop (s : ChaosState) (steps : Nat) (iter : Nat) : ChaosState × Nat :=
|
||
match steps with
|
||
| 0 => (s, iter)
|
||
| steps' + 1 =>
|
||
let sNext := chaosContractionStep s anchor r
|
||
let diff := chaosStateDiff sNext s
|
||
if Q16_16.lt diff threshold then (sNext, iter + 1)
|
||
else loop sNext steps' (iter + 1)
|
||
loop initial maxSteps 0
|
||
|
||
def spectralRefineChaos
|
||
(window : List Q16_16)
|
||
(r : Q16_16)
|
||
(threshold : Q16_16)
|
||
(maxSteps : Nat)
|
||
: QuadraticPacket × Nat :=
|
||
let lsqPkt := spectralRefineLsq window
|
||
let anchor := packetToChaosState lsqPkt
|
||
let perturbScale := Q16_16.ofRatio 1 10
|
||
let initial := {
|
||
position := Q16_16.add anchor.position (Q16_16.mul perturbScale Q16_16.one),
|
||
height := Q16_16.add anchor.height (Q16_16.mul perturbScale Q16_16.one),
|
||
width := Q16_16.add anchor.width (Q16_16.mul perturbScale Q16_16.one)
|
||
}
|
||
let (refined, steps) := chaosConverge initial anchor window r threshold maxSteps
|
||
let pkt := chaosStateToPacket refined window
|
||
(pkt, steps)
|
||
|
||
def spectralWindowToRegimeChaos (window : List Q16_16) : MagneticRegime :=
|
||
let (pkt, _steps) := spectralRefineChaos window (Q16_16.ofRatio 1 2) (Q16_16.ofRatio 1 100) 20
|
||
let domain := spectralPacketToDomain pkt 0
|
||
domainRegime [domain]
|
||
|
||
-- ════════════════════════════════════════════════════════════
|
||
-- §6 Verification Examples
|
||
-- ════════════════════════════════════════════════════════════
|
||
|
||
-- expect: true
|
||
#eval predictViability (Q16_16.ofInt 4) (Q16_16.ofInt 5) (Q16_16.ofInt 10)
|
||
-- expect: false
|
||
#eval predictViability (Q16_16.ofInt 1) (Q16_16.ofInt 20) (Q16_16.ofInt 10)
|
||
-- expect: (270531, 322240)
|
||
#eval picardBlitStep (Q16_16.ofInt 4) (Q16_16.ofInt 5) (Q16_16.div (Q16_16.ofInt 1) (Q16_16.ofInt 10))
|
||
-- expect: #[{ id := 0, a := 303810, b := 300899, confidence := 1310620, isActive := true }]
|
||
#eval executePipeline #[(Q16_16.ofInt 4, Q16_16.ofInt 5, Q16_16.ofInt 20)] 5
|
||
|
||
-- ════════════════════════════════════════════════════════════
|
||
-- §6b Spectral Refinement Verification
|
||
-- ════════════════════════════════════════════════════════════
|
||
|
||
def fixtureSpectralWindow : List Q16_16 := [
|
||
Q16_16.ofRawInt 655360, Q16_16.ofRawInt 1310720, Q16_16.ofRawInt 6553600, Q16_16.ofRawInt 2621440,
|
||
Q16_16.ofRawInt 1310720, Q16_16.ofRawInt 655360, Q16_16.ofRawInt 327680, Q16_16.ofRawInt 327680
|
||
]
|
||
|
||
-- expect: #[1488172, 1012395, -189193]
|
||
#eval! quadraticFitCoeffs fixtureSpectralWindow
|
||
-- expect: { a := -189193, h := 175346, k := 2842534, delta := 32823952, width := 311698, chiSq := 403453585 }
|
||
#eval! spectralRefineLsq fixtureSpectralWindow
|
||
-- expect: { position := 175346, height := 2842534, width := 311698 }
|
||
#eval! packetToPeak (spectralRefineLsq fixtureSpectralWindow)
|
||
-- expect: 403453585
|
||
#eval! (spectralRefineLsq fixtureSpectralWindow).chiSq
|
||
-- expect: 32823952
|
||
#eval! (spectralRefineLsq fixtureSpectralWindow).delta
|
||
-- expect: Semantics.PistSimulation.MagneticRegime.bloch
|
||
#eval! spectralWindowToRegime fixtureSpectralWindow
|
||
-- expect: Semantics.PistSimulation.MagneticRegime.bloch
|
||
#eval! packetToRegime (spectralRefineLsq fixtureSpectralWindow)
|
||
-- expect: { position := 175346, height := 2842534, width := 311698 }
|
||
#eval! packetToChaosState (spectralRefineLsq fixtureSpectralWindow)
|
||
-- expect: { position := 175346, height := 2842534, width := 311698 }
|
||
#eval! chaosContractionStep
|
||
(packetToChaosState (spectralRefineLsq fixtureSpectralWindow))
|
||
(packetToChaosState (spectralRefineLsq fixtureSpectralWindow))
|
||
(Q16_16.ofRatio 1 2)
|
||
-- expect: ({ a := 2893, h := 175550, k := 2842738, delta := -501955, width := 311902, chiSq := 403774829 }, 5)
|
||
#eval! spectralRefineChaos fixtureSpectralWindow (Q16_16.ofRatio 1 2) (Q16_16.ofRatio 1 100) 20
|
||
-- expect: Semantics.PistSimulation.MagneticRegime.neel
|
||
#eval! spectralWindowToRegimeChaos fixtureSpectralWindow
|
||
-- expect: { gram := #[#[524288, 1835008, 9175040], ...], sol := ... }
|
||
#eval! buildMatrixPacket fixtureSpectralWindow
|
||
|
||
-- ════════════════════════════════════════════════════════════
|
||
-- §7 TreeDIAT — Tree-to-Shell Coordinate Transform
|
||
-- ════════════════════════════════════════════════════════════
|
||
--
|
||
-- TreeDIAT follows the same evolution path as DIAT in DynamicCanal.lean:
|
||
-- DIAT → LanePayload → Lane → NodeState → N-DAG → Dynamic Canal → Throat
|
||
--
|
||
-- TreeDIAT mirrors each layer for tree-structured data (search traces,
|
||
-- n-gram tries, FAMM Delta-DAGs) so they can participate in the same
|
||
-- spectral refinement and regime classification as integer-shell data.
|
||
|
||
-- ── 7a. Tree Node (canonical input) ─────────────────────────
|
||
|
||
inductive TreeNode
|
||
| leaf (label : Nat)
|
||
| node (label : Nat) (left right : TreeNode)
|
||
deriving Repr
|
||
|
||
def treeMetrics (t : TreeNode) : Nat × Nat × Nat × Nat :=
|
||
let rec go (t : TreeNode) (depth : Nat) : Nat × Nat × Nat × Nat :=
|
||
match t with
|
||
| TreeNode.leaf lbl =>
|
||
(depth + 1, 1, 1, lbl)
|
||
| TreeNode.node lbl l r =>
|
||
let (dL, leafL, nodeL, maxL) := go l (depth + 1)
|
||
let (dR, leafR, nodeR, maxR) := go r (depth + 1)
|
||
(max dL dR, leafL + leafR, nodeL + nodeR + 1, max (max lbl maxL) maxR)
|
||
go t 0
|
||
|
||
-- ── 7b. TreeDIAT (structural feature vector) ──────────────
|
||
|
||
structure TreeDIAT where
|
||
depth : Nat
|
||
leafCount : Nat
|
||
nodeCount : Nat
|
||
labelCount : Nat
|
||
embeddingScore : Q16_16
|
||
deriving Repr, Inhabited
|
||
|
||
def treeDIATEmbeddingScore (depth leafCount labelCount nodeCount : Nat) : Q16_16 :=
|
||
if nodeCount = 0 then Q16_16.zero
|
||
else
|
||
let num := Q16_16.ofNat leafCount
|
||
let den := Q16_16.ofNat (depth * labelCount + 1)
|
||
Q16_16.div num den
|
||
|
||
def treeToDIAT (t : TreeNode) : TreeDIAT :=
|
||
let (d, leafC, nodeC, maxLbl) := treeMetrics t
|
||
let lblC := maxLbl + 1
|
||
let score := treeDIATEmbeddingScore d leafC lblC nodeC
|
||
{ depth := d, leafCount := leafC, nodeCount := nodeC, labelCount := lblC, embeddingScore := score }
|
||
|
||
/-- Normalised embedding score = score / (1 + score), in [0,1]. -/
|
||
def treeDIATNormEmbedding (td : TreeDIAT) : Q16_16 :=
|
||
let s := td.embeddingScore
|
||
let one := Q16_16.one
|
||
Q16_16.div s (Q16_16.add one s)
|
||
|
||
-- ── 7c. TreeLanePayload (analogous to LanePayload) ─────────
|
||
|
||
structure TreeLanePayload where
|
||
diat : TreeDIAT
|
||
codonWindow : UInt32
|
||
metadata : Q16_16
|
||
deriving Repr, Inhabited
|
||
|
||
-- ── 7d. TreeLane (physics state, analogous to Lane) ────────
|
||
|
||
structure TreeLane where
|
||
active : Bool
|
||
node : UInt32
|
||
pos : Q16_16 × Q16_16 × Q16_16
|
||
phase : Q16_16
|
||
stress : Q16_16
|
||
pressure : Q16_16
|
||
lambdaEff : Q16_16
|
||
energy : Q16_16
|
||
mismatch : Q16_16
|
||
regime : MagneticRegime
|
||
payload : TreeLanePayload
|
||
deriving Repr
|
||
|
||
-- ── 7e. TreeNodeState (analogous to NodeState) ─────────────
|
||
|
||
structure TreeNodeState where
|
||
diatState : Q16_16
|
||
waveState : Q16_16
|
||
timeState : Q16_16
|
||
deriving Repr
|
||
|
||
-- ── 7f. TreeEdge / TreeN-DAG (graph topology) ────────────
|
||
|
||
structure TreeEdge where
|
||
src : UInt32
|
||
dst : UInt32
|
||
torsion : Q16_16 -- parent-child rotation measure
|
||
loss : Q16_16 -- embedding cost of this edge
|
||
deriving Repr
|
||
|
||
structure TreeNDAG where
|
||
nodes : Array TreeNodeState
|
||
edges : Array TreeEdge
|
||
deriving Repr
|
||
|
||
-- ── 7g. TreeDynamicCanal (pressure-adaptive transport) ─────
|
||
|
||
/-- Effective resistance for tree-structured flow.
|
||
λ_eff(P) = λ₀ / (1 + ξ · P · depth)
|
||
Deep trees with high pressure become bottlenecks. -/
|
||
def treeDynamicCanalLambda (lambda0 xi pressure : Q16_16) (depth : Nat) : Q16_16 :=
|
||
let depthQ := Q16_16.ofNat depth
|
||
let xiP := Q16_16.mul (Q16_16.mul xi pressure) depthQ
|
||
let denom := Q16_16.add Q16_16.one xiP
|
||
Q16_16.div lambda0 denom
|
||
|
||
/-- Tree canal compliance = 1 / λ_eff. -/
|
||
def treeCanalCompliance (lambda0 xi pressure : Q16_16) (depth : Nat) : Q16_16 :=
|
||
Q16_16.recip (treeDynamicCanalLambda lambda0 xi pressure depth)
|
||
|
||
-- ── 7h. TreeThroat (regime transition classifier) ──────────
|
||
|
||
inductive TreeThroatClass
|
||
| stableBridge -- bushy, low pressure, high embeddability
|
||
| lossyChannel -- moderate, some pressure loss
|
||
| rupture -- stringy, high pressure, low embeddability
|
||
deriving Repr, BEq
|
||
|
||
/-- Classify a TreeLane by its physics state. -/
|
||
def classifyTreeThroat (lane : TreeLane) : TreeThroatClass :=
|
||
let td := lane.payload.diat
|
||
let normEmbed := treeDIATNormEmbedding td
|
||
if Q16_16.gt normEmbed (Q16_16.ofRatio 3 4) && Q16_16.lt lane.pressure (Q16_16.ofRatio 1 2) then
|
||
TreeThroatClass.stableBridge
|
||
else if Q16_16.lt normEmbed (Q16_16.ofRatio 1 4) || Q16_16.gt lane.pressure (Q16_16.ofRatio 3 2) then
|
||
TreeThroatClass.rupture
|
||
else
|
||
TreeThroatClass.lossyChannel
|
||
|
||
-- ── 7i. TreeSequenceRegime (meta-classifier) ──────────────
|
||
|
||
def treeSequenceRegime (seq : List TreeDIAT) : MagneticRegime :=
|
||
if seq.isEmpty then MagneticRegime.uglyAsymmetricPruning
|
||
else
|
||
let len := seq.length
|
||
if len > 1000000 then MagneticRegime.horribleManifoldTearing
|
||
else
|
||
let avgScore := Q16_16.div (seq.foldl (λ acc td => Q16_16.add acc td.embeddingScore) Q16_16.zero) (Q16_16.ofNat len)
|
||
if Q16_16.lt avgScore (Q16_16.ofRatio 1 10) then MagneticRegime.uglyAsymmetricPruning
|
||
else MagneticRegime.bloch
|
||
|
||
-- ── 7j. Projection into chaos-game space ────────────────────
|
||
|
||
def treeDIATToChaosState (td : TreeDIAT) : ChaosState :=
|
||
{ position := td.embeddingScore
|
||
, height := Q16_16.ofNat td.depth
|
||
, width := Q16_16.ofNat td.nodeCount }
|
||
|
||
-- ════════════════════════════════════════════════════════════
|
||
-- §7k Verification Fixtures
|
||
-- ════════════════════════════════════════════════════════════
|
||
|
||
def fixtureBushyTree : TreeNode :=
|
||
TreeNode.node 0
|
||
(TreeNode.node 1 (TreeNode.leaf 0) (TreeNode.leaf 1))
|
||
(TreeNode.node 0 (TreeNode.leaf 1) (TreeNode.leaf 0))
|
||
|
||
def fixtureStringyTree : TreeNode :=
|
||
TreeNode.node 0
|
||
(TreeNode.node 1
|
||
(TreeNode.node 2
|
||
(TreeNode.node 3 (TreeNode.leaf 0) (TreeNode.leaf 0))
|
||
(TreeNode.leaf 0))
|
||
(TreeNode.leaf 0))
|
||
(TreeNode.leaf 0)
|
||
|
||
def fixtureBalancedTree : TreeNode :=
|
||
TreeNode.node 2
|
||
(TreeNode.node 1 (TreeNode.leaf 0) (TreeNode.leaf 2))
|
||
(TreeNode.node 0 (TreeNode.leaf 1) (TreeNode.leaf 2))
|
||
|
||
/- Tree metrics and DIAT encoding. -/
|
||
-- expect: (3, 4, 7, 1)
|
||
#eval! treeMetrics fixtureBushyTree
|
||
-- expect: { depth := 3, leafCount := 4, nodeCount := 7, labelCount := 2, embeddingScore := 37449 }
|
||
#eval! treeToDIAT fixtureBushyTree
|
||
-- expect: (5, 5, 9, 3)
|
||
#eval! treeMetrics fixtureStringyTree
|
||
-- expect: { depth := 5, leafCount := 5, nodeCount := 9, labelCount := 4, embeddingScore := 15603 }
|
||
#eval! treeToDIAT fixtureStringyTree
|
||
-- expect: (3, 4, 7, 2)
|
||
#eval! treeMetrics fixtureBalancedTree
|
||
-- expect: { depth := 3, leafCount := 4, nodeCount := 7, labelCount := 3, embeddingScore := 26214 }
|
||
#eval! treeToDIAT fixtureBalancedTree
|
||
|
||
/- Embedding scores: bushy > balanced > stringy. -/
|
||
-- expect: 37449
|
||
#eval! (treeToDIAT fixtureBushyTree).embeddingScore
|
||
-- expect: 15603
|
||
#eval! (treeToDIAT fixtureStringyTree).embeddingScore
|
||
-- expect: 26214
|
||
#eval! (treeToDIAT fixtureBalancedTree).embeddingScore
|
||
|
||
/- Normalised embedding scores. -/
|
||
-- expect: 23831
|
||
#eval! treeDIATNormEmbedding (treeToDIAT fixtureBushyTree)
|
||
-- expect: 12602
|
||
#eval! treeDIATNormEmbedding (treeToDIAT fixtureStringyTree)
|
||
-- expect: 18724
|
||
#eval! treeDIATNormEmbedding (treeToDIAT fixtureBalancedTree)
|
||
|
||
/- Chaos-state projection. -/
|
||
-- expect: { position := 37449, height := 196608, width := 458752 }
|
||
#eval! treeDIATToChaosState (treeToDIAT fixtureBushyTree)
|
||
|
||
/- Regime classification. -/
|
||
-- expect: Semantics.PistSimulation.MagneticRegime.bloch
|
||
#eval! treeSequenceRegime [treeToDIAT fixtureBushyTree]
|
||
-- expect: Semantics.PistSimulation.MagneticRegime.bloch
|
||
#eval! treeSequenceRegime [treeToDIAT fixtureStringyTree]
|
||
-- expect: Semantics.PistSimulation.MagneticRegime.bloch
|
||
#eval! treeSequenceRegime [treeToDIAT fixtureBushyTree, treeToDIAT fixtureBalancedTree, treeToDIAT fixtureStringyTree]
|
||
|
||
/- Chaos-game contraction on tree DIAT anchor. -/
|
||
-- expect: ({ position := 37652, height := 196812, width := 458956 }, 5)
|
||
#eval! let td := treeToDIAT fixtureBushyTree;
|
||
let anchor := treeDIATToChaosState td;
|
||
let perturb := { position := Q16_16.add anchor.position (Q16_16.ofRatio 1 10)
|
||
, height := Q16_16.add anchor.height (Q16_16.ofRatio 1 10)
|
||
, width := Q16_16.add anchor.width (Q16_16.ofRatio 1 10) };
|
||
chaosConverge perturb anchor [] (Q16_16.ofRatio 1 2) (Q16_16.ofRatio 1 100) 20
|
||
|
||
/- TreeLane construction and throat classification. -/
|
||
-- expect: Semantics.PistSimulation.TreeThroatClass.lossyChannel
|
||
#eval! let td := treeToDIAT fixtureBushyTree;
|
||
let lane : TreeLane := {
|
||
active := true, node := 0,
|
||
pos := (Q16_16.ofNat td.depth, Q16_16.ofNat td.leafCount, Q16_16.ofNat td.nodeCount),
|
||
phase := td.embeddingScore, stress := Q16_16.ofRatio 1 10,
|
||
pressure := Q16_16.ofRatio 1 4, lambdaEff := Q16_16.one,
|
||
energy := Q16_16.ofNat td.nodeCount, mismatch := Q16_16.zero,
|
||
regime := MagneticRegime.bloch,
|
||
payload := { diat := td, codonWindow := 0, metadata := Q16_16.zero }
|
||
};
|
||
classifyTreeThroat lane
|
||
|
||
/- Stringy tree lane → rupture throat. -/
|
||
-- expect: Semantics.PistSimulation.TreeThroatClass.rupture
|
||
#eval! let td := treeToDIAT fixtureStringyTree;
|
||
let lane : TreeLane := {
|
||
active := true, node := 1,
|
||
pos := (Q16_16.ofNat td.depth, Q16_16.ofNat td.leafCount, Q16_16.ofNat td.nodeCount),
|
||
phase := td.embeddingScore, stress := Q16_16.ofRatio 3 10,
|
||
pressure := Q16_16.ofRatio 2 1, lambdaEff := Q16_16.ofRatio 1 2,
|
||
energy := Q16_16.ofNat td.nodeCount, mismatch := Q16_16.ofRatio 1 5,
|
||
regime := MagneticRegime.uglyAsymmetricPruning,
|
||
payload := { diat := td, codonWindow := 0, metadata := Q16_16.zero }
|
||
};
|
||
classifyTreeThroat lane
|
||
|
||
/- TreeDynamicCanal: λ_eff for bushy vs stringy at same pressure. -/
|
||
-- expect: 56989
|
||
#eval! treeDynamicCanalLambda Q16_16.one (Q16_16.ofRatio 1 10) (Q16_16.ofRatio 1 2)
|
||
(treeToDIAT fixtureBushyTree).depth
|
||
-- expect: 52431
|
||
#eval! treeDynamicCanalLambda Q16_16.one (Q16_16.ofRatio 1 10) (Q16_16.ofRatio 1 2)
|
||
(treeToDIAT fixtureStringyTree).depth
|
||
|
||
/- TreeN-DAG witness: 2-node, 1-edge graph. -/
|
||
-- expect: { nodes := #[...], edges := #[...] }
|
||
#eval! let n1 : TreeNodeState := { diatState := (treeToDIAT fixtureBushyTree).embeddingScore, waveState := Q16_16.zero, timeState := Q16_16.zero };
|
||
let n2 : TreeNodeState := { diatState := (treeToDIAT fixtureStringyTree).embeddingScore, waveState := Q16_16.zero, timeState := Q16_16.one };
|
||
let e1 : TreeEdge := { src := 0, dst := 1, torsion := Q16_16.ofRatio 1 4, loss := Q16_16.ofRatio 1 10 };
|
||
TreeNDAG.mk #[n1, n2] #[e1]
|
||
|
||
-- ════════════════════════════════════════════════════════════
|
||
-- §8 PhiNUVMAP — Golden-Ratio Fractal 16D Coordinate System
|
||
-- ════════════════════════════════════════════════════════════
|
||
--
|
||
-- PhiNUVMAP lifts NUVMAP into a 16D golden-ratio-scaled fractal space.
|
||
--
|
||
-- Core insight: φ = (1+√5)/2 is the unique number where φ^2 = φ + 1.
|
||
-- This gives the Fibonacci recurrence, which yields self-similar tilings
|
||
-- at all scales — the definition of a fractal.
|
||
--
|
||
-- In PhiNUVMAP:
|
||
-- • Coordinates scale by φ (zoom in) or φ^(-1) (zoom out / contract)
|
||
-- • The space is naturally non-uniform: denser near the center
|
||
-- • The golden contraction law s' = c + φ^(-1)·(s-c) is exact
|
||
-- • TreeDIAT states project into this 16D space and contract toward
|
||
-- their anchor at the golden rate
|
||
|
||
-- ── 8a. Golden ratio in Q16_16 ──────────────────────────────
|
||
|
||
/-- φ ≈ 4181/2584 = 1.6180339887… (error < 10⁻⁹).
|
||
Both 4181 and 2584 are Fibonacci numbers, so this is the canonical
|
||
rational approximation for fixed-point golden ratio work. -/
|
||
def phiQ16_16 : Q16_16 := Q16_16.ofRatio 4181 2584
|
||
|
||
/-- φ⁻¹ ≈ 2584/4181 = 0.6180339887…
|
||
Satisfies φ · φ⁻¹ = 1 in real arithmetic; in Q16_16 the product is
|
||
within 1 ULP of one. -/
|
||
def phiInvQ16_16 : Q16_16 := Q16_16.ofRatio 2584 4181
|
||
|
||
/-- φ² = φ + 1 (the defining identity) approximated in Q16_16.
|
||
Used for fractal self-similarity checks. -/
|
||
def phiSqQ16_16 : Q16_16 := Q16_16.add phiQ16_16 Q16_16.one
|
||
|
||
-- ── 8b. 16D vector utilities ───────────────────────────────
|
||
|
||
/-- Component-wise subtraction of two 16D vectors. -/
|
||
def vec16Sub (a b : Array Q16_16) : Array Q16_16 :=
|
||
a.zip b |>.map (λ (x, y) => Q16_16.sub x y)
|
||
|
||
/-- Component-wise addition of two 16D vectors. -/
|
||
def vec16Add (a b : Array Q16_16) : Array Q16_16 :=
|
||
a.zip b |>.map (λ (x, y) => Q16_16.add x y)
|
||
|
||
/-- Component-wise scalar multiplication of a 16D vector. -/
|
||
def vec16Scale (s : Q16_16) (v : Array Q16_16) : Array Q16_16 :=
|
||
v.map (λ x => Q16_16.mul s x)
|
||
|
||
/-- 16D zero vector. -/
|
||
def vec16Zero : Array Q16_16 :=
|
||
#[Q16_16.zero, Q16_16.zero, Q16_16.zero, Q16_16.zero,
|
||
Q16_16.zero, Q16_16.zero, Q16_16.zero, Q16_16.zero,
|
||
Q16_16.zero, Q16_16.zero, Q16_16.zero, Q16_16.zero,
|
||
Q16_16.zero, Q16_16.zero, Q16_16.zero, Q16_16.zero]
|
||
|
||
-- ── 8c. PhiNUVMAP spectral mode (local, mirrors NUVMAP) ───
|
||
|
||
inductive PhiSpectralMode
|
||
| dc
|
||
| lowFreq
|
||
| midFreq
|
||
| highFreq
|
||
| ultraFreq
|
||
| transient
|
||
deriving Repr, BEq
|
||
|
||
-- ── 8d. PhiNUVMAP structure ───────────────────────────────
|
||
|
||
/-- PhiNUVMAP: a 16D golden-ratio fractal coordinate system.
|
||
Fields:
|
||
center — shared 16D attractor point (the "golden center")
|
||
coords — list of 16D coordinates (tree states, anchors, etc.)
|
||
scaleLevel— fractal zoom level k (coordinates conceptually scaled by φ^k)
|
||
spectralMode — dc / low / mid / high / ultra / transient -/
|
||
structure PhiNUVMAP where
|
||
center : Array Q16_16 -- length 16
|
||
coords : Array (Array Q16_16) -- each length 16
|
||
scaleLevel : Nat -- zoom level k
|
||
spectralMode : PhiSpectralMode
|
||
deriving Repr
|
||
|
||
-- ── 8d. Golden contraction law ────────────────────────────
|
||
|
||
/-- Golden contraction: s' = center + φ⁻¹ · (s - center).
|
||
After t iterations: ||s(t) - c|| = φ⁻ᵗ · ||s(0) - c||.
|
||
This is the fractal self-similarity engine. -/
|
||
def phiContract (state center : Array Q16_16) : Array Q16_16 :=
|
||
let diff := vec16Sub state center
|
||
let scaled := vec16Scale phiInvQ16_16 diff
|
||
vec16Add center scaled
|
||
|
||
/-- Multi-step golden contraction.
|
||
Returns (final_state, number_of_steps). -/
|
||
def phiContractN (state center : Array Q16_16) (steps : Nat) : Array Q16_16 :=
|
||
let rec loop (s : Array Q16_16) (n : Nat) : Array Q16_16 :=
|
||
match n with
|
||
| 0 => s
|
||
| n' + 1 => loop (phiContract s center) n'
|
||
loop state steps
|
||
|
||
-- ── 8e. Fractal zoom operations ────────────────────────────
|
||
|
||
/-- Zoom IN by one fractal level: multiply coordinates by φ.
|
||
Conceptually: coord' = φ · coord. -/
|
||
def phiZoomIn (coord : Array Q16_16) : Array Q16_16 :=
|
||
vec16Scale phiQ16_16 coord
|
||
|
||
/-- Zoom OUT by one fractal level: multiply coordinates by φ⁻¹.
|
||
This is the same as one golden contraction step toward origin. -/
|
||
def phiZoomOut (coord : Array Q16_16) : Array Q16_16 :=
|
||
vec16Scale phiInvQ16_16 coord
|
||
|
||
/-- Scale a PhiNUVMAP coordinate by φ^k for arbitrary integer k.
|
||
Positive k = zoom in (enlarge); negative k = zoom out (shrink). -/
|
||
def phiScaleBy (coord : Array Q16_16) (k : Int) : Array Q16_16 :=
|
||
if k >= 0 then
|
||
let rec zoomIn (c : Array Q16_16) (n : Nat) : Array Q16_16 :=
|
||
match n with
|
||
| 0 => c
|
||
| n' + 1 => zoomIn (phiZoomIn c) n'
|
||
zoomIn coord k.toNat
|
||
else
|
||
let rec zoomOut (c : Array Q16_16) (n : Nat) : Array Q16_16 :=
|
||
match n with
|
||
| 0 => c
|
||
| n' + 1 => zoomOut (phiZoomOut c) n'
|
||
zoomOut coord (-k).toNat
|
||
|
||
-- ── 8f. Tree-to-PhiNUVMAP projection ──────────────────────
|
||
|
||
/-- Project a TreeDIAT into the 16D φ-NUVMAP space.
|
||
Maps tree metrics into dimensions 0-5, derived features into 6-11,
|
||
and pads with zeros for 12-15. The 16th position is filled with
|
||
the embedding score as the "attractor weight". -/
|
||
def treeDIATToPhiNUVMAP (td : TreeDIAT) : Array Q16_16 :=
|
||
let d := Q16_16.ofNat td.depth
|
||
let lc := Q16_16.ofNat td.leafCount
|
||
let nc := Q16_16.ofNat td.nodeCount
|
||
let lbl := Q16_16.ofNat td.labelCount
|
||
let score := td.embeddingScore
|
||
let normScore := treeDIATNormEmbedding td
|
||
let d_times_lbl := Q16_16.mul d lbl
|
||
let lc_over_nc := if td.nodeCount = 0 then Q16_16.zero else Q16_16.div lc nc
|
||
let score_times_d := Q16_16.mul score d
|
||
-- Dimensions 0-7: primary tree features
|
||
-- Dimensions 8-15: structural pressure, normalized features, padding
|
||
#[score, d, lc, nc, lbl, d_times_lbl, lc_over_nc, score_times_d,
|
||
normScore, Q16_16.sub Q16_16.one normScore, -- embeddability + residual
|
||
Q16_16.div d (Q16_16.ofNat 10), -- depth/10 (scale proxy)
|
||
Q16_16.div nc (Q16_16.ofNat 100), -- nodeCount/100 (mass proxy)
|
||
Q16_16.zero, Q16_16.zero, Q16_16.zero, Q16_16.zero]
|
||
|
||
/-- Build a PhiNUVMAP from a TreeDIAT with a given center and scale level.
|
||
The tree state becomes the single coordinate; the center is provided
|
||
by the caller (typically an anchor tree or the global golden center). -/
|
||
def treeDIATToPhiNUVMAPState (td : TreeDIAT) (center : Array Q16_16)
|
||
(scaleLevel : Nat) (mode : PhiSpectralMode) : PhiNUVMAP :=
|
||
{ center := center
|
||
, coords := #[treeDIATToPhiNUVMAP td]
|
||
, scaleLevel := scaleLevel
|
||
, spectralMode := mode }
|
||
|
||
-- ── 8g. 16D chaos game with φ-contraction ─────────────────
|
||
|
||
/-- One step of the 16D φ-NUVMAP chaos game.
|
||
X_{t+1} = anchor + φ⁻¹ · (X_t - anchor) + ε
|
||
where ε is a small perturbation (simulates exploration).
|
||
In this formal version, ε is deterministic (tests stability). -/
|
||
def phiNUVMAPChaosStep (state anchor : Array Q16_16) (epsilon : Array Q16_16)
|
||
: Array Q16_16 :=
|
||
let contracted := phiContract state anchor
|
||
vec16Add contracted epsilon
|
||
|
||
/-- Run the 16D φ-NUVMAP chaos game for N steps.
|
||
Returns the final state. -/
|
||
def phiNUVMAPChaosRun (initial anchor : Array Q16_16) (epsilon : Array Q16_16)
|
||
(steps : Nat) : Array Q16_16 :=
|
||
let rec loop (s : Array Q16_16) (n : Nat) : Array Q16_16 :=
|
||
match n with
|
||
| 0 => s
|
||
| n' + 1 => loop (phiNUVMAPChaosStep s anchor epsilon) n'
|
||
loop initial steps
|
||
|
||
-- ── 8h. Verification witnesses ────────────────────────────
|
||
|
||
/- Golden ratio approximation witness: φ · φ⁻¹ ≈ 1. -/
|
||
-- expect: 65534
|
||
#eval! Q16_16.mul phiQ16_16 phiInvQ16_16
|
||
|
||
/- φ² = φ + 1 witness. -/
|
||
-- expect: 171573
|
||
#eval! Q16_16.mul phiQ16_16 phiQ16_16
|
||
-- expect: 171575
|
||
#eval! phiSqQ16_16
|
||
|
||
/- Golden contraction of a 16D state toward origin.
|
||
After one step, each component should be ≈ 0.618 × original. -/
|
||
-- expect: #[405030, 810060, 1215090, 1620120, 2025150, 2430180, 2835210, 3240240, 3645270, 4050300, 4455330, 4860360, 5265390, ...]
|
||
#eval! let s := #[Q16_16.ofNat 10, Q16_16.ofNat 20, Q16_16.ofNat 30, Q16_16.ofNat 40,
|
||
Q16_16.ofNat 50, Q16_16.ofNat 60, Q16_16.ofNat 70, Q16_16.ofNat 80,
|
||
Q16_16.ofNat 90, Q16_16.ofNat 100, Q16_16.ofNat 110, Q16_16.ofNat 120,
|
||
Q16_16.ofNat 130, Q16_16.ofNat 140, Q16_16.ofNat 150, Q16_16.ofNat 160];
|
||
phiContract s vec16Zero
|
||
|
||
/- After 5 contraction steps toward origin: state ≈ φ⁻⁵ · initial.
|
||
φ⁻⁵ ≈ 0.090, so 160 → ≈ 14.4. -/
|
||
-- expect: #[59089, 118179, 177270, 236360, 295450, 354541, 413631, 472721, 531812, 590902, 649992, 709083, 768173, 827263, 886354, ...]
|
||
#eval! let s := #[Q16_16.ofNat 10, Q16_16.ofNat 20, Q16_16.ofNat 30, Q16_16.ofNat 40,
|
||
Q16_16.ofNat 50, Q16_16.ofNat 60, Q16_16.ofNat 70, Q16_16.ofNat 80,
|
||
Q16_16.ofNat 90, Q16_16.ofNat 100, Q16_16.ofNat 110, Q16_16.ofNat 120,
|
||
Q16_16.ofNat 130, Q16_16.ofNat 140, Q16_16.ofNat 150, Q16_16.ofNat 160];
|
||
phiContractN s vec16Zero 5
|
||
|
||
/- Fractal zoom: zoom in ×1 then out ×1 = identity (up to rounding). -/
|
||
-- expect: #[6553493, 13106987, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
|
||
#eval! let c := #[Q16_16.ofNat 100, Q16_16.ofNat 200, Q16_16.zero, Q16_16.zero,
|
||
Q16_16.zero, Q16_16.zero, Q16_16.zero, Q16_16.zero,
|
||
Q16_16.zero, Q16_16.zero, Q16_16.zero, Q16_16.zero,
|
||
Q16_16.zero, Q16_16.zero, Q16_16.zero, Q16_16.zero];
|
||
phiZoomOut (phiZoomIn c)
|
||
|
||
/- TreeDIAT projected into 16D φ-NUVMAP space. -/
|
||
-- expect: #[37449, 196608, 262144, 458752, 131072, 393216, 37449, 112347, 23831, 41705, 19660, 4587, 0, 0, 0, 0]
|
||
#eval! treeDIATToPhiNUVMAP (treeToDIAT fixtureBushyTree)
|
||
|
||
/- TreeDIAT projected into 16D φ-NUVMAP space (stringy). -/
|
||
-- expect: #[15603, 327680, 327680, 589824, 262144, 1310720, 36408, 78015, 12602, 52934, 32768, 5898, 0, 0, 0, 0]
|
||
#eval! treeDIATToPhiNUVMAP (treeToDIAT fixtureStringyTree)
|
||
|
||
/- Golden contraction of bushy-tree 16D state toward stringy-tree 16D state.
|
||
The bushy tree should contract toward the stringy-tree anchor. -/
|
||
-- expect: #[29104, 246674, 287177, 508818, 181138, 743678, 37051, 99233, 19541, 45994, 24666, 5087, 0, 0, 0, 0]
|
||
#eval! let bushy16 := treeDIATToPhiNUVMAP (treeToDIAT fixtureBushyTree);
|
||
let stringy16 := treeDIATToPhiNUVMAP (treeToDIAT fixtureStringyTree);
|
||
phiContract bushy16 stringy16
|
||
|
||
/- 16D φ-NUVMAP chaos game: bushy tree contracts toward stringy-tree anchor
|
||
with small deterministic perturbation, 10 steps. -/
|
||
-- expect: #[15779, 326613, 327146, 588757, 261077, 1303259, 36415, 78293, 12692, 52841, 32660, 5886, 0, 0, 0, 0]
|
||
#eval! let bushy16 := treeDIATToPhiNUVMAP (treeToDIAT fixtureBushyTree);
|
||
let stringy16 := treeDIATToPhiNUVMAP (treeToDIAT fixtureStringyTree);
|
||
let eps := vec16Scale (Q16_16.ofRatio 1 100) vec16Zero; -- zero perturbation for stability
|
||
phiNUVMAPChaosRun bushy16 stringy16 eps 10
|
||
|
||
/- Scale level witness: bushy tree at scale level 0. -/
|
||
-- expect: { center := #[0, 0, ...], coords := #[#[37449, 196608, ...]], scaleLevel := 0, mode := dc }
|
||
#eval! treeDIATToPhiNUVMAPState (treeToDIAT fixtureBushyTree) vec16Zero 0 PhiSpectralMode.dc
|
||
|
||
/- Scale level witness: stringy tree at scale level 3 (zoomed in). -/
|
||
-- expect: { center := #[0, 0, ...], coords := #[#[15603, 327680, ...]], scaleLevel := 3, mode := transient }
|
||
#eval! treeDIATToPhiNUVMAPState (treeToDIAT fixtureStringyTree) vec16Zero 3 PhiSpectralMode.transient
|
||
|
||
-- ════════════════════════════════════════════════════════════
|
||
-- §9 Burgers-PhiNUVMAP Bridge
|
||
-- ════════════════════════════════════════════════════════════
|
||
--
|
||
-- The Burgers equation u_t + u·u_x = ν·u_xx exhibits two regimes:
|
||
-- • ν large → smooth, diffusive, oscillatory (NÉEL / Δ < 0)
|
||
-- • ν → 0 → shock formation, discontinuous (BLOCH / Δ ≥ 0)
|
||
--
|
||
-- The spectral pipeline in PistSimulation.lean (8-bin window, quadratic
|
||
-- fit, discriminant gate) classifies Burgers solutions directly.
|
||
--
|
||
-- PhiNUVMAP adds a 16D golden-ratio fractal parameter space where:
|
||
-- • Each BurgersState maps to a 16D vector
|
||
-- • Golden contraction s' = c + φ⁻¹·(s-c) models viscous dissipation
|
||
-- • Shock detection = regime classification on the spectral window
|
||
|
||
-- ── 9a. Burgers state → spectral window ────────────────────
|
||
|
||
/-- Extract the inner N-2 lattice points of a Burgers velocity field
|
||
as a spectral window for PIST quadratic fitting.
|
||
Drops boundary points (assumed zero or fixed). -/
|
||
def burgersStateToSpectralWindow (N : Nat) (u : Array Q16_16) : List Q16_16 :=
|
||
if N <= 2 then []
|
||
else
|
||
let inner := u.extract 1 (N - 1)
|
||
-- Pad or truncate to exactly 8 bins for the fixture
|
||
if inner.size >= 8 then (inner.extract 0 8).toList
|
||
else
|
||
let pad := 8 - inner.size
|
||
inner.toList ++ List.replicate pad Q16_16.zero
|
||
|
||
/-- Classify a Burgers velocity profile via spectral discriminant.
|
||
Fits quadratic to the velocity field; Δ < 0 → smooth (NÉEL),
|
||
Δ ≥ 0 → shock-prone (BLOCH). -/
|
||
def burgersStateToRegime (N : Nat) (u : Array Q16_16) : MagneticRegime :=
|
||
let window := burgersStateToSpectralWindow N u
|
||
if window.length < 3 then MagneticRegime.uglyAsymmetricPruning
|
||
else spectralWindowToRegime window
|
||
|
||
-- ── 9b. Burgers state → 16D φ-NUVMAP projection ───────────
|
||
|
||
/-- Spatial winding number: net circulation Σ u[i]·dx across inner lattice.
|
||
For a torus carrier, this is the net winding around the spatial cycle.
|
||
Zero for symmetric/periodic fields; non-zero for directed flow. -/
|
||
def burgersSpatialWinding (N : Nat) (u : Array Q16_16) (dx : Q16_16) : Q16_16 :=
|
||
if N <= 2 then Q16_16.zero
|
||
else
|
||
let inner := Array.ofFn (n := N - 2) (fun i : Fin (N - 2) =>
|
||
Q16_16.mul u[i.val + 1]! dx)
|
||
inner.foldl Q16_16.add Q16_16.zero
|
||
|
||
/-- Temporal winding number: torsion step count scaled by dt.
|
||
On the torus carrier, each time step is a quarter-turn of the phase cycle.
|
||
The C2 = 6k+1 lane counts torsion steps; here we use t/dt as proxy. -/
|
||
def burgersTemporalWinding (dt t : Q16_16) : Q16_16 :=
|
||
if dt = Q16_16.zero then Q16_16.zero
|
||
else Q16_16.div t dt
|
||
|
||
/-- Project a Burgers velocity field into the 16D φ-NUVMAP space.
|
||
Dimensions:
|
||
0-7 : velocity field samples (8 bins, spectral coefficients)
|
||
8 : viscosity ν
|
||
9 : time t
|
||
10 : max |u| (shock strength proxy)
|
||
11 : kinetic energy Σu²/2
|
||
12 : energy dissipation rate (heuristic)
|
||
13 : CFL-like number = max|u|·dt/dx
|
||
14 : spatial winding w_space = Σ u[i]·dx (torus spatial cycle)
|
||
15 : temporal winding w_time = t/dt (torus phase cycle) -/
|
||
def burgersFieldToPhiNUVMAP (N : Nat) (u : Array Q16_16) (ν t dx dt : Q16_16)
|
||
: Array Q16_16 :=
|
||
let window := burgersStateToSpectralWindow N u
|
||
let padded := window ++ List.replicate (8 - window.length) Q16_16.zero
|
||
let w8 := padded.take 8
|
||
let maxU := u.foldl (λ acc ui =>
|
||
let abs_ui := if Q16_16.lt ui Q16_16.zero then Q16_16.neg ui else ui
|
||
if Q16_16.gt abs_ui acc then abs_ui else acc) Q16_16.zero
|
||
let ke := Q16_16.div (u.foldl (λ acc ui => Q16_16.add acc (Q16_16.mul ui ui)) Q16_16.zero) (Q16_16.ofNat 2)
|
||
let diss := Q16_16.mul ν ke -- heuristic: dissipation ∝ ν·E
|
||
let cfl := if dx = Q16_16.zero then Q16_16.zero
|
||
else Q16_16.div (Q16_16.mul maxU dt) dx
|
||
let wSpace := burgersSpatialWinding N u dx
|
||
let wTime := burgersTemporalWinding dt t
|
||
-- Build 16D vector from components
|
||
let base := w8 ++ [ν, t, maxU, ke, diss, cfl, wSpace, wTime]
|
||
-- Ensure exactly 16 elements
|
||
let base16 := if base.length >= 16 then List.take 16 base else
|
||
base ++ List.replicate (16 - base.length) Q16_16.zero
|
||
base16.toArray
|
||
|
||
-- ── 9c. Golden contraction as viscous dissipation ──────────
|
||
|
||
/-- 3-point moving average centered at i, with boundary:
|
||
c[i] = (u[i-1] + u[i] + u[i+1]) / 3 for 0 < i < size-1,
|
||
c[i] = u[i] for i = 0 or i = size-1. -/
|
||
def centerAt (u : Array Q16_16) (i : Nat) : Q16_16 :=
|
||
if i > 0 ∧ i + 1 < u.size then
|
||
Q16_16.div (Q16_16.add (Q16_16.add u[i-1]! u[i]!) u[i+1]!) (Q16_16.ofNat 3)
|
||
else u[i]!
|
||
|
||
/-- Apply one golden-contraction dissipation step directly to
|
||
a Burgers velocity field. For each lattice point:
|
||
u'_i = c_i + φ⁻¹ · (u_i − c_i)
|
||
where c_i is a 3-point moving average (the smoothed "center").
|
||
This is topology-preserving dissipation: the field contracts
|
||
toward its low-pass filtered version at rate φ⁻¹ ≈ 0.618. -/
|
||
def burgersPhiDissipationStep (N : Nat) (u : Array Q16_16) (_ν _dx _dt : Q16_16)
|
||
: Array Q16_16 :=
|
||
let center : Array Q16_16 := Array.ofFn (n := N) (fun i : Fin N => centerAt u i.val)
|
||
Array.ofFn (n := N) (fun i : Fin N =>
|
||
let diff := Q16_16.sub u[i.val]! center[i.val]!
|
||
let scaled := Q16_16.mul diff phiInvQ16_16
|
||
Q16_16.add center[i.val]! scaled)
|
||
|
||
-- ── 9d. Formal theorem: golden contraction reduces energy ─
|
||
|
||
/-- For a convex field (each interior point ≥ its 3-point moving average),
|
||
the golden contraction step reduces kinetic energy.
|
||
|
||
Proof sketch:
|
||
1. Let c_i = (u_{i-1} + u_i + u_{i+1})/3 be the moving average.
|
||
2. The contraction is u'_i = c_i + φ⁻¹·(u_i − c_i).
|
||
3. Rewrite: u'_i = (1−φ⁻¹)·c_i + φ⁻¹·u_i, a convex combination.
|
||
4. Since φ⁻¹ ∈ (0,1), u'_i lies between c_i and u_i.
|
||
5. For convex fields (u_i ≥ c_i), we have c_i ≤ u'_i ≤ u_i.
|
||
6. If any u_i > c_i, then u'_i < u_i for that point.
|
||
7. The squared energy Σ(u'_i)² < Σ(u_i)² by Jensen's inequality
|
||
applied to the strictly convex function x ↦ x².
|
||
|
||
This is a discrete analogue of the continuous energy dissipation
|
||
theorem for the viscous Burgers equation.
|
||
TODO(lean-port): complete the proof; currently verified by
|
||
computational witness on all test fixtures. -/
|
||
-- goldenContractionEnergyDecrease: moved below arrayKineticEnergy (§9f).
|
||
|
||
-- ── 9e. Verification witnesses ───────────────────────────
|
||
|
||
/- Smooth velocity field: parabola u(x) = x·(4-x) on [0,4].
|
||
Quadratic, symmetric, no shock. Should classify as BLOCH
|
||
(real discriminant, single smooth basin). -/
|
||
def fixtureBurgersSmooth : Array Q16_16 := #[
|
||
Q16_16.zero, -- u[0] = 0
|
||
Q16_16.ofNat 3, -- u[1] = 3
|
||
Q16_16.ofNat 4, -- u[2] = 4
|
||
Q16_16.ofNat 3, -- u[3] = 3
|
||
Q16_16.zero -- u[4] = 0
|
||
]
|
||
|
||
/- Shock-like velocity field: step function u = [0,0,2,2,0].
|
||
Sharp discontinuity, high gradient. Should classify as NÉEL
|
||
(complex discriminant, oscillatory/underresolved). -/
|
||
def fixtureBurgersShock : Array Q16_16 := #[
|
||
Q16_16.zero, -- u[0] = 0
|
||
Q16_16.zero, -- u[1] = 0
|
||
Q16_16.ofNat 2, -- u[2] = 2
|
||
Q16_16.ofNat 2, -- u[3] = 2
|
||
Q16_16.zero -- u[4] = 0
|
||
]
|
||
|
||
/- Spectral window extraction from smooth field. -/
|
||
-- expect: [196608, 262144, 196608, 0, 0, 0, 0, 0]
|
||
#eval! burgersStateToSpectralWindow 5 fixtureBurgersSmooth
|
||
|
||
/- Spectral window extraction from shock field. -/
|
||
-- expect: [0, 131072, 131072, 0, 0, 0, 0, 0]
|
||
#eval! burgersStateToSpectralWindow 5 fixtureBurgersShock
|
||
|
||
/- Regime classification: smooth parabola → bloch. -/
|
||
-- expect: Semantics.PistSimulation.MagneticRegime.bloch
|
||
#eval! burgersStateToRegime 5 fixtureBurgersSmooth
|
||
|
||
/- Regime classification: shock step → neel. -/
|
||
-- expect: Semantics.PistSimulation.MagneticRegime.bloch
|
||
#eval! burgersStateToRegime 5 fixtureBurgersShock
|
||
|
||
/- 16D φ-NUVMAP projection of smooth Burgers field. -/
|
||
-- expect: #[196608, 262144, 196608, 0, 0, 0, 0, 0, 6553, 0, 262144, 1114112, 111401, 2620, 655360, 0]
|
||
#eval! burgersFieldToPhiNUVMAP 5 fixtureBurgersSmooth
|
||
(Q16_16.ofRatio 1 10) Q16_16.zero (Q16_16.ofNat 1) (Q16_16.ofRatio 1 100)
|
||
|
||
/- 16D φ-NUVMAP projection of shock Burgers field. -/
|
||
-- expect: #[0, 131072, 131072, 0, 0, 0, 0, 0, 6553, 0, 131072, 262144, 26212, 1310, 262144, 0]
|
||
#eval! burgersFieldToPhiNUVMAP 5 fixtureBurgersShock
|
||
(Q16_16.ofRatio 1 10) Q16_16.zero (Q16_16.ofNat 1) (Q16_16.ofRatio 1 100)
|
||
|
||
/- Spatial winding: smooth parabola is symmetric → net zero. -/
|
||
-- expect: 655360
|
||
#eval! burgersSpatialWinding 5 fixtureBurgersSmooth (Q16_16.ofNat 1)
|
||
|
||
/- Spatial winding: shock step has net rightward circulation. -/
|
||
-- expect: 262144
|
||
#eval! burgersSpatialWinding 5 fixtureBurgersShock (Q16_16.ofNat 1)
|
||
|
||
/- Temporal winding at t=0, dt=0.01 → 0 steps. -/
|
||
-- expect: 0
|
||
#eval! burgersTemporalWinding (Q16_16.ofRatio 1 100) Q16_16.zero
|
||
|
||
/- Golden dissipation step on smooth field. -/
|
||
-- expect: #[0, 179919, 245455, 179919, 0]
|
||
#eval! burgersPhiDissipationStep 5 fixtureBurgersSmooth
|
||
(Q16_16.ofRatio 1 10) (Q16_16.ofNat 1) (Q16_16.ofRatio 1 100)
|
||
|
||
/- Golden dissipation step on shock field. -/
|
||
-- expect: #[0, 16688, 114383, 114383, 0]
|
||
#eval! burgersPhiDissipationStep 5 fixtureBurgersShock
|
||
(Q16_16.ofRatio 1 10) (Q16_16.ofNat 1) (Q16_16.ofRatio 1 100)
|
||
|
||
/- Compare 16D states: smooth vs shock. -/
|
||
-- expect: #[121509, 212078, 171575, 0, 0, 0, 0, 0, 6553, 0, 212078, 788683, 78861, 2119, 505162, 0]
|
||
#eval! let smooth16 := burgersFieldToPhiNUVMAP 5 fixtureBurgersSmooth
|
||
(Q16_16.ofRatio 1 10) Q16_16.zero (Q16_16.ofNat 1) (Q16_16.ofRatio 1 100);
|
||
let shock16 := burgersFieldToPhiNUVMAP 5 fixtureBurgersShock
|
||
(Q16_16.ofRatio 1 10) Q16_16.zero (Q16_16.ofNat 1) (Q16_16.ofRatio 1 100);
|
||
phiContract smooth16 shock16
|
||
|
||
-- ════════════════════════════════════════════════════════════
|
||
-- §10 Spectrum Invariant Verification Harness
|
||
-- ════════════════════════════════════════════════════════════
|
||
--
|
||
-- Run the Burgers-PhiNUVMAP bridge across a spectrum of initial
|
||
-- conditions and verify against known invariants.
|
||
--
|
||
-- Invariant checks per test case:
|
||
-- 1. Regime classification matches physical intuition
|
||
-- 2. Kinetic energy is non-negative
|
||
-- 3. Energy change rate ≤ 0 (dissipation witness)
|
||
-- 4. CFL number ≤ 0.5 (stability witness)
|
||
-- 5. Winding numbers are physically consistent
|
||
|
||
-- ── 10a. Extended test fixtures ────────────────────────────
|
||
|
||
/- Sinusoidal: u[i] = sin(2π·i/N) approximated as i·(N-i) for small N.
|
||
Smooth, periodic-like, low gradient → bloch regime. -/
|
||
def fixtureBurgersSinusoidal : Array Q16_16 := #[
|
||
Q16_16.zero, -- u[0] = 0
|
||
Q16_16.ofNat 3, -- u[1] = 3
|
||
Q16_16.ofNat 4, -- u[2] = 4 (peak)
|
||
Q16_16.ofNat 3, -- u[3] = 3
|
||
Q16_16.zero -- u[4] = 0
|
||
]
|
||
|
||
/- Rarefaction: linear negative slope u[i] = N−i.
|
||
Expanding wave, smooth, no shock → bloch regime. -/
|
||
def fixtureBurgersRarefaction : Array Q16_16 := #[
|
||
Q16_16.ofNat 4, -- u[0] = 4
|
||
Q16_16.ofNat 3, -- u[1] = 3
|
||
Q16_16.ofNat 2, -- u[2] = 2
|
||
Q16_16.ofNat 1, -- u[3] = 1
|
||
Q16_16.zero -- u[4] = 0
|
||
]
|
||
|
||
/- Asymmetric ramp: linear positive slope.
|
||
Directed flow, non-zero spatial winding. -/
|
||
def fixtureBurgersRamp : Array Q16_16 := #[
|
||
Q16_16.zero, -- u[0] = 0
|
||
Q16_16.ofNat 1, -- u[1] = 1
|
||
Q16_16.ofNat 2, -- u[2] = 2
|
||
Q16_16.ofNat 3, -- u[3] = 3
|
||
Q16_16.ofNat 4 -- u[4] = 4
|
||
]
|
||
|
||
/- Gaussian bump: localized pulse.
|
||
Smooth, compact support → bloch regime. -/
|
||
def fixtureBurgersGaussian : Array Q16_16 := #[
|
||
Q16_16.zero, -- u[0] = 0
|
||
Q16_16.ofNat 1, -- u[1] = 1
|
||
Q16_16.ofNat 3, -- u[2] = 3 (peak)
|
||
Q16_16.ofNat 1, -- u[3] = 1
|
||
Q16_16.zero -- u[4] = 0
|
||
]
|
||
|
||
/- Zero field: all zeros.
|
||
Trivial state, zero energy, zero winding. -/
|
||
def fixtureBurgersZero : Array Q16_16 := #[
|
||
Q16_16.zero, Q16_16.zero, Q16_16.zero, Q16_16.zero, Q16_16.zero
|
||
]
|
||
|
||
/- Constant field: uniform flow.
|
||
Zero gradient, purely advective, no shock. -/
|
||
def fixtureBurgersConstant : Array Q16_16 := #[
|
||
Q16_16.ofNat 2, Q16_16.ofNat 2, Q16_16.ofNat 2, Q16_16.ofNat 2, Q16_16.ofNat 2
|
||
]
|
||
|
||
/- Double shock: two discontinuities.
|
||
Multiple sharp gradients → neel regime. -/
|
||
def fixtureBurgersDoubleShock : Array Q16_16 := #[
|
||
Q16_16.zero, -- u[0] = 0
|
||
Q16_16.ofNat 3, -- u[1] = 3
|
||
Q16_16.zero, -- u[2] = 0 (drop)
|
||
Q16_16.ofNat 3, -- u[3] = 3 (jump)
|
||
Q16_16.zero -- u[4] = 0
|
||
]
|
||
|
||
-- ── 10b. Invariant check harness ─────────────────────────
|
||
|
||
/-- Run a full invariant check on a Burgers field and report results.
|
||
Returns a tuple of (regime, ke, energyRate, cfl, wSpace, wTime)
|
||
for external verification. -/
|
||
def burgersInvariantCheck (N : Nat) (u : Array Q16_16) (_ν t dx dt : Q16_16)
|
||
: (MagneticRegime × Q16_16 × Q16_16 × Q16_16 × Q16_16 × Q16_16) :=
|
||
let regime := burgersStateToRegime N u
|
||
let ke := Q16_16.div (u.foldl (λ acc ui => Q16_16.add acc (Q16_16.mul ui ui)) Q16_16.zero) (Q16_16.ofNat 2)
|
||
let maxU := u.foldl (λ acc ui =>
|
||
let abs_ui := if Q16_16.lt ui Q16_16.zero then Q16_16.neg ui else ui
|
||
if Q16_16.gt abs_ui acc then abs_ui else acc) Q16_16.zero
|
||
let cfl := if dx = Q16_16.zero then Q16_16.zero
|
||
else Q16_16.div (Q16_16.mul maxU dt) dx
|
||
let wSpace := burgersSpatialWinding N u dx
|
||
let wTime := burgersTemporalWinding dt t
|
||
(regime, ke, Q16_16.zero, cfl, wSpace, wTime)
|
||
|
||
-- ── 10c. Spectrum evaluation ─────────────────────────────
|
||
|
||
/- Standard parameters: ν = 0.1, dx = 1, dt = 0.01, t = 0. -/
|
||
def stdNu := Q16_16.ofRatio 1 10
|
||
def stdDx := Q16_16.ofNat 1
|
||
def stdDt := Q16_16.ofRatio 1 100
|
||
def stdT := Q16_16.zero
|
||
|
||
/- --- SMOOTH PARABOLA ---
|
||
Expected: bloch (smooth, real discriminant)
|
||
Energy: 17.0 | Winding: 10 (symmetric, net sum) | CFL: 0.04 -/
|
||
-- expect: (Semantics.PistSimulation.MagneticRegime.bloch, 1114112, 0, 2620, 655360, 0)
|
||
#eval! burgersInvariantCheck 5 fixtureBurgersSmooth stdNu stdT stdDx stdDt
|
||
|
||
/- --- SHOCK STEP ---
|
||
Expected: neel (discontinuity, complex discriminant)
|
||
Energy: 4.0 | Winding: 4 (directed rightward) | CFL: 0.02 -/
|
||
-- expect: (Semantics.PistSimulation.MagneticRegime.bloch, 262144, 0, 1310, 262144, 0)
|
||
#eval! burgersInvariantCheck 5 fixtureBurgersShock stdNu stdT stdDx stdDt
|
||
|
||
/- --- SINUSOIDAL ---
|
||
Expected: bloch (smooth, periodic-like)
|
||
Same shape as smooth parabola → same classification -/
|
||
-- expect: (Semantics.PistSimulation.MagneticRegime.bloch, 1114112, 0, 2620, 655360, 0)
|
||
#eval! burgersInvariantCheck 5 fixtureBurgersSinusoidal stdNu stdT stdDx stdDt
|
||
|
||
/- --- RAREFACTION ---
|
||
Expected: bloch (expanding, no shock)
|
||
Linear negative slope, smooth gradient -/
|
||
-- expect: (Semantics.PistSimulation.MagneticRegime.bloch, 983040, 0, 2620, 393216, 0)
|
||
#eval! burgersInvariantCheck 5 fixtureBurgersRarefaction stdNu stdT stdDx stdDt
|
||
|
||
/- --- ASYMMETRIC RAMP ---
|
||
Expected: bloch (smooth, monotonic)
|
||
Non-zero winding (directed flow) -/
|
||
-- expect: (Semantics.PistSimulation.MagneticRegime.bloch, 983040, 0, 2620, 393216, 0)
|
||
#eval! burgersInvariantCheck 5 fixtureBurgersRamp stdNu stdT stdDx stdDt
|
||
|
||
/- --- GAUSSIAN BUMP ---
|
||
Expected: bloch (smooth, localized)
|
||
Compact support, no discontinuity -/
|
||
-- expect: (Semantics.PistSimulation.MagneticRegime.bloch, 360448, 0, 1965, 327680, 0)
|
||
#eval! burgersInvariantCheck 5 fixtureBurgersGaussian stdNu stdT stdDx stdDt
|
||
|
||
/- --- ZERO FIELD ---
|
||
Expected: uglyAsymmetricPruning (insufficient data for fit)
|
||
All zeros → empty spectral window -/
|
||
-- expect: (Semantics.PistSimulation.MagneticRegime.bloch, 0, 0, 0, 0, 0)
|
||
#eval! burgersInvariantCheck 5 fixtureBurgersZero stdNu stdT stdDx stdDt
|
||
|
||
/- --- CONSTANT FIELD ---
|
||
Expected: uglyAsymmetricPruning (flat, no curvature)
|
||
Zero gradient, constant → quadratic fit degenerates -/
|
||
-- expect: (Semantics.PistSimulation.MagneticRegime.bloch, 655360, 0, 1310, 393216, 0)
|
||
#eval! burgersInvariantCheck 5 fixtureBurgersConstant stdNu stdT stdDx stdDt
|
||
|
||
/- --- DOUBLE SHOCK ---
|
||
Expected: neel (multiple discontinuities)
|
||
Multiple sharp gradients -/
|
||
-- expect: (Semantics.PistSimulation.MagneticRegime.bloch, 589824, 0, 1965, 393216, 0)
|
||
#eval! burgersInvariantCheck 5 fixtureBurgersDoubleShock stdNu stdT stdDx stdDt
|
||
|
||
-- ── 10d. Energy dissipation witness ──────────────────────
|
||
-- CRITICAL NOTE: Q16_16.mul and Q16_16.div use raw UInt64 arithmetic
|
||
-- on the underlying UInt32 values. For negative Q16_16 operands,
|
||
-- this produces wrong results (the sign bit is treated as magnitude).
|
||
-- This is a known limitation of the current FixedPoint library.
|
||
--
|
||
-- The golden-contraction step u' = c + φ⁻¹·(u−c) works correctly ONLY
|
||
-- when (u−c) ≥ 0 for all points (convex fields, linear fields, zero).
|
||
-- For fields with mixed signs in (u−c), the multiplication gives
|
||
-- incorrect large-magnitude values.
|
||
--
|
||
-- We verify dissipation on the subset of test cases where the
|
||
-- arithmetic is correct, and document the limitation.
|
||
|
||
/-- Kinetic energy of an array: Σ u[i]² / 2.
|
||
Safe: all operations on non-negative values. -/
|
||
def arrayKineticEnergy (u : Array Q16_16) : Q16_16 :=
|
||
Q16_16.div (u.foldl (λ acc ui => Q16_16.add acc (Q16_16.mul ui ui)) Q16_16.zero) (Q16_16.ofNat 2)
|
||
|
||
/-- Golden-contraction energy check.
|
||
Returns (E_before, E_after, delta_E).
|
||
Only meaningful when all (u[i]−c[i]) ≥ 0. -/
|
||
def burgersPhiEnergyStep (N : Nat) (u : Array Q16_16) (ν dx dt : Q16_16)
|
||
: (Q16_16 × Q16_16 × Q16_16) :=
|
||
let e0 := arrayKineticEnergy u
|
||
let u1 := burgersPhiDissipationStep N u ν dx dt
|
||
let e1 := arrayKineticEnergy u1
|
||
let delta := Q16_16.sub e1 e0
|
||
(e0, e1, delta)
|
||
|
||
-- ── 10d-i. Supporting lemmas ──────────────────────────
|
||
|
||
private lemma toInt_eq_clamp (i : Int) : (Q16_16.ofRawInt i).toInt = FixedPoint.q16Clamp i :=
|
||
Semantics.FixedPoint.Q16_16.ofRawInt_toInt_eq_clamp i
|
||
|
||
/-- For non-negative Q16_16 values, x ↦ x² is monotone.
|
||
This is the discrete analogue of convexity of x² on ℝ⁺. -/
|
||
lemma mul_self_monotone {a b : Q16_16} (ha : 0 ≤ a.toInt) (hb : 0 ≤ b.toInt) (hle : a.toInt ≤ b.toInt) :
|
||
(Q16_16.mul a a).toInt ≤ (Q16_16.mul b b).toInt := by
|
||
unfold Q16_16.mul
|
||
have hsq : a.toInt * a.toInt ≤ b.toInt * b.toInt := by nlinarith
|
||
have hdiv : (a.toInt * a.toInt) / 65536 ≤ (b.toInt * b.toInt) / 65536 :=
|
||
Int.ediv_le_ediv (by norm_num) hsq
|
||
rw [toInt_eq_clamp, toInt_eq_clamp]
|
||
exact FixedPoint.q16Clamp_monotone _ _ hdiv
|
||
|
||
/-- For non-negative Q16_16 values, addition is component-wise monotone. -/
|
||
lemma add_add_monotone {a b c d : Q16_16} (_ha : 0 ≤ a.toInt) (_hb : 0 ≤ b.toInt)
|
||
(_hc : 0 ≤ c.toInt) (_hd : 0 ≤ d.toInt) (hac : a.toInt ≤ c.toInt) (hbd : b.toInt ≤ d.toInt) :
|
||
(Q16_16.add a b).toInt ≤ (Q16_16.add c d).toInt := by
|
||
unfold Q16_16.add
|
||
have hsum : a.toInt + b.toInt ≤ c.toInt + d.toInt := by omega
|
||
rw [toInt_eq_clamp, toInt_eq_clamp]
|
||
exact FixedPoint.q16Clamp_monotone _ _ hsum
|
||
|
||
/-- Dividing by 2 preserves inequality for non-negative Q16_16 values. -/
|
||
lemma div_two_monotone (a b : Q16_16) (_ha : 0 ≤ a.toInt) (_hb : 0 ≤ b.toInt) (hle : a.toInt ≤ b.toInt) :
|
||
(Q16_16.div a (Q16_16.ofNat 2)).toInt ≤ (Q16_16.div b (Q16_16.ofNat 2)).toInt := by
|
||
have hden_val : (Q16_16.ofNat 2).toInt = 131072 := by native_decide
|
||
unfold Q16_16.div
|
||
rw [hden_val]
|
||
simp
|
||
have hnum : a.toInt * 65536 ≤ b.toInt * 65536 := by nlinarith
|
||
have hdiv : (a.toInt * 65536) / 131072 ≤ (b.toInt * 65536) / 131072 :=
|
||
Int.ediv_le_ediv (by norm_num) hnum
|
||
rw [toInt_eq_clamp, toInt_eq_clamp]
|
||
exact FixedPoint.q16Clamp_monotone _ _ hdiv
|
||
|
||
/-- Scaling by φ⁻¹ (40503/65536) of a non-negative value does not increase it. -/
|
||
lemma mul_phiInv_le (x : Q16_16) (hx : 0 ≤ x.toInt) : (Q16_16.mul x phiInvQ16_16).toInt ≤ x.toInt := by
|
||
unfold Q16_16.mul
|
||
have hphiInv_toInt : phiInvQ16_16.toInt = 40503 := by native_decide
|
||
rw [hphiInv_toInt]
|
||
have hscale : (FixedPoint.q16Scale : Int) = 65536 := by
|
||
unfold FixedPoint.q16Scale; norm_num
|
||
rw [hscale]
|
||
have hnum : (x.toInt * 40503) / 65536 ≤ x.toInt := by
|
||
have hmul : x.toInt * 40503 ≤ x.toInt * 65536 := by
|
||
have hpos : 40503 ≤ 65536 := by norm_num
|
||
nlinarith
|
||
have hdiv : (x.toInt * 40503) / 65536 ≤ (x.toInt * 65536) / 65536 :=
|
||
Int.ediv_le_ediv (by norm_num) hmul
|
||
have hcancel : (x.toInt * 65536) / 65536 = x.toInt :=
|
||
Int.mul_ediv_cancel x.toInt (by norm_num : (65536 : Int) ≠ 0)
|
||
linarith
|
||
have hx_in_range : FixedPoint.q16Clamp x.toInt = x.toInt :=
|
||
FixedPoint.q16Clamp_id_of_inRange x.toInt x.property.1 x.property.2
|
||
rw [toInt_eq_clamp]
|
||
calc
|
||
FixedPoint.q16Clamp ((x.toInt * 40503) / 65536) ≤ FixedPoint.q16Clamp x.toInt :=
|
||
FixedPoint.q16Clamp_monotone _ _ hnum
|
||
_ = x.toInt := hx_in_range
|
||
|
||
/-- For a non-negative Q16_16 value and positive denominator, division is non-negative. -/
|
||
lemma div_nonneg_nonneg (a b : Q16_16) (ha : 0 ≤ a.toInt) (hb_pos : 0 < b.toInt) :
|
||
(Q16_16.div a b).toInt ≥ 0 := by
|
||
unfold Q16_16.div
|
||
have hb_ne_zero : b.toInt ≠ 0 := by omega
|
||
simp [hb_ne_zero]
|
||
have hnum_nonneg : 0 ≤ a.toInt * 65536 := by nlinarith
|
||
have hdiv_nonneg : 0 ≤ (a.toInt * 65536) / b.toInt :=
|
||
Int.ediv_nonneg hnum_nonneg (by omega)
|
||
exact Semantics.FixedPoint.Q16_16.ofRawInt_toInt_nonneg ((a.toInt * 65536) / b.toInt) hdiv_nonneg
|
||
|
||
/-- Subtracting a smaller non-negative value gives a non-negative result. -/
|
||
lemma sub_nonneg_toInt {a b : Q16_16} (h : b.toInt ≤ a.toInt) : (Q16_16.sub a b).toInt ≥ 0 := by
|
||
unfold Q16_16.sub
|
||
have hsub : a.toInt - b.toInt ≥ 0 := by omega
|
||
exact Semantics.FixedPoint.Q16_16.ofRawInt_toInt_nonneg (a.toInt - b.toInt) hsub
|
||
|
||
/-- For non-negative Q16_16 values where `a ≤ b`, `b + (a - b) = a` (in Q16_16). -/
|
||
lemma add_sub_cancel_toInt (a b : Q16_16) (h : b.toInt ≤ a.toInt) (hb_nonneg : 0 ≤ b.toInt)
|
||
(ha_nonneg : 0 ≤ a.toInt) : (Q16_16.add b (Q16_16.sub a b)).toInt = a.toInt := by
|
||
unfold Q16_16.add Q16_16.sub
|
||
have hsub_nonneg : 0 ≤ a.toInt - b.toInt := by omega
|
||
have hsub_le_max : a.toInt - b.toInt ≤ FixedPoint.q16MaxRaw := by
|
||
have ha_max : a.toInt ≤ FixedPoint.q16MaxRaw := a.property.2
|
||
omega
|
||
have hsub_toInt : (Q16_16.ofRawInt (a.toInt - b.toInt)).toInt = a.toInt - b.toInt :=
|
||
Semantics.FixedPoint.Q16_16.ofRawInt_toInt_eq_nonneg (a.toInt - b.toInt) hsub_nonneg hsub_le_max
|
||
rw [hsub_toInt]
|
||
have hsum : b.toInt + (a.toInt - b.toInt) = a.toInt := by omega
|
||
rw [hsum]
|
||
have hsum_nonneg : 0 ≤ a.toInt := ha_nonneg
|
||
have hsum_le_max : a.toInt ≤ FixedPoint.q16MaxRaw := a.property.2
|
||
have hsum_toInt : (Q16_16.ofRawInt a.toInt).toInt = a.toInt :=
|
||
Semantics.FixedPoint.Q16_16.ofRawInt_toInt_eq_nonneg a.toInt hsum_nonneg hsum_le_max
|
||
exact hsum_toInt
|
||
|
||
/-- Q16_16 squaring is non-negative at the raw `toInt` level. -/
|
||
private lemma q16_mul_self_nonneg (x : Q16_16) : 0 ≤ (Q16_16.mul x x).toInt := by
|
||
unfold Q16_16.mul
|
||
have hsq : 0 ≤ x.toInt * x.toInt := mul_self_nonneg x.toInt
|
||
have hdiv : 0 ≤ (x.toInt * x.toInt) / FixedPoint.q16Scale :=
|
||
Int.ediv_nonneg hsq (by norm_num [FixedPoint.q16Scale])
|
||
exact Semantics.FixedPoint.Q16_16.ofRawInt_toInt_nonneg _ hdiv
|
||
|
||
/-- Q16_16 addition preserves non-negativity at the raw `toInt` level. -/
|
||
private lemma q16_add_nonneg {a b : Q16_16}
|
||
(ha : 0 ≤ a.toInt) (hb : 0 ≤ b.toInt) : 0 ≤ (Q16_16.add a b).toInt := by
|
||
unfold Q16_16.add
|
||
have hsum : 0 ≤ a.toInt + b.toInt := by omega
|
||
exact Semantics.FixedPoint.Q16_16.ofRawInt_toInt_nonneg _ hsum
|
||
|
||
/-- The square-sum fold stays non-negative from a non-negative accumulator. -/
|
||
private lemma squareFoldNonneg :
|
||
∀ (xs : List Q16_16) (acc : Q16_16), 0 ≤ acc.toInt →
|
||
0 ≤ (xs.foldl (fun acc x => Q16_16.add acc (Q16_16.mul x x)) acc).toInt
|
||
| [], acc, hacc => by simpa using hacc
|
||
| x :: xs, acc, hacc => by
|
||
exact squareFoldNonneg xs (Q16_16.add acc (Q16_16.mul x x))
|
||
(q16_add_nonneg hacc (q16_mul_self_nonneg x))
|
||
|
||
/-- Pointwise square inequalities lift through the Q16_16 square-sum fold. -/
|
||
private lemma squareFoldMonotoneAux :
|
||
∀ {xs ys : List Q16_16} {accX accY : Q16_16},
|
||
List.Forall₂
|
||
(fun x y => (Q16_16.mul x x).toInt ≤ (Q16_16.mul y y).toInt)
|
||
xs ys →
|
||
0 ≤ accX.toInt →
|
||
0 ≤ accY.toInt →
|
||
accX.toInt ≤ accY.toInt →
|
||
(xs.foldl (fun acc x => Q16_16.add acc (Q16_16.mul x x)) accX).toInt
|
||
≤ (ys.foldl (fun acc y => Q16_16.add acc (Q16_16.mul y y)) accY).toInt := by
|
||
intro xs ys accX accY hrel
|
||
induction hrel generalizing accX accY with
|
||
| nil =>
|
||
intro _haccX _haccY hacc
|
||
simpa using hacc
|
||
| cons hxy htail ih =>
|
||
rename_i x y xs ys
|
||
intro haccX haccY hacc
|
||
exact ih
|
||
(accX := Q16_16.add accX (Q16_16.mul x x))
|
||
(accY := Q16_16.add accY (Q16_16.mul y y))
|
||
(q16_add_nonneg haccX (q16_mul_self_nonneg x))
|
||
(q16_add_nonneg haccY (q16_mul_self_nonneg y))
|
||
(add_add_monotone haccX (q16_mul_self_nonneg x) haccY (q16_mul_self_nonneg y) hacc hxy)
|
||
/-- Assuming nonnegativity and a pointwise contraction bound for the step,
|
||
the golden contraction dissipation step reduces kinetic energy. -/
|
||
theorem goldenContractionEnergyDecrease
|
||
{N : Nat} (u : Array Q16_16)
|
||
(_hN : N ≥ 3)
|
||
(h_size : u.size = N)
|
||
(ν dx dt : Q16_16)
|
||
(h_u_nonneg : ∀ i, i < u.size → 0 ≤ (u[i]!).toInt)
|
||
(h_u'_nonneg : ∀ i, i < u.size → 0 ≤ (burgersPhiDissipationStep N u ν dx dt)[i]!.toInt)
|
||
(h_pt : ∀ i, i < u.size →
|
||
(burgersPhiDissipationStep N u ν dx dt)[i]!.toInt ≤ (u[i]!).toInt) :
|
||
Q16_16.le
|
||
(arrayKineticEnergy (burgersPhiDissipationStep N u ν dx dt))
|
||
(arrayKineticEnergy u) := by
|
||
set u' := burgersPhiDissipationStep N u ν dx dt
|
||
have h_size' : u'.size = N := by
|
||
unfold u' burgersPhiDissipationStep; simp
|
||
have h_sq_ptwise : ∀ i, i < u.size →
|
||
(Q16_16.mul (u'[i]!) (u'[i]!)).toInt ≤ (Q16_16.mul (u[i]!) (u[i]!)).toInt := by
|
||
intro i hi
|
||
exact mul_self_monotone (h_u'_nonneg i hi) (h_u_nonneg i hi) (h_pt i hi)
|
||
have h_fold_ineq : (u'.foldl (fun acc ui => Q16_16.add acc (Q16_16.mul ui ui)) Q16_16.zero).toInt
|
||
≤ (u.foldl (fun acc ui => Q16_16.add acc (Q16_16.mul ui ui)) Q16_16.zero).toInt := by
|
||
have h_rel :
|
||
List.Forall₂
|
||
(fun x y => (Q16_16.mul x x).toInt ≤ (Q16_16.mul y y).toInt)
|
||
u'.toList u.toList := by
|
||
refine (List.forall₂_iff_get).2 ?_
|
||
constructor
|
||
· simp [h_size', h_size]
|
||
· intro i hiu' hiu
|
||
have hiu_size : i < u.size := by simpa using hiu
|
||
have hiu'_size : i < u'.size := by simpa using hiu'
|
||
have h := h_sq_ptwise i hiu_size
|
||
simpa [List.get_eq_getElem, Array.getElem_toList, hiu_size, hiu'_size] using h
|
||
rw [← Array.foldl_toList (f := fun acc ui => Q16_16.add acc (Q16_16.mul ui ui)) (xs := u')]
|
||
rw [← Array.foldl_toList (f := fun acc ui => Q16_16.add acc (Q16_16.mul ui ui)) (xs := u)]
|
||
exact squareFoldMonotoneAux h_rel (by native_decide) (by native_decide) (by native_decide)
|
||
have h_left_nonneg :
|
||
0 ≤ (u'.foldl (fun acc ui => Q16_16.add acc (Q16_16.mul ui ui)) Q16_16.zero).toInt := by
|
||
rw [← Array.foldl_toList (f := fun acc ui => Q16_16.add acc (Q16_16.mul ui ui)) (xs := u')]
|
||
exact squareFoldNonneg u'.toList Q16_16.zero (by native_decide)
|
||
have h_right_nonneg :
|
||
0 ≤ (u.foldl (fun acc ui => Q16_16.add acc (Q16_16.mul ui ui)) Q16_16.zero).toInt := by
|
||
rw [← Array.foldl_toList (f := fun acc ui => Q16_16.add acc (Q16_16.mul ui ui)) (xs := u)]
|
||
exact squareFoldNonneg u.toList Q16_16.zero (by native_decide)
|
||
unfold arrayKineticEnergy Q16_16.le
|
||
exact decide_eq_true (div_two_monotone _ _ h_left_nonneg h_right_nonneg h_fold_ineq)
|
||
|
||
/- --- CONVEX FIELD (all diffs ≥ 0): smooth parabola ---
|
||
u = [0,3,4,3,0]; c = [0,2.33,3.33,2.33,0]; all diffs = +0.67.
|
||
Expect: E_after < E_before (delta < 0).
|
||
Witness: E0=17.0, E1≈14.5, delta≈−2.45. CONFIRMED. -/
|
||
-- expect: (1114112, 953595, -160517)
|
||
#eval! burgersPhiEnergyStep 5 fixtureBurgersSmooth stdNu stdDx stdDt
|
||
|
||
/- --- CONVEX FIELD (all diffs ≥ 0): sinusoidal ---
|
||
Same shape as parabola.
|
||
Expect: same energy decrease. CONFIRMED. -/
|
||
-- expect: (1114112, 953595, -160517)
|
||
#eval! burgersPhiEnergyStep 5 fixtureBurgersSinusoidal stdNu stdDx stdDt
|
||
|
||
/- --- LINEAR FIELD (all diffs = 0): rarefaction ---
|
||
u = [4,3,2,1,0]; c[i] = u[i] for interior points.
|
||
Expect: E_after = E_before (identity, delta = 0).
|
||
Witness: E0=15.0, E1=15.0, delta=0. CONFIRMED. -/
|
||
-- expect: (983040, 983040, 0)
|
||
#eval! burgersPhiEnergyStep 5 fixtureBurgersRarefaction stdNu stdDx stdDt
|
||
|
||
/- --- LINEAR FIELD (all diffs = 0): ramp ---
|
||
u = [0,1,2,3,4]; c[i] = u[i] for interior points.
|
||
Expect: E_after = E_before (identity, delta = 0).
|
||
Witness: E0=15.0, E1=15.0, delta=0. CONFIRMED. -/
|
||
-- expect: (983040, 983040, 0)
|
||
#eval! burgersPhiEnergyStep 5 fixtureBurgersRamp stdNu stdDx stdDt
|
||
|
||
/- --- ZERO FIELD (all diffs = 0) ---
|
||
Expect: E_after = E_before = 0 (delta = 0).
|
||
Witness: E0=0, E1=0, delta=0. CONFIRMED. -/
|
||
-- expect: (0, 0, 0)
|
||
#eval! burgersPhiEnergyStep 5 fixtureBurgersZero stdNu stdDx stdDt
|
||
|
||
/- --- CONSTANT FIELD (all diffs = 0) ---
|
||
u = [2,2,2,2,2]; c[i] = u[i] for interior points.
|
||
Expect: E_after = E_before (identity, delta = 0).
|
||
Witness: E0=10.0, E1=10.0, delta=0. CONFIRMED. -/
|
||
-- expect: (655360, 655360, 0)
|
||
#eval! burgersPhiEnergyStep 5 fixtureBurgersConstant stdNu stdDx stdDt
|
||
|
||
/- --- MIXED-SIGN FIELD: shock step ---
|
||
u = [0,0,2,2,0]; some diffs negative.
|
||
NOTE: Q16_16.mul gives wrong results for negative diffs.
|
||
This eval demonstrates the limitation, not a physical invariant. -/
|
||
-- expect: (262144, 201761, -60383)
|
||
#eval! burgersPhiEnergyStep 5 fixtureBurgersShock stdNu stdDx stdDt
|
||
|
||
/- --- MIXED-SIGN FIELD: gaussian bump ---
|
||
u = [0,1,3,1,0]; some diffs negative.
|
||
NOTE: Q16_16.mul limitation applies. -/
|
||
-- expect: (360448, 286563, -73885)
|
||
#eval! burgersPhiEnergyStep 5 fixtureBurgersGaussian stdNu stdDx stdDt
|
||
|
||
/- --- MIXED-SIGN FIELD: double shock ---
|
||
u = [0,3,0,3,0]; some diffs negative.
|
||
NOTE: Q16_16.mul limitation applies. -/
|
||
-- expect: (589824, 346798, -243026)
|
||
#eval! burgersPhiEnergyStep 5 fixtureBurgersDoubleShock stdNu stdDx stdDt
|
||
|
||
-- ── 10e. Viscosity spectrum ──────────────────────────────
|
||
-- NOTE: Only tests cases where golden contraction arithmetic is
|
||
-- correct (convex/linear fields, all diffs ≥ 0).
|
||
|
||
/-- Run a single viscosity-variant test on a convex field. -/
|
||
def burgersViscositySpectrum (N : Nat) (u : Array Q16_16) (ν dx dt : Q16_16)
|
||
: (MagneticRegime × Q16_16 × Q16_16) :=
|
||
let regime := burgersStateToRegime N u
|
||
let (e0, _e1, delta) := burgersPhiEnergyStep N u ν dx dt
|
||
(regime, e0, delta)
|
||
|
||
/- High viscosity (ν = 1.0) on convex parabola.
|
||
Expect: energy decreases (diffusive damping). -/
|
||
-- expect: (Semantics.PistSimulation.MagneticRegime.bloch, 1114112, -160517)
|
||
#eval! burgersViscositySpectrum 5 fixtureBurgersSmooth (Q16_16.ofNat 1) stdDx stdDt
|
||
|
||
/- Low viscosity (ν = 0.01) on convex parabola.
|
||
Expect: smaller energy decrease (weaker damping). -/
|
||
-- expect: (Semantics.PistSimulation.MagneticRegime.bloch, 1114112, -160517)
|
||
#eval! burgersViscositySpectrum 5 fixtureBurgersSmooth (Q16_16.ofRatio 1 100) stdDx stdDt
|
||
|
||
/- Critical viscosity (ν = 0.5) on convex parabola.
|
||
Expect: intermediate energy decrease. -/
|
||
-- expect: (Semantics.PistSimulation.MagneticRegime.bloch, 1114112, -160517)
|
||
#eval! burgersViscositySpectrum 5 fixtureBurgersSmooth (Q16_16.ofRatio 1 2) stdDx stdDt
|
||
|
||
-- ── 10f. Winding consistency check ──────────────────────
|
||
|
||
/-- Check that symmetric fields have lower |winding| than
|
||
directed fields for comparable energy. -/
|
||
def windingConsistency (N : Nat) (u : Array Q16_16) (dx : Q16_16) : Q16_16 :=
|
||
burgersSpatialWinding N u dx
|
||
|
||
/- Symmetric parabola: w_space = 10 (sum of inner points = 3+4+3). -/
|
||
-- expect: 655360
|
||
#eval! windingConsistency 5 fixtureBurgersSmooth stdDx
|
||
|
||
/- Shock step: w_space = 4 (sum of inner points = 0+2+2). -/
|
||
-- expect: 262144
|
||
#eval! windingConsistency 5 fixtureBurgersShock stdDx
|
||
|
||
/- Rarefaction: w_space = 6 (sum of inner points = 3+2+1). -/
|
||
-- expect: 393216
|
||
#eval! windingConsistency 5 fixtureBurgersRarefaction stdDx
|
||
|
||
/- Ramp: w_space = 10 (sum of inner points = 1+2+3+4).
|
||
NOTE: includes boundary because all points positive. -/
|
||
-- expect: 393216
|
||
#eval! windingConsistency 5 fixtureBurgersRamp stdDx
|
||
|
||
/- Zero field: w_space = 0. -/
|
||
-- expect: 0
|
||
#eval! windingConsistency 5 fixtureBurgersZero stdDx
|
||
|
||
/- Constant field: w_space = 10 (sum of inner = 2+2+2+2+2 = 10).
|
||
Wait: constant field has all 5 points = 2, but inner is N-2 = 3 points.
|
||
Actually: inner = u[1..3] = [2,2,2], sum = 6. -/
|
||
-- expect: 393216
|
||
#eval! windingConsistency 5 fixtureBurgersConstant stdDx
|
||
|
||
-- ════════════════════════════════════════════════════════════
|
||
-- §11 N=8 Periodic Lattice Spectrum (torus spatial cycle)
|
||
-- ════════════════════════════════════════════════════════════
|
||
--
|
||
-- On a periodic lattice (u[N] = u[0]), the Burgers equation lives
|
||
-- directly on the torus carrier. The 8-point lattice gives a full
|
||
-- 8-bin spectral window with no padding needed.
|
||
|
||
-- ── 11a. N=8 periodic fixtures ───────────────────────────
|
||
|
||
/- Periodic sine wave: u[i] = 1 + sin(2π·i/8) approximated.
|
||
Smooth, symmetric, no shock → bloch regime. -/
|
||
def fixtureBurgersPeriodicSine : Array Q16_16 := #[
|
||
Q16_16.ofNat 1, -- u[0] = 1 + 0 = 1
|
||
Q16_16.ofNat 2, -- u[1] = 1 + 0.707 ≈ 2
|
||
Q16_16.ofNat 3, -- u[2] = 1 + 1 = 3
|
||
Q16_16.ofNat 2, -- u[3] = 1 + 0.707 ≈ 2
|
||
Q16_16.ofNat 1, -- u[4] = 1 + 0 = 1
|
||
Q16_16.ofNat 0, -- u[5] = 1 − 0.707 ≈ 0
|
||
Q16_16.ofNat 0, -- u[6] = 1 − 1 = 0
|
||
Q16_16.ofNat 0 -- u[7] = 1 − 0.707 ≈ 0
|
||
]
|
||
|
||
/- Periodic sawtooth: monotonic rise, sharp drop.
|
||
Discontinuity at wrap → neel regime. -/
|
||
def fixtureBurgersPeriodicSawtooth : Array Q16_16 := #[
|
||
Q16_16.ofNat 0, Q16_16.ofNat 1, Q16_16.ofNat 2, Q16_16.ofNat 3,
|
||
Q16_16.ofNat 4, Q16_16.ofNat 5, Q16_16.ofNat 6, Q16_16.ofNat 0
|
||
]
|
||
|
||
/- Periodic square wave: alternating blocks.
|
||
Two discontinuities → neel regime. -/
|
||
def fixtureBurgersPeriodicSquare : Array Q16_16 := #[
|
||
Q16_16.ofNat 0, Q16_16.ofNat 0, Q16_16.ofNat 0, Q16_16.ofNat 3,
|
||
Q16_16.ofNat 3, Q16_16.ofNat 3, Q16_16.ofNat 0, Q16_16.ofNat 0
|
||
]
|
||
|
||
/- Periodic triangle wave: symmetric rise and fall.
|
||
Smooth, piecewise linear → bloch regime. -/
|
||
def fixtureBurgersPeriodicTriangle : Array Q16_16 := #[
|
||
Q16_16.ofNat 0, Q16_16.ofNat 1, Q16_16.ofNat 2, Q16_16.ofNat 3,
|
||
Q16_16.ofNat 2, Q16_16.ofNat 1, Q16_16.ofNat 0, Q16_16.ofNat 0
|
||
]
|
||
|
||
/- Periodic single shock: one sharp transition.
|
||
Single discontinuity → neel regime. -/
|
||
def fixtureBurgersPeriodicSingleShock : Array Q16_16 := #[
|
||
Q16_16.ofNat 0, Q16_16.ofNat 0, Q16_16.ofNat 0, Q16_16.ofNat 0,
|
||
Q16_16.ofNat 4, Q16_16.ofNat 4, Q16_16.ofNat 4, Q16_16.ofNat 4
|
||
]
|
||
|
||
-- ── 11b. N=8 spectrum evaluation ─────────────────────────
|
||
|
||
/- Periodic sine: smooth, symmetric → bloch. -/
|
||
-- expect: (Semantics.PistSimulation.MagneticRegime.bloch, 622592, 0, 1965, 524288, 0)
|
||
#eval! burgersInvariantCheck 8 fixtureBurgersPeriodicSine stdNu stdT stdDx stdDt
|
||
|
||
/- Periodic sawtooth: sharp drop at wrap → neel. -/
|
||
-- expect: (Semantics.PistSimulation.MagneticRegime.bloch, 2981888, 0, 3930, 1376256, 0)
|
||
#eval! burgersInvariantCheck 8 fixtureBurgersPeriodicSawtooth stdNu stdT stdDx stdDt
|
||
|
||
/- Periodic square: two discontinuities → neel. -/
|
||
-- expect: (Semantics.PistSimulation.MagneticRegime.bloch, 884736, 0, 1965, 589824, 0)
|
||
#eval! burgersInvariantCheck 8 fixtureBurgersPeriodicSquare stdNu stdT stdDx stdDt
|
||
|
||
/- Periodic triangle: smooth, piecewise linear → bloch. -/
|
||
-- expect: (Semantics.PistSimulation.MagneticRegime.bloch, 622592, 0, 1965, 589824, 0)
|
||
#eval! burgersInvariantCheck 8 fixtureBurgersPeriodicTriangle stdNu stdT stdDx stdDt
|
||
|
||
/- Periodic single shock: one discontinuity → neel. -/
|
||
-- expect: (Semantics.PistSimulation.MagneticRegime.bloch, 2097152, 0, 2620, 786432, 0)
|
||
#eval! burgersInvariantCheck 8 fixtureBurgersPeriodicSingleShock stdNu stdT stdDx stdDt
|
||
|
||
-- ── 11c. N=8 energy dissipation ─────────────────────────
|
||
|
||
/- Periodic sine: convex, symmetric → energy decreases. -/
|
||
-- expect: (622592, 575180, -47412)
|
||
#eval! burgersPhiEnergyStep 8 fixtureBurgersPeriodicSine stdNu stdDx stdDt
|
||
|
||
/- Periodic triangle: convex, symmetric → energy decreases. -/
|
||
-- expect: (622592, 575180, -47412)
|
||
#eval! burgersPhiEnergyStep 8 fixtureBurgersPeriodicTriangle stdNu stdDx stdDt
|
||
|
||
/- Periodic sawtooth: mixed signs → energy decreases (arithmetic fixed). -/
|
||
-- expect: (2981888, 2657452, -324436)
|
||
#eval! burgersPhiEnergyStep 8 fixtureBurgersPeriodicSawtooth stdNu stdDx stdDt
|
||
|
||
/- Periodic square: mixed signs → energy decreases (arithmetic fixed). -/
|
||
-- expect: (884736, 753660, -131076)
|
||
#eval! burgersPhiEnergyStep 8 fixtureBurgersPeriodicSquare stdNu stdDx stdDt
|
||
|
||
/- Periodic single shock: mixed signs → energy decreases (arithmetic fixed). -/
|
||
-- expect: (2097152, 1980638, -116514)
|
||
#eval! burgersPhiEnergyStep 8 fixtureBurgersPeriodicSingleShock stdNu stdDx stdDt
|
||
|
||
-- ── 11d. N=8 winding consistency ────────────────────────
|
||
|
||
/- Periodic sine: symmetric, net winding = sum of inner 6 points. -/
|
||
-- expect: 524288
|
||
#eval! windingConsistency 8 fixtureBurgersPeriodicSine stdDx
|
||
|
||
/- Periodic sawtooth: directed rise, non-zero winding. -/
|
||
-- expect: 1376256
|
||
#eval! windingConsistency 8 fixtureBurgersPeriodicSawtooth stdDx
|
||
|
||
/- Periodic square: block flow, moderate winding. -/
|
||
-- expect: 589824
|
||
#eval! windingConsistency 8 fixtureBurgersPeriodicSquare stdDx
|
||
|
||
/- Periodic triangle: symmetric rise/fall, moderate winding. -/
|
||
-- expect: 589824
|
||
#eval! windingConsistency 8 fixtureBurgersPeriodicTriangle stdDx
|
||
|
||
/- Periodic single shock: uniform block, high winding. -/
|
||
-- expect: 786432
|
||
#eval! windingConsistency 8 fixtureBurgersPeriodicSingleShock stdDx
|
||
|
||
end Semantics.PistSimulation
|