#!/usr/bin/env python3
"""
pipeline_core.py — Module-swappable six-stage search engine.

Uses the ACTUAL SilverSight chiral implementation:
  - ChiralLabel: 4 types (achiral_stable, chiral_scarred, left_handed, right_handed)
  - Phase: Z/360Z at 45° steps (0,45,90,135,180,225,270,315)
  - Rossby drift: weights left=+65536, right=-65536, scarred=+32768, achiral=0
  - Quaternion basis: achiral=1, left=i, right=j, scarred=k (HopfFibration)
  - Golden angle: 25042 Q16_16 (2π/φ²), winding mod 28
  - Fisher-Rao sphere: p → 2√p maps simplex to S²

No floats. No native_decide. All Q16_16 integer arithmetic.

Usage:
  from pipeline_core import Pipeline, BraidStorm, TreeBraid, AngrySphinx,
                           MultisurfacePacker, COUCHFilter, SidonFilter
  pipe = Pipeline([...])
  result = pipe.run(labels, S, moduli)
"""

import sys, math, json, time, hashlib, random
from abc import ABC, abstractmethod
from dataclasses import dataclass, field
from pathlib import Path
from itertools import product
from collections import Counter
from typing import Any

REPO_ROOT = Path(__file__).resolve().parent.parent
ARTIFACTS_DIR = REPO_ROOT / ".openresearch" / "artifacts"
ARTIFACTS_DIR.mkdir(parents=True, exist_ok=True)

# ── Q16_16 constants (no floats) ──────────────────────────────────────
Q16_ONE = 65536
Q16_HALF = 32768          # 0.5 in Q16_16 (chiral_scarred weight)
Q16_THRESHOLD_COUCH = 49152  # 0.75
Q16_SUBLEQ_SELFLOOP = 53908  # 0.823
Q16_AVX_SELFLOOP = 57942     # 0.885
Q16_RING_SELFLOOP = 0        # 0.0
Q16_GOLDEN_ANGLE = 25042     # 2π/φ² in Q16_16
EXOTIC_CLASSES = 28          # Durán exotic sphere classes
Q16_TWO_PI = 411775          # 2π in Q16_16 (≈ 6.28318)

# ── ChiralLabel (from BraidStateN.lean) ───────────────────────────────
CHIRAL_LABELS = ["achiral_stable", "chiral_scarred",
                 "left_handed_mass_bias", "right_handed_vector_bias"]

# Rossby drift weights (from rossbyDriftFromChirality)
ROSSBY_WEIGHTS = {
    "achiral_stable": 0,           # Q16_16.zero
    "chiral_scarred": Q16_HALF,    # 32768 = 0.5
    "left_handed_mass_bias": Q16_ONE,   # +1
    "right_handed_vector_bias": -Q16_ONE,  # -1
}

# Phase angles (from HachimojiBase.lean) — Z/360Z at 45° steps
PHASE_ANGLES = [0, 45, 90, 135, 180, 225, 270, 315]

# Chirality from phase (from HachimojiBase.chirality)
def chirality_from_phase(phase):
    """Map phase angle to chirality (HachimojiBase.lean line 176)."""
    if phase in (0, 90, 180): return "ambidextrous"
    if phase in (45, 135): return "left"
    if phase in (225, 270, 315): return "right"
    return "right"

# Quaternion basis mapping (from HopfFibration.lean ofChiralLabel)
# achiral_stable → 1 (scalar), left → i, right → j, scarred → k
QUATERNION_BASIS = {
    "achiral_stable": (Q16_ONE, 0, 0, 0),       # 1
    "left_handed_mass_bias": (0, Q16_ONE, 0, 0),  # i
    "right_handed_vector_bias": (0, 0, Q16_ONE, 0),  # j
    "chiral_scarred": (0, 0, 0, Q16_ONE),       # k
}

# ChiralLabel from chirality (mapping Hachimoji chirality → BraidStateN ChiralLabel)
def chiral_label_from_phase(phase):
    """Map phase to ChiralLabel (combining HachimojiBase + BraidStateN)."""
    chi = chirality_from_phase(phase)
    if chi == "ambidextrous": return "achiral_stable"
    if chi == "left": return "left_handed_mass_bias"
    if chi == "right": return "right_handed_vector_bias"
    return "achiral_stable"

