-- SilverSight.PIST.Spectral — minimal fixed-point spectral feature extraction
--
-- Ports only the surface needed for the RRC 250-equation classifier:
--   • integer square root
--   • power iteration for dominant eigenvalue
--   • SpectralProfile record
--   • computeSpectral over an 8×8 Int matrix
--
-- This is a stripped-down port of Research Stack's Semantics.PIST.Spectral.
-- Photonic extensions, tactic-family heuristics, and Kubelka-Munk material
-- are intentionally omitted to keep SilverSight focused on the RRC gate.

import SilverSight.FixedPoint

namespace SilverSight.PIST.Spectral

open SilverSight.FixedPoint
open SilverSight.FixedPoint.Q16_16

-- ─────────────────────────────────────────────────────────────────────────────
-- §1  Integer square root
-- ─────────────────────────────────────────────────────────────────────────────

/-- Integer square root via Newton's method. Returns floor(√n). -/
def isqrt (n : Int) : Int :=
  if n ≤ 0 then 0
  else
    let rec loop (x : Int) (fuel : Nat) : Int :=
      match fuel with
      | 0     => x
      | f + 1 =>
        let x' := (x + n / x) / 2
        if x' ≥ x then x else loop x' f
    loop (n / 2 + 1) 64

-- ─────────────────────────────────────────────────────────────────────────────
-- §2  Matrix-vector helpers (pure Int arithmetic)
-- ─────────────────────────────────────────────────────────────────────────────

@[inline]
private def getEntry (mat : Array (Array Int)) (i j : Nat) : Int :=
  mat.getD i #[] |>.getD j 0

private def rowSum (mat : Array (Array Int)) (i n : Nat) : Int :=
  (List.range n).foldl (fun acc j => acc + getEntry mat i j) 0

private def symmetrize (mat : Array (Array Int)) (n : Nat) : Array (Array Int) :=
  Array.ofFn (n := n) fun i =>
    Array.ofFn (n := n) fun j =>
      (getEntry mat i.val j.val + getEntry mat j.val i.val) / 2

private def buildLaplacian (sym : Array (Array Int)) (n : Nat) : Array (Array Int) :=
  Array.ofFn (n := n) fun i =>
    let deg := rowSum sym i.val n
    Array.ofFn (n := n) fun j =>
      if i.val = j.val then deg else -(getEntry sym i.val j.val)

private def buildATA (mat : Array (Array Int)) (n : Nat) : Array (Array Int) :=
  Array.ofFn (n := n) fun i =>
    Array.ofFn (n := n) fun j =>
      (List.range n).foldl (fun acc k =>
        acc + getEntry mat k i.val * getEntry mat k j.val) 0

-- ─────────────────────────────────────────────────────────────────────────────
-- §3  Power iteration (Q16_16 fixed-point)
-- ─────────────────────────────────────────────────────────────────────────────

private def normSqRaw (v : Array Q16_16) : Int :=
  v.foldl (fun acc x => acc + x.toInt * x.toInt) 0

private def matVecMul (mat : Array (Array Int)) (n : Nat) (v : Array Q16_16) : Array Q16_16 :=
  Array.ofFn (n := n) fun i =>
    let s : Int := (List.range n).foldl (fun acc j =>
      acc + getEntry mat i.val j * (v.getD j zero).toInt) 0
    ofRawInt s

/-- Dominant eigenvalue of a square Int matrix via power iteration. -/
def powerIteration (mat : Array (Array Int)) (maxIter : Nat := 100) : Q16_16 :=
  let n := mat.size
  if n = 0 then zero
  else
    let initVal := ofRawInt ((q16Scale : Int) / (n : Int))
    let v₀ : Array Q16_16 := Array.replicate n initVal
    let rec iterate (v : Array Q16_16) (fuel : Nat) : Array Q16_16 :=
      match fuel with
      | 0     => v
      | f + 1 =>
        let mv  := matVecMul mat n v
        let nm2 := normSqRaw mv
        if nm2 ≤ 0 then v
        else
          let nm := isqrt nm2
          if nm = 0 then v
          else
            let vn := mv.map (fun x => ofRawInt (x.toInt * q16Scale / nm))
            iterate vn f
    let vFinal := iterate v₀ maxIter
    let mv  := matVecMul mat n vFinal
    let num : Int := (List.range n).foldl (fun acc i =>
      acc + (vFinal.getD i zero).toInt * (mv.getD i zero).toInt) 0
    let den : Int := normSqRaw vFinal
    if den ≤ 0 then zero
    else ofRawInt (num * q16Scale / den)

