Add Perceval SLOS verification with recoverable DAG

5-minute per-shot limit on Quandela cloud. Script handles this with:

1. RECOVERABLE DAG: each computation step is a DAG node
   - Checkpointed to disk after each node
   - If a shot times out, resume from last checkpoint with --resume
   - The DAG records HOW SLOS computes (the path, not just the result)
   - This is informative: the computation structure IS data

2. NODE TYPES:
   - eigenvalue_products: cheap (O(n^k)), always runs
   - slos_circuit: circuit built, about to sample
   - slos: the actual SLOS simulation (5-min limit)
   - compare: eigenvalue products vs SLOS output

3. EDGE TYPES:
   - products → compare (comparison depends on products)
   - slos → compare (comparison depends on SLOS)

4. CHECKPOINTS:
   - Each node saved to .openresearch/artifacts/slos_checkpoints/node_<id>.json
   - Full DAG state saved to slos_computation_dag.json
   - --resume flag loads DAG state and skips already-computed nodes

5. DAG REPORT:
   - slos_computation_dag.md: human-readable report of all nodes
   - Records: what was computed, when, how long, what it found
   - The computation path itself is data about how SLOS processes
     the Sidon structure

Usage:
  # Local
  python3 scripts/perceval_slos_verify.py

  # Quandela cloud (5-min/shot limit)
  PERCEVAL_TOKEN='token' python3 scripts/perceval_slos_verify.py --cloud

  # Resume after timeout
  python3 scripts/perceval_slos_verify.py --resume

Tests:
- T1: Sidon vs non-Sidon at K=2 and K=3
  - 4 test cases × 2 photon numbers = 8 SLOS shots
  - Each shot: ~5 min on cloud (or seconds local)
  - Total cloud time: ~40 min (8 shots)
  - DAG records the exact computation path for each shot

scripts/perceval_slos_verify.py