# 8-strand default chiral assignments (from BraidStateN.lean rossbyLabels8)
# Cycles through phases 0..315 for strands 0..7
DEFAULT_PHASES_8 = PHASE_ANGLES  # [0, 45, 90, 135, 180, 225, 270, 315]


def rossby_drift(chiral_labels):
    """Compute Rossby drift from per-strand chiral labels.
    (Python port of rossbyDriftFromChirality from BraidStateN.lean)
    Returns (asymmetry_q16, is_active)."""
    total = 0
    for label in chiral_labels:
        total += ROSSBY_WEIGHTS.get(label, 0)
    return total, total != 0

def helical_residue(k):
    """Helical boundary residue at step k: ⌊k·ψ⌋ mod 28.
    (Python port of helicalResidue from HopfFibration.lean)"""
    return ((k * Q16_GOLDEN_ANGLE) // Q16_ONE) % EXOTIC_CLASSES

def quaternion_multiply(q1, q2):
    """Multiply two Q16_16 quaternions (a,b,c,d) where q = a + bi + cj + dk.
    All integer arithmetic, result scaled by Q16_ONE (divide at end)."""
    a1, b1, c1, d1 = q1
    a2, b2, c2, d2 = q2
    # Hamilton product: (a1a2 - b1b2 - c1c2 - d1d2) + ...i + ...j + ...k
    a = (a1*a2 - b1*b2 - c1*c2 - d1*d2) // Q16_ONE
    b = (a1*b2 + b1*a2 + c1*d2 - d1*c2) // Q16_ONE
    c = (a1*c2 - b1*d2 + c1*a2 + d1*b2) // Q16_ONE
    d = (a1*d2 + b1*c2 - c1*b2 + d1*a2) // Q16_ONE
    return (a, b, c, d)


# ── Config: the unit that flows through the pipeline ──────────────────

@dataclass
class Config:
    """A single configuration flowing through the pipeline."""
    chiral: tuple             # per-strand ChiralLabel string tuple
    phases: tuple             # per-strand phase angles (0,45,...,315)
    labels: tuple             # Sidon label integers
    S: int                    # reflection point
    moduli: tuple             # CRT moduli
    cost: int = 0
    self_loop: int = 0        # Q16_16 raw
    sidon_score: int = 0      # Q16_16 raw
    collisions: int = 0
    rossby_drift: int = 0     # Q16_16 raw
    helical_residue: int = 0  # 0..27
    metadata: dict = field(default_factory=dict)


# ── Pipeline Context ──────────────────────────────────────────────────

@dataclass
class PipelineContext:
    crossing_pairs: tuple = ()
    groups: tuple = ()
    seed: int = 0
    stage_timings: dict = field(default_factory=dict)


# ── Filter interface ──────────────────────────────────────────────────

class Filter(ABC):
    @abstractmethod
    def apply(self, configs: list[Config], ctx: PipelineContext) -> list[Config]: ...
    @property
    @abstractmethod
    def name(self) -> str: ...
    def run_stage(self, configs, ctx):
        t0 = time.time()
        result = self.apply(configs, ctx)
        ctx.stage_timings[self.name] = {
            "input": len(configs), "output": len(result),
            "time_s": round(time.time() - t0, 4)}
        return result


# ── Stage 1: BraidStorm — Generate chiral configurations ─────────────

class BraidStorm(Filter):
    """Generates chiral configurations by assigning ChiralLabels to strands.

    Each strand gets a phase angle from Z/360Z (45° steps). The chiral
    label is DERIVED from the phase (HachimojiBase.lean). Different phase
    assignments = different chiral configurations.

    For 8 strands with 8 phases: 8! = 40320 permutations (too many).
    Use k swap positions: 2^k configurations from positional swaps.
    """
    def __init__(self, k=8):
        self.k = k
    @property
    def name(self): return f"BraidStorm(k={self.k})"
    def apply(self, configs, ctx):
        if not configs: return []
        template = configs[0]
        # Generate 2^k chiral configs via positional swaps
        all_swaps = list(product([0, 1], repeat=self.k))
        result = []
        for swap_config in all_swaps:
            phases = list(template.phases)
            # Apply swaps: swap[j]=1 swaps phases[j] and phases[j+1]
            for j in range(min(self.k, len(phases)-1)):
                if swap_config[j] == 1:
                    phases[j], phases[j+1] = phases[j+1], phases[j]
            chiral = tuple(chiral_label_from_phase(p) for p in phases)
            result.append(Config(
                chiral=chiral, phases=tuple(phases),
                labels=template.labels, S=template.S, moduli=template.moduli))
        return result


# ── Stage 2: TreeBraid — Factorize via braid relations ────────────────

class TreeBraid(Filter):
    """Factorize crossing space. σ_i σ_j = σ_j σ_i when |i-j| >= 2."""
    @property
    def name(self): return "TreeBraid"
    def apply(self, configs, ctx):
        if not configs: return []
        k = self.k = len(configs[0].phases) - 1
        pairs = tuple((i, i+1) for i in range(k))
        ctx.crossing_pairs = pairs
        groups = self._factorize(k, pairs)
        ctx.groups = groups
        for c in configs:
            c.metadata["groups"] = groups
        return configs
    def _factorize(self, k, pairs):
        groups, remaining = [], list(range(k))
        while remaining:
            group = [remaining[0]]
            for idx in remaining[1:]:
                si, sj = pairs[idx]
                if all(abs(si - pairs[g][0]) >= 2 and abs(si - pairs[g][1]) >= 2 and
                       abs(sj - pairs[g][0]) >= 2 and abs(sj - pairs[g][1]) >= 2
                       for g in group):
                    group.append(idx)
            for g in group: remaining.remove(g)
            groups.append(tuple(group))
        return tuple(groups)


# ── Stage 3: AngrySphinx — Resource budget ───────────────────────────

class AngrySphinx(Filter):
    """Filter by compute budget. Cost = 2^(scarred+left+right count)."""
    def __init__(self, budget=128): self.budget = budget
    @property
    def name(self): return f"AngrySphinx(budget={self.budget})"
    def apply(self, configs, ctx):
        result = []
        for c in configs:
            active = sum(1 for cl in c.chiral if cl != "achiral_stable")
            c.cost = 1 << active
            if c.cost <= self.budget: result.append(c)
        return result


# ── Stage 4: MultisurfacePacker ───────────────────────────────────────

class MultisurfacePacker(Filter):
    def __init__(self, max_surfaces=64): self.max_surfaces = max_surfaces
    @property
    def name(self): return f"MultisurfacePacker(max={self.max_surfaces})"
    def apply(self, configs, ctx):
        if len(configs) <= self.max_surfaces: return configs
        return sorted(configs, key=lambda c: c.cost)[:self.max_surfaces]


# ── Stage 5: COUCH — Geometric stability via Rossby drift ─────────────

class COUCHFilter(Filter):
    """COUCH gate: two-stage stability check.

    Stage A — Rossby/Kelvin regime check (from BraidStateN.lean):
      - Rossby drift ≠ 0 → Rossby regime (dispersive, active mixing) → PASS
      - Rossby drift = 0 → Kelvin regime (boundary-trapped, NO mixing) → FAIL

      The Kelvin regime is the degenerate case: perfectly balanced chiral
      distribution (left+right cancel, no scarred). Energy dissipation rate
      is ZERO (proven: rossby_energy_dissipation_rate requires isActive).
      This is the "q=1 degenerate" / "rational surface" / "AVX-512 stuck"
      case — the system can't mix.

      rossbyLabels8 (BraidStateN.lean line 327): alternating left/right
      → drift ≠ 0 → Rossby → dispersive → COUCH passes
      kelvinLabels8 (BraidStateN.lean line 338): all achiral_stable
      → drift = 0 → Kelvin → boundary-trapped → COUCH fails

    Stage B — Scarred contention check (from HCMR):
      - scarred_count → self_loop proxy → threshold check
      - High scarred count = high contention = AVX-512 regime

    A config passes COUCH iff BOTH stages pass.
    """
    def __init__(self, threshold=Q16_THRESHOLD_COUCH): self.threshold = threshold
    @property
    def name(self): return f"COUCH(threshold={self.threshold})"
    def apply(self, configs, ctx):
        result = []
        for c in configs:
            # Stage A: Rossby/Kelvin regime
            drift, is_active = rossby_drift(c.chiral)
            c.rossby_drift = drift
            if not is_active:
                # Kelvin regime: boundary-trapped, no mixing → COUCH FAILS
                c.metadata["regime"] = "kelvin"
                c.self_loop = Q16_AVX_SELFLOOP  # stuck, like AVX-512
                continue
            c.metadata["regime"] = "rossby"

            # Stage B: scarred contention
            scarred_count = sum(1 for cl in c.chiral if cl == "chiral_scarred")
            c.self_loop = (Q16_AVX_SELFLOOP * scarred_count) // max(len(c.chiral), 1)
            if c.self_loop >= self.threshold:
                continue  # too much contention

            # Helical residue (winding mod 28, from HopfFibration.lean)
            step = c.metadata.get("step", 0)
            c.helical_residue = helical_residue(step)
            result.append(c)
        return result


# ── Stage 6: Sidon Filter — Algebraic uniqueness ──────────────────────

class SidonFilter(Filter):
    """Sidon filter using CRT reconstruction.

    The chiral configuration determines which phase (and thus which
    ChiralLabel) is at each strand position. Different permutations
    pair different labels with different moduli (position-dependent).

    This is NOT a ring automorphism — it's a positional permutation
    on the sphere (S² via Fisher-Rao embedding p → 2√p).
    """
    @property
    def name(self): return "SidonFilter"
    def apply(self, configs, ctx):
        result = []
        for c in configs:
            embedded = self._embed(c)
            collisions = self._sidon_check(embedded, c.moduli)
            c.collisions = collisions
            total = len(c.labels) * (len(c.labels) + 1) // 2
            c.sidon_score = Q16_ONE - (Q16_ONE * collisions) // max(total, 1)
            if collisions == 0: result.append(c)
        return result
    def _embed(self, c):
        """CRT embed with positional chirality.
        The phase at each position determines the ChiralLabel, which
        determines the quaternion basis (1,i,j,k) for the Hopf fibration.
        The CRT modulus at each position encodes the geometric constraint
        at that spherical position."""
        embedded = []
        for pos, (label, phase) in enumerate(zip(c.labels, c.phases)):
            chiral_label = chiral_label_from_phase(phase)
            quat = QUATERNION_BASIS[chiral_label]
            # CRT: identity = label % L0, reflection = (S-label) % Lj
            row = [label % c.moduli[0]]
            for j in range(1, len(c.moduli)):
                row.append((c.S - label) % c.moduli[j])
            # Store quaternion alongside CRT residues
            embedded.append({"residues": row, "quaternion": quat,
                             "chiral": chiral_label, "phase": phase})
        return embedded
    def _sidon_check(self, embedded, moduli):
        """Check Sidon on CRT sums (proven chiral-invariant for negation,
        but positional permutation of phases changes which label gets
        which modulus, which CAN break Sidon)."""
        M = 1
        for m in moduli: M *= m
        vals = [self._crt_reconstruct(e["residues"], moduli) for e in embedded]
        sums = []
        for i in range(len(vals)):
            for j in range(i, len(vals)):
                sums.append((vals[i] + vals[j]) % M)
        counts = Counter(sums)
        return sum(cnt - 1 for cnt in counts.values())
    def _crt_reconstruct(self, residues, moduli):
        M = 1
        for m in moduli: M *= m
        x = 0
        for r, m in zip(residues, moduli):
            Mi = M // m
            inv = self._modinv(Mi % m, m)
            if inv is None: return 0
            x = (x + r * Mi * inv) % M
        return x
    def _egcd(self, a, b):
        if b == 0: return a, 1, 0
        g, x, y = self._egcd(b, a % b)
        return g, y, x - (a // b) * y
    def _modinv(self, a, m):
        g, x, _ = self._egcd(a % m, m)
        return x % m if g == 1 else None


# ── Swappable: Quaternion Product Sidon Filter ────────────────────────

class QuaternionSidonFilter(SidonFilter):
    """Sidon filter using quaternion products (not CRT sums).

    Uses the ACTUAL quaternion basis mapping from HopfFibration.lean:
      achiral_stable → 1, left → i, right → j, scarred → k

    The quaternion product q_i · q_j is NOT invariant under positional
    permutation (because different positions have different chiral labels
    → different quaternion basis elements → different products).

    This is the filter that actually discriminates chiral configurations.
    """
    @property
    def name(self): return "QuaternionSidonFilter"
    def _sidon_check(self, embedded, moduli):
        """Check Sidon on quaternion products (Hamilton product).
        Two pairs (i,j) and (k,l) collide if q_i·q_j = q_k·q_l."""
        n = len(embedded)
        products = []
        for i in range(n):
            for j in range(i, n):
                qi = embedded[i]["quaternion"]
                qj = embedded[j]["quaternion"]
                prod = quaternion_multiply(qi, qj)
                products.append(prod)
        counts = Counter(products)
        return sum(cnt - 1 for cnt in counts.values())


# ── Pipeline ──────────────────────────────────────────────────────────

class Pipeline:
    def __init__(self, filters): self.filters = filters
    def run(self, labels, S, moduli, phases=None):
        t0 = time.time()
        ctx = PipelineContext(seed=hash((tuple(labels), S, tuple(moduli))) % (2**31))
        if phases is None:
            phases = tuple(DEFAULT_PHASES_8[:len(labels)])
        configs = [Config(
            chiral=tuple(chiral_label_from_phase(p) for p in phases),
            phases=phases, labels=tuple(labels), S=S, moduli=tuple(moduli))]
        for f in self.filters:
            configs = f.run_stage(configs, ctx)
        elapsed = time.time() - t0
        result = {
            "experiment": "pipeline_core",
            "timestamp": time.strftime("%Y-%m-%dT%H:%M:%SZ", time.gmtime()),
            "stages": [f.name for f in self.filters],
            "labels": list(labels), "S": S, "moduli": list(moduli),
            "phases": list(phases),
            "stage_timings": ctx.stage_timings,
            "total_output": len(configs),
            "elapsed_s": round(elapsed, 4),
            "final_configs": [{
                "phases": list(c.phases),
                "chiral": list(c.chiral),
                "rossby_drift": c.rossby_drift,
                "helical_residue": c.helical_residue,
                "self_loop": c.self_loop,
                "sidon_score": c.sidon_score,
                "collisions": c.collisions,
            } for c in configs[:20]],  # first 20 for brevity
        }
        content = json.dumps(result, indent=2, sort_keys=True, default=str)
        result["sha256"] = hashlib.sha256(content.encode()).hexdigest()
        print(f"\n{'='*60}")
        print(f"  PIPELINE: {' → '.join(f.name for f in self.filters)}")
        print(f"{'='*60}")
        for s, t in ctx.stage_timings.items():
            print(f"  {s:35s} {t['input']:6d} → {t['output']:6d}  ({t['time_s']:.4f}s)")
        print(f"{'='*60}")
        print(f"  Total output: {len(configs)}  ({elapsed:.2f}s)")
        print(f"{'='*60}")
        return result


if __name__ == "__main__":
    import argparse
    parser = argparse.ArgumentParser(description="Module-swappable pipeline")
    parser.add_argument("--strands", type=int, default=8)
    parser.add_argument("--budget", type=int, default=128)
    parser.add_argument("--surfaces", type=int, default=64)
    parser.add_argument("--filter", choices=["crt", "quat"], default="crt")
    parser.add_argument("--output", default="pipeline_result.json")
    args = parser.parse_args()

    labels = [1, 2, 4, 8, 16, 32, 64, 128]
    S = 128
    moduli = [7, 3, 5, 11, 13, 17, 19, 23, 29]
    phases = tuple(DEFAULT_PHASES_8[:len(labels)])

    sidon = SidonFilter() if args.filter == "crt" else QuaternionSidonFilter()

    pipe = Pipeline([
        BraidStorm(k=args.strands),
        TreeBraid(),
        AngrySphinx(budget=args.budget),
        MultisurfacePacker(max_surfaces=args.surfaces),
        COUCHFilter(),
        sidon,
    ])
    result = pipe.run(labels=labels, S=S, moduli=moduli, phases=phases)
    (ARTIFACTS_DIR / "pipeline_result.json").write_text(
        json.dumps(result, indent=2, default=str))
    print(f"\nResults → {ARTIFACTS_DIR / 'pipeline_result.json'}")