-- ─────────────────────────────────────────────────────────────────────────────
-- §4  Spectral profile
-- ─────────────────────────────────────────────────────────────────────────────

structure SpectralProfile where
  matrix_size             : Q16_16
  rank                    : Q16_16
  spectral_gap            : Q16_16
  density                 : Q16_16
  trace_val               : Q16_16
  frobenius_norm          : Q16_16
  laplacian_zero_count    : Q16_16
  adjacency_eigenvalue_max : Q16_16
  laplacian_eigenvalue_max : Q16_16
  singular_value_max      : Q16_16
  deriving Repr

def emptyProfile : SpectralProfile :=
  { matrix_size             := zero
    rank                    := zero
    spectral_gap            := zero
    density                 := zero
    trace_val               := zero
    frobenius_norm          := zero
    laplacian_zero_count    := zero
    adjacency_eigenvalue_max := zero
    laplacian_eigenvalue_max := zero
    singular_value_max      := zero }

-- ─────────────────────────────────────────────────────────────────────────────
-- §5  computeSpectral
-- ─────────────────────────────────────────────────────────────────────────────

/-- Compute the full spectral profile of a proof transition matrix.
    Pure Int + Q16_16 arithmetic; no Float. -/
def computeSpectral (mat : Array (Array Int)) : SpectralProfile :=
  let n := mat.size
  if n = 0 then emptyProfile
  else
    let sym := symmetrize mat n
    let lap := buildLaplacian sym n

    let evMax := powerIteration sym

    -- Second eigenvalue via shift-deflation (0.9 * ev_max on the diagonal)
    let shiftAmt : Int := (evMax.toInt * 58982) / q16Scale
    let shiftedMat : Array (Array Int) :=
      Array.ofFn (n := n) fun i =>
        Array.ofFn (n := n) fun j =>
          let base := getEntry sym i.val j.val
          if i.val = j.val then base - shiftAmt else base
    let evShift := powerIteration shiftedMat
    let evSecond : Q16_16 :=
      if evShift.toInt < evMax.toInt
      then ofRawInt (Int.natAbs (evMax.toInt - evShift.toInt) : Int)
      else evMax
    let spectralGap := sub evMax evSecond

    let lapMax := powerIteration lap

    let total : Int := mat.foldl (fun acc row => acc + row.foldl (· + ·) 0) 0
    let nSqI : Int := (n * n : Nat)
    let density := if nSqI ≤ 0 then zero else ofRawInt (total * q16Scale / nSqI)

    let traceInt : Int := (List.range n).foldl (fun acc i => acc + getEntry mat i i) 0
    let traceVal := ofRawInt traceInt

    let frobSq : Int := mat.foldl (fun acc row => acc + row.foldl (fun a c => a + c * c) 0) 0
    let frobNorm := ofRawInt (isqrt frobSq)

    let rankVal : Int := mat.foldl (fun acc row =>
      acc + if row.any (· ≠ 0) then 1 else 0) 0

    let lapZero : Int := (List.range n).foldl (fun acc i =>
      let rs := rowSum mat i n
      let d  := getEntry mat i i
      acc + if rs - d = 0 then 1 else 0) 0

    let ataM := buildATA mat n
    let ataEv := powerIteration ataM
    let svMax := if ataEv.toInt > 0 then sqrt ataEv else zero

    { matrix_size             := ofRawInt (n : Int)
      rank                    := ofRawInt rankVal
      spectral_gap            := spectralGap
      density                 := density
      trace_val               := traceVal
      frobenius_norm          := frobNorm
      laplacian_zero_count    := ofRawInt lapZero
      adjacency_eigenvalue_max := evMax
      laplacian_eigenvalue_max := lapMax
      singular_value_max      := svMax }

-- ─────────────────────────────────────────────────────────────────────────────
-- §6  Witnesses
-- ─────────────────────────────────────────────────────────────────────────────

#eval isqrt 9    -- expect 3
#eval isqrt 16   -- expect 4
#eval isqrt 200  -- expect 14
#eval isqrt 0    -- expect 0

#eval! (powerIteration #[#[1, 0], #[0, 1]]).toInt  -- expect 65536

#eval! (computeSpectral #[#[1, 1], #[0, 1]]).matrix_size.toInt  -- expect 2
#eval! (computeSpectral #[#[1, 1], #[0, 1]]).rank.toInt          -- expect 2
#eval! (computeSpectral #[#[1, 1], #[0, 1]]).trace_val.toInt     -- expect 2

end SilverSight.PIST.Spectral