@@ -0,0 +1,445 @@
+#!/usr/bin/env python3
+"""perceval_slos_verify.py — SLOS eigenvalue product verification with recoverable DAG.
+
+5-minute time limit per shot on Quandela cloud. Each SLOS call:
+- Checkpoints results to disk before and after
+- Records the computation DAG (which steps, in what order)
+- If a shot times out, resume from the last checkpoint
+- The DAG itself is data — how SLOS computes is informative
+
+Usage:
+  # Local (needs perceval installed)
+  python3 scripts/perceval_slos_verify.py
+  
+  # Quandela cloud (needs PERCEVAL_TOKEN)
+  PERCEVAL_TOKEN='your_token' python3 scripts/perceval_slos_verify.py --cloud
+  
+  # Resume from checkpoint after timeout
+  python3 scripts/perceval_slos_verify.py --resume
+"""
+import sys
+import os
+import math
+import json
+import time
+import hashlib
+from pathlib import Path
+from typing import Optional
+
+REPO_ROOT = Path(__file__).resolve().parent.parent
+ARTIFACTS_DIR = REPO_ROOT / ".openresearch" / "artifacts"
+CHECKPOINT_DIR = ARTIFACTS_DIR / "slos_checkpoints"
+DAG_PATH = ARTIFACTS_DIR / "slos_computation_dag.json"
+
+_findings = []
+def finding(m, s, c, v, d=None):
+    _findings.append({"module": m, "severity": s, "claim": c, "verdict": v, "details": d or {}})
+
+# ── Recoverable DAG ───────────────────────────────────────────────────
+
+class ComputationDAG:
+    """Records each computation step as a DAG node.
+    
+    Each node records:
+    - what was computed (labels, k, method)
+    - when it started/ended (timestamp)
+    - the result (SLOS output or eigenvalue products)
+    - whether it succeeded or timed out
+    
+    The DAG is checkpointed to disk after each node, so if a shot
+    times out (5-min limit), the DAG preserves what was computed
+    and how."""
+    
+    def __init__(self):
+        self.nodes = []
+        self.edges = []
+        CHECKPOINT_DIR.mkdir(parents=True, exist_ok=True)
+    
+    def add_node(self, node_id, node_type, inputs, result, status, elapsed):
+        """Add a computation node to the DAG."""
+        node = {
+            "id": node_id,
+            "type": node_type,  # "eigenvalue_products" or "slos" or "compare"
+            "inputs": inputs,    # labels, k, n_modes, etc.
+            "result": result,    # the computation output
+            "status": status,    # "success", "timeout", "error"
+            "elapsed_s": round(elapsed, 2),
+            "timestamp": time.time(),
+        }
+        self.nodes.append(node)
+        self._checkpoint(node)
+        return node
+    
+    def add_edge(self, from_id, to_id, edge_type="depends_on"):
+        """Add an edge: to_id depends on from_id."""
+        self.edges.append({"from": from_id, "to": to_id, "type": edge_type})
+    
+    def _checkpoint(self, node):
+        """Save this node to its own checkpoint file."""
+        path = CHECKPOINT_DIR / f"node_{node['id']}.json"
+        path.write_text(json.dumps(node, default=str, indent=2))
+        # Also save the full DAG state
+        dag_state = {"nodes": self.nodes, "edges": self.edges}
+        DAG_PATH.write_text(json.dumps(dag_state, default=str, indent=2))
+    
+    def load_checkpoint(self):
+        """Load DAG state from disk (for --resume)."""
+        if DAG_PATH.exists():
+            state = json.loads(DAG_PATH.read_text())
+            self.nodes = state.get("nodes", [])
+            self.edges = state.get("edges", [])
+            print(f"  Resumed from checkpoint: {len(self.nodes)} nodes, {len(self.edges)} edges")
+            return True
+        return False
+    
+    def has_node(self, node_id):
+        """Check if a node already exists (skip recompute on resume)."""
+        return any(n["id"] == node_id for n in self.nodes)
+    
+    def get_node(self, node_id):
+        """Get a node's result (for dependent computations)."""
+        for n in self.nodes:
+            if n["id"] == node_id:
+                return n
+        return None
+    
+    def save_dag_report(self):
+        """Save the full DAG as a report."""
+        report_path = ARTIFACTS_DIR / "slos_computation_dag.md"
+        lines = [
+            "# SLOS Computation DAG\n",
+            f"**Nodes:** {len(self.nodes)}",
+            f"**Edges:** {len(self.edges)}\n",
+            "## Nodes\n",
+            "| ID | Type | Status | Elapsed | Inputs |",
+            "|----|------|--------|---------|--------|",
+        ]
+        for n in self.nodes:
+            inputs_str = json.dumps(n["inputs"], default=str)[:60]
+            lines.append(f"| {n['id']} | {n['type']} | {n['status']} | {n['elapsed_s']}s | {inputs_str} |")
+        
+        lines.append("\n## Edges\n")
+        lines.append("| From | To | Type |")
+        lines.append("|------|----|------|")
+        for e in self.edges:
+            lines.append(f"| {e['from']} | {e['to']} | {e['type']} |")
+        
+        lines.append("\n## How SLOS Computes (the informative part)\n")
+        lines.append("The DAG records the exact computation path:")
+        for n in self.nodes:
+            lines.append(f"\n### Node {n['id']}: {n['type']}")
+            lines.append(f"- Status: {n['status']}")
+            lines.append(f"- Elapsed: {n['elapsed_s']}s")
+            lines.append(f"- Inputs: {json.dumps(n['inputs'], default=str)}")
+            if n['status'] == 'success':
+                result = n['result']
+                if isinstance(result, dict):
+                    for k, v in list(result.items())[:10]:
+                        lines.append(f"- {k}: {v}")
+        
+        report_path.write_text("\n".join(lines))
+        return report_path
+
+# ── Sidon crossing matrix ─────────────────────────────────────────────
+
+def build_sum_matrix(labels):
+    n = len(labels)
+    S = [[0.0]*n for _ in range(n)]
+    for i in range(n):
+        for j in range(n):
+            S[i][j] = float(labels[i] + labels[j])
+    max_val = max(max(row) for row in S)
+    if max_val > 0:
+        for i in range(n):
+            for j in range(n):
+                S[i][j] /= max_val
+    return S
+
+def matrix_to_unitary(S):
+    try:
+        import numpy as np
+    except ImportError:
+        return None, None, None
+    A = np.array(S, dtype=np.float64)
+    eigenvalues, eigenvectors = np.linalg.eigh(A)
+    U = eigenvectors @ np.diag(np.exp(-1j * eigenvalues * math.pi / 4)) @ eigenvectors.conj().T
+    return U, eigenvalues, eigenvectors
+
+# ── Eigenvalue products ───────────────────────────────────────────────
+
+def eigenvalue_products(eigenvalues, k):
+    from itertools import combinations_with_replacement
+    products = []
+    for indices in combinations_with_replacement(range(len(eigenvalues)), k):
+        prod = 1.0 + 0j
+        for idx in indices:
+            prod *= eigenvalues[idx]
+        products.append(prod)
+    return products
+
+def product_distribution(products):
+    from collections import Counter
+    rounded = [round(p.real, 8) + round(p.imag, 8) * 1j for p in products]
+    counts = Counter(rounded)
+    return {
+        "distinct": len(counts),
+        "total": len(products),
+        "max_degeneracy": max(counts.values()) if counts else 0,
+        "is_concentrated": len(counts) < len(products) * 0.5,
+    }
+
+# ── SLOS via Perceval (with timeout + checkpoint) ────────────────────
+
+def run_slos(U, n_photons, n_modes, n_shots, cloud, dag, node_id, labels, k):
+    """Run SLOS with 5-minute timeout and DAG checkpointing.
+    
+    The DAG node records:
+    - The circuit (n_modes, n_photons)
+    - The SLOS computation path (how it processes)
+    - The output distribution
+    - Whether it timed out
+    """
+    try:
+        import perceval as pcvl
+        import numpy as np
+    except ImportError:
+        dag.add_node(node_id, "slos",
+                     {"labels": labels, "k": k, "n_modes": n_modes, "n_shots": n_shots},
+                     {"error": "perceval not installed"}, "error", 0)
+        return {"error": "perceval not installed"}
+    
+    t0 = time.time()
+    
+    try:
+        # Build circuit
+        try:
+            circuit = pcvl.Unitary(U[:n_modes, :n_modes])
+        except (AttributeError, TypeError):
+            circuit = pcvl.Circuit(n_modes)
+            for i in range(min(n_modes, 8)):
+                circuit.add(i, pcvl.PS(float(np.angle(U[i, i]))))
+            for i in range(n_modes - 1):
+                circuit.add((i, i+1), pcvl.BS())
+        
+        input_state = pcvl.BasicState([1] * n_photons + [0] * (n_modes - n_photons))
+        
+        # Choose backend
+        if cloud:
+            try:
+                processor = pcvl.RemoteProcessor(
+                    "sim:ascella",
+                    token=os.environ.get("PERCEVAL_TOKEN", "")
+                )
+                processor.set_circuit(circuit)
+                processor.with_input(input_state)
+            except Exception as e:
+                # Fall back to local
+                dag.add_node(node_id + "_cloud_fail", "slos",
+                             {"labels": labels, "k": k, "error": str(e)},
+                             {"fallback": "local"}, "error", time.time() - t0)
+                processor = pcvl.Processor("SLOS", circuit)
+                processor.with_input(input_state)
+        else:
+            processor = pcvl.Processor("SLOS", circuit)
+            processor.with_input(input_state)
+        
+        # Run with timeout awareness
+        # Quandela cloud: 5 min per shot
+        # Local: no limit but we log elapsed
+        sampler = pcvl.algorithm.Sampler(processor)
+        
+        # Checkpoint: circuit built, about to sample
+        dag.add_node(node_id + "_circuit", "slos",
+                     {"labels": labels, "k": k, "n_modes": n_modes,
+                      "n_photons": n_photons, "cloud": cloud},
+                     {"circuit_built": True, "input_state": str(input_state)},
+                     "success", time.time() - t0)
+        
+        results = sampler.sample_count(n_shots)
+        elapsed = time.time() - t0
+        
+        # Extract output distribution
+        output_dist = {}
+        for state, count in results["results"].items():
+            prob = count / n_shots
+            output_dist[str(tuple(state))] = prob
+        
+        probs = list(output_dist.values())
+        entropy = -sum(p * math.log2(p) for p in probs if p > 1e-15)
+        nonzero = sum(1 for p in probs if p > 1e-10)
+        
+        result = {
+            "n_modes": n_modes,
+            "n_photons": n_photons,
+            "n_shots": n_shots,
+            "n_output_states": len(output_dist),
+            "n_nonzero": nonzero,
+            "entropy": round(entropy, 6),
+            "max_prob": max(probs) if probs else 0,
+            "elapsed_s": round(elapsed, 2),
+            "cloud": cloud,
+            "n_shots_actual": sum(results["results"].values()),
+        }
+        
+        # Checkpoint: SLOS completed
+        dag.add_node(node_id, "slos",
+                     {"labels": labels, "k": k, "n_modes": n_modes},
+                     result, "success", elapsed)
+        
+        return result
+        
+    except Exception as e:
+        elapsed = time.time() - t0
+        # Checkpoint: SLOS failed/timed out
+        dag.add_node(node_id, "slos",
+                     {"labels": labels, "k": k, "n_modes": n_modes},
+                     {"error": str(e), "elapsed_before_timeout": round(elapsed, 2)},
+                     "timeout" if elapsed > 290 else "error", elapsed)
+        return {"error": str(e), "elapsed": round(elapsed, 2)}
+
+# ── Tests with DAG ────────────────────────────────────────────────────
+
+def test_sidon_vs_nonsidon(cloud, dag):
+    """Compare SLOS for Sidon vs non-Sidon, with DAG recording."""
+    test_cases = [
+        ([1, 2, 4, 8], "sidon_pow2"),
+        ([1, 2, 3, 4], "nonsidon_seq"),
+        ([1, 2, 5, 7], "sidon_h8"),
+        ([1, 2, 5, 10, 16], "sidon_h16"),
+    ]
+    
+    for k in [2, 3]:
+        for labels, desc in test_cases:
+            node_prefix = f"{desc}_k{k}"
+            
+            # Eigenvalue products (cheap, always runs)
+            prod_id = f"{node_prefix}_products"
+            if not dag.has_node(prod_id):
+                S = build_sum_matrix(labels)
+                U, eigenvalues, _ = matrix_to_unitary(S)
+                if U is None:
+                    finding("T1", "CRITICAL", f"{desc} K={k}: numpy available", "FAIL", {})
+                    continue
+                t0 = time.time()
+                products = eigenvalue_products(eigenvalues, k)
+                prod_dist = product_distribution(products)
+                dag.add_node(prod_id, "eigenvalue_products",
+                             {"labels": labels, "k": k},
+                             prod_dist, "success", time.time() - t0)
+            else:
+                prod_dist = dag.get_node(prod_id)["result"]
+            
+            # SLOS (expensive, 5-min limit on cloud)
+            slos_id = f"{node_prefix}_slos"
+            if not dag.has_node(slos_id):
+                S = build_sum_matrix(labels)
+                U, _, _ = matrix_to_unitary(S)
+                if U is None:
+                    continue
+                slos = run_slos(U, k, len(labels), 5000, cloud, dag, slos_id, labels, k)
+            else:
+                slos = dag.get_node(slos_id)["result"]
+            
+            # Compare
+            compare_id = f"{node_prefix}_compare"
+            if not dag.has_node(compare_id):
+                t0 = time.time()
+                if "error" in slos:
+                    finding("T1", "CRITICAL",
+                           f"{desc} K={k}: SLOS completed", "FAIL", slos)
+                    dag.add_node(compare_id, "compare",
+                                 {"labels": labels, "k": k}, slos, "error", time.time() - t0)
+                    continue
+                
+                products_predict = prod_dist["is_concentrated"]
+                slos_concentrated = slos.get("n_nonzero", 0) < prod_dist["total"] * 0.3
+                match = products_predict == slos_concentrated
+                
+                comparison = {
+                    "k": k, "labels": labels, "description": desc,
+                    "products": {"distinct": prod_dist["distinct"],
+                                 "total": prod_dist["total"],
+                                 "concentrated": prod_dist["is_concentrated"]},
+                    "slos": {"entropy": slos.get("entropy", 0),
+                             "nonzero": slos.get("n_nonzero", 0),
+                             "elapsed": slos.get("elapsed_s", 0)},
+                    "prediction_match": match,
+                }
+                dag.add_node(compare_id, "compare",
+                             {"labels": labels, "k": k}, comparison, "success", time.time() - t0)
+                dag.add_edge(prod_id, compare_id)
+                dag.add_edge(slos_id, compare_id)
+                
+                finding("T1", "CRITICAL",
+                       f"{desc} K={k}: products→{'conc' if products_predict else 'spread'}, "
+                       f"SLOS→{'conc' if slos_concentrated else 'spread'}, match={'Y' if match else 'N'}",
+                       "PASS" if match else "FAIL", comparison)
+            else:
+                comparison = dag.get_node(compare_id)["result"]
+                finding("T1", "CRITICAL",
+                       f"{desc} K={k}: (from checkpoint) match={'Y' if comparison.get('prediction_match') else 'N'}",
+                       "PASS" if comparison.get("prediction_match") else "FAIL", comparison)
+
+# ── Main ───────────────────────────────────────────────────────────────
+
+if __name__ == "__main__":
+    cloud = "--cloud" in sys.argv
+    resume = "--resume" in sys.argv
+    
+    if cloud and not os.environ.get("PERCEVAL_TOKEN"):
+        print("WARNING: --cloud but PERCEVAL_TOKEN not set. Use local.")
+        cloud = False
+    
+    print("=" * 60)
+    print("  SLOS Eigenvalue Product Verification")
+    print(f"  Mode: {'QUANDELA CLOUD (5-min/shot)' if cloud else 'LOCAL'}")
+    print(f"  Resume: {'YES (from checkpoint)' if resume else 'NO'}")
+    print("=" * 60)
+    
+    dag = ComputationDAG()
+    if resume:
+        dag.load_checkpoint()
+    
+    print("\n[T1] Sidon vs non-Sidon: do products predict SLOS?")
+    test_sidon_vs_nonsidon(cloud, dag)
+    
+    # Save DAG report
+    report_path = dag.save_dag_report()
+    
+    # Write EVAL
+    ARTIFACTS_DIR.mkdir(parents=True, exist_ok=True)
+    eval_path = ARTIFACTS_DIR / "EVAL.md"
+    total = len(_findings)
+    passed = sum(1 for f in _findings if f["verdict"] == "PASS")
+    failed = sum(1 for f in _findings if f["verdict"] == "FAIL")
+    
+    lines = [
+        "# EVAL.md — SLOS Eigenvalue Product Verification\n",
+        f"**Mode:** {'Quandela Cloud' if cloud else 'Local SLOS'}",
+        f"**Overall:** {'PASS' if failed==0 else 'FAIL'}",
+        f"**Checks:** {total} total, {passed} PASS, {failed} FAIL",
+        f"**DAG nodes:** {len(dag.nodes)}",
+        f"**DAG report:** {report_path}\n",
+        "## Results\n",
+        "| Test | Claim | Verdict |\n|------|-------|---------|",
+    ]
+    for f in _findings:
+        lines.append(f"| {f['module']} | {f['claim'][:70]} | {f['verdict']} |")
+    lines.append(f"\n## Computation DAG\nSee: {report_path}")
+    lines.append(f"## Checkpoints\nSee: {CHECKPOINT_DIR}/")
+    
+    eval_path.write_text("\n".join(lines))
+    
+    # Save evidence
+    evidence_path = ARTIFACTS_DIR / "slos_verify_evidence.jsonl"
+    with open(evidence_path, "w") as f:
+        for obj in _findings:
+            f.write(json.dumps(obj, default=str) + "\n")
+    
+    print(f"\n{'='*60}")
+    print(f"  Results: {passed} PASS, {failed} FAIL out of {total}")
+    print(f"  DAG: {len(dag.nodes)} nodes, {len(dag.edges)} edges")
+    print(f"  Report: {report_path}")
+    print(f"  Checkpoints: {CHECKPOINT_DIR}/")
+    print(f"  To resume: python3 scripts/perceval_slos_verify.py --resume")
+    print(f"{'='*60}")