From 7b0c0a363b10d40eda3ec58968d45f9dbbfa8e63 Mon Sep 17 00:00:00 2001 From: allaun Date: Sun, 21 Jun 2026 01:02:07 -0500 Subject: [PATCH] feat(infra): citation crawler and spectral cross-validation research shims Add exploratory Python shims for: - OpenAlex citation crawling - Semantic Scholar citation crawling - PageRank eigenvalue survey - SNAP PIST spectral cross-validation - Fiedler vector validation --- .../shim/openalex_citation_crawler.py | 310 ++++++++++++ .../shim/pagerank_eigenvalue_survey.py | 229 +++++++++ 4-Infrastructure/shim/s2_citation_crawler.py | 333 +++++++++++++ .../shim/snap_pist_spectral_crossval.py | 444 ++++++++++++++++++ 4-Infrastructure/shim/validate_fiedler.py | 438 +++++++++++++++++ 5 files changed, 1754 insertions(+) create mode 100644 4-Infrastructure/shim/openalex_citation_crawler.py create mode 100644 4-Infrastructure/shim/pagerank_eigenvalue_survey.py create mode 100644 4-Infrastructure/shim/s2_citation_crawler.py create mode 100644 4-Infrastructure/shim/snap_pist_spectral_crossval.py create mode 100644 4-Infrastructure/shim/validate_fiedler.py diff --git a/4-Infrastructure/shim/openalex_citation_crawler.py b/4-Infrastructure/shim/openalex_citation_crawler.py new file mode 100644 index 00000000..ad3d3eae --- /dev/null +++ b/4-Infrastructure/shim/openalex_citation_crawler.py @@ -0,0 +1,310 @@ +#!/usr/bin/env python3 +""" +openalex_citation_crawler.py + +BFS citation graph crawler using the OpenAlex API — completely free, no key, +50 req/sec in the polite pool (just pass your email as mailto= param). + +Also provides a Barabási-Albert synthetic graph for instant offline testing +when no network is available (or when you just want a known structure). + +Feeds into pagerank_eigenvalue_survey.run_survey() for spectral bin analysis. + +Usage: + # Real citation graph (online): + python openalex_citation_crawler.py --seed "PageRank Brin Page 1998" --hops 2 + + # Synthetic scale-free graph (offline, instant): + python openalex_citation_crawler.py --synthetic --nodes 300 + + # Save and survey: + python openalex_citation_crawler.py --seed "HITS Kleinberg 1999" --hops 2 --out hits_graph.json + +Requires: requests, networkx, scipy, numpy + uv pip install requests networkx scipy numpy +""" + +import argparse +import json +import sys +import time +from collections import deque +from dataclasses import dataclass, field, asdict +from pathlib import Path +from typing import List, Dict, Optional, Set + +import requests +import networkx as nx +import numpy as np + +sys.path.insert(0, str(Path(__file__).parent)) +try: + from pagerank_eigenvalue_survey import run_survey + HAS_SURVEY = True +except ImportError: + HAS_SURVEY = False + print("[warn] pagerank_eigenvalue_survey not found — will print basic PageRank only") + +# ── OpenAlex API ────────────────────────────────────────────────────────────── + +OA_BASE = "https://api.openalex.org" +# Using the project email puts requests into the polite pool (50 req/sec). +OA_EMAIL = "allaunjsilverfox@gmail.com" + +@dataclass +class OAPaper: + work_id: str + title: str + year: Optional[int] + cite_count: int + doi: Optional[str] + references: List[str] = field(default_factory=list) # OpenAlex work IDs + citations: List[str] = field(default_factory=list) # OpenAlex work IDs + +class OAClient: + def __init__(self, rate_delay: float = 0.05): # 20 req/sec, well under 50 + self._last = 0.0 + self._delay = rate_delay + self._stats = {"requests": 0, "errors": 0} + + def _get(self, path: str, params: Dict = None) -> Optional[dict]: + elapsed = time.monotonic() - self._last + if elapsed < self._delay: + time.sleep(self._delay - elapsed) + self._last = time.monotonic() + self._stats["requests"] += 1 + + p = {"mailto": OA_EMAIL, **(params or {})} + try: + resp = requests.get(f"{OA_BASE}{path}", params=p, timeout=15) + if resp.status_code == 429: + print(" [rate-limit] waiting 30s") + time.sleep(30) + return self._get(path, params) + resp.raise_for_status() + return resp.json() + except Exception as exc: + self._stats["errors"] += 1 + print(f" [oa error] {path}: {exc}") + return None + + def search(self, query: str, limit: int = 5) -> List[dict]: + data = self._get("/works", { + "search": query, + "per-page": limit, + "select": "id,title,publication_year,cited_by_count,doi", + }) + return data.get("results", []) if data else [] + + def paper(self, work_id: str) -> Optional[OAPaper]: + # work_id may be full URL or short form like W2741809807 + short_id = work_id.replace("https://openalex.org/", "") + data = self._get(f"/works/{short_id}", { + "select": "id,title,publication_year,cited_by_count,doi,referenced_works" + }) + if not data or "id" not in data: + return None + return OAPaper( + work_id=data["id"], + title=data.get("title") or "", + year=data.get("publication_year"), + cite_count=data.get("cited_by_count", 0), + doi=data.get("doi"), + references=[r for r in data.get("referenced_works", []) if r], + ) + + def citing(self, work_id: str, limit: int = 50) -> List[str]: + """IDs of papers that cite work_id.""" + short_id = work_id.replace("https://openalex.org/", "") + data = self._get("/works", { + "filter": f"cites:{short_id}", + "per-page": limit, + "select": "id", + "sort": "cited_by_count:desc", # take the most-cited citers first + }) + if not data: + return [] + return [r["id"] for r in data.get("results", []) if r.get("id")] + + def stats(self) -> str: + return f"requests={self._stats['requests']} errors={self._stats['errors']}" + + +# ── BFS ─────────────────────────────────────────────────────────────────────── + +def crawl_openalex(seed_id: str, hops: int, max_refs: int, max_cites: int, + max_papers: int) -> Dict[str, OAPaper]: + client = OAClient() + visited: Set[str] = set() + papers: Dict[str, OAPaper] = {} + queue: deque = deque([(seed_id, 0)]) + + print(f"\nOpenAlex BFS: seed={seed_id} hops={hops} max={max_papers}") + + while queue and len(papers) < max_papers: + wid, depth = queue.popleft() + if wid in visited: + continue + visited.add(wid) + + paper = client.paper(wid) + if not paper: + continue + papers[wid] = paper + print(f" [{depth}/{hops}] {len(papers):>4} cite={paper.cite_count:>6}" + f" {paper.title[:65]}") + + if depth < hops: + # References are already in the paper response + ref_ids = paper.references[:max_refs] + cite_ids = client.citing(wid, limit=max_cites) + paper.citations = cite_ids + + for nid in ref_ids + cite_ids: + if nid and nid not in visited: + queue.append((nid, depth + 1)) + + print(f"\nDone: {len(papers)} papers | {client.stats()}") + return papers + + +# ── Synthetic scale-free graph ──────────────────────────────────────────────── + +def build_synthetic_webgraph(n: int = 300, m: int = 3, spam_nodes: int = 15, + seed: int = 42) -> nx.DiGraph: + """ + Barabási-Albert preferential attachment graph → directed citation graph. + + Properties matching the real web: + - Power-law degree distribution (few hubs, many leaves) + - Small-world connectivity + - Natural community structure emerges at n≥200 + + Spam injection: `spam_nodes` low-degree nodes → hub node 0 (artificial + authority boost, should trigger bad-link gate on merger). + """ + rng = np.random.default_rng(seed) + ug = nx.barabasi_albert_graph(n, m, seed=seed) + g = ug.to_directed() + + # Prune ~30% of back-edges to make it more citation-like (mostly one-way) + edges_to_remove = [(u, v) for u, v in list(g.edges()) + if u > v and rng.random() < 0.7] + g.remove_edges_from(edges_to_remove) + + # Add spam link farm → node 0 (highest PageRank target) + for i in range(n, n + spam_nodes): + g.add_node(i, spam=True) + g.add_edge(i, 0) + if i > n: + g.add_edge(i, i - 1) # farm cross-links + + print(f"Synthetic graph: {g.number_of_nodes()} nodes, " + f"{g.number_of_edges()} edges (BA n={n} m={m} + {spam_nodes} spam)") + return g + + +# ── Build NetworkX graph from OA papers ────────────────────────────────────── + +def oa_to_networkx(papers: Dict[str, OAPaper]) -> nx.DiGraph: + g = nx.DiGraph() + id_to_int = {wid: i for i, wid in enumerate(papers)} + + for wid, p in papers.items(): + nid = id_to_int[wid] + g.add_node(nid, work_id=wid, title=p.title[:80], + year=p.year or 0, cite_count=p.cite_count) + + for wid, p in papers.items(): + src = id_to_int[wid] + for ref_id in p.references: + if ref_id in id_to_int: + g.add_edge(src, id_to_int[ref_id]) # src cites ref + + print(f"Graph: {g.number_of_nodes()} nodes, {g.number_of_edges()} edges") + return g + + +# ── Lightweight survey fallback ─────────────────────────────────────────────── + +def basic_survey(g: nx.DiGraph): + """Minimal survey when pagerank_eigenvalue_survey is not importable.""" + pr = nx.pagerank(g, alpha=0.85) + try: + hubs, auths = nx.hits(g, max_iter=300) + except Exception: + hubs = auths = {n: 1 / len(g) for n in g} + + top = sorted(pr.items(), key=lambda x: -x[1])[:15] + print("\n[PageRank top 15]") + for nid, score in top: + title = g.nodes[nid].get("title", f"node-{nid}") + year = g.nodes[nid].get("year", "") + print(f" {score:.5f} [{year}] {title}") + + top_auth = sorted(auths.items(), key=lambda x: -x[1])[:10] + print("\n[HITS authority top 10]") + for nid, score in top_auth: + title = g.nodes[nid].get("title", f"node-{nid}") + print(f" {score:.5f} {title}") + + +# ── Main ────────────────────────────────────────────────────────────────────── + +def main(): + ap = argparse.ArgumentParser( + description="OpenAlex BFS citation crawler + spectral survey (no API key needed)") + mode = ap.add_mutually_exclusive_group() + mode.add_argument("--seed", default="PageRank Brin Page 1998", + help="Title/author query (default: PageRank paper)") + mode.add_argument("--seed-id", help="OpenAlex work ID e.g. W2741809807") + mode.add_argument("--synthetic", action="store_true", + help="Use synthetic Barabási-Albert graph (offline, instant)") + + ap.add_argument("--nodes", type=int, default=300, help="Synthetic graph node count") + ap.add_argument("--hops", type=int, default=2) + ap.add_argument("--max-refs", type=int, default=40) + ap.add_argument("--max-cites", type=int, default=30) + ap.add_argument("--max-papers", type=int, default=300) + ap.add_argument("--out", default=None, help="Save paper JSON here") + ap.add_argument("--no-survey", action="store_true") + args = ap.parse_args() + + if args.synthetic: + g = build_synthetic_webgraph(n=args.nodes) + else: + client = OAClient() + if args.seed_id: + seed_id = args.seed_id + if not seed_id.startswith("https://"): + seed_id = f"https://openalex.org/{seed_id}" + else: + print(f"Searching OpenAlex: {args.seed}") + results = client.search(args.seed, limit=3) + if not results: + print("ERROR: no results") + sys.exit(1) + best = results[0] + seed_id = best["id"] + print(f" seed: [{best.get('publication_year','')}] {best.get('title','?')}") + print(f" id: {seed_id}") + + papers = crawl_openalex( + seed_id, args.hops, args.max_refs, args.max_cites, args.max_papers) + + if args.out: + Path(args.out).write_text( + json.dumps({wid: asdict(p) for wid, p in papers.items()}, indent=2)) + print(f"Saved → {args.out}") + + g = oa_to_networkx(papers) + + if not args.no_survey: + if HAS_SURVEY: + run_survey(g) + else: + basic_survey(g) + + +if __name__ == "__main__": + main() diff --git a/4-Infrastructure/shim/pagerank_eigenvalue_survey.py b/4-Infrastructure/shim/pagerank_eigenvalue_survey.py new file mode 100644 index 00000000..0e831322 --- /dev/null +++ b/4-Infrastructure/shim/pagerank_eigenvalue_survey.py @@ -0,0 +1,229 @@ +#!/usr/bin/env python3 +""" +pagerank_eigenvalue_survey.py + +Surveys the eigenvalue sorts produced by PageRank, HITS, Fiedler, and PPR +on a synthetic social graph, then maps each sort onto the 8-bin Sidon spectral +basis from Semantics.Spectrum / Semantics.GraphRank. + +Goal: show that each classical ranking algorithm "observes" a different subset +of the 8 spectral bins, and that the bad-link gate (verifySpectralGap) is the +Sidon-filter analog of the cut identified by the Fiedler eigenvector. + +Requires: networkx, scipy, numpy + uv pip install networkx scipy numpy +""" + +import numpy as np +import networkx as nx +from scipy.sparse.linalg import eigsh +from typing import List, Dict, Tuple + +# ── Sidon spectral primitives (mirror Semantics.Spectrum) ──────────────────── + +BIN_COUNT = 8 + +def verify_spectral_gap(active: List[int]) -> bool: + """No two active bin indices may be adjacent. Mirrors verifySpectralGap.""" + return all(abs(i - j) > 1 for i in active for j in active if i != j) + +def piecewise_merge(a: List[float], b: List[float]) -> List[float]: + """Saturating superposition. Mirrors SpectralSignature.piecewiseMerge.""" + return [min(1.0, x + y) for x, y in zip(a, b)] + +def spectral_overlap(a: List[float], b: List[float]) -> float: + """Inner product. Mirrors SpectralSignature.spectralOverlap.""" + return sum(x * y for x, y in zip(a, b)) + +def bad_link(src_sig: List[float], edge_sig: List[float], + dst_sig: List[float], threshold: float = 0.05) -> bool: + merged = piecewise_merge(piecewise_merge(src_sig, edge_sig), dst_sig) + active = [i for i, v in enumerate(merged) if v > threshold] + return not verify_spectral_gap(active) + +# ── Rank mode → bin mapping (mirrors Semantics.GraphRank.RankMode) ──────────── + +def pagerank_sig(score: float, max_score: float) -> List[float]: + """PageRank mass pools into bin 0 (DC / stationary distribution).""" + sig = [0.0] * BIN_COUNT + sig[0] = score / max_score if max_score > 0 else 0.0 + return sig + +def hits_sig(auth: float, hub: float) -> List[float]: + """HITS: authority → bin 0, hub → bin 2. Gap of 2 satisfies Sidon.""" + sig = [0.0] * BIN_COUNT + sig[0] = auth + sig[2] = hub # bin 1 is skipped — maintains verifySpectralGap + return sig + +def fiedler_sig(fiedler_val: float) -> List[float]: + """Fiedler component → bin 1 (community boundary = second eigenvector).""" + sig = [0.0] * BIN_COUNT + sig[1] = min(1.0, abs(fiedler_val)) + return sig + +def ppr_sig(score: float, seed_bin: int) -> List[float]: + """PPR seeded at bin k: propagates energy near that bin.""" + sig = [0.0] * BIN_COUNT + sig[seed_bin] = score + return sig + +# ── Build test graph ────────────────────────────────────────────────────────── + +def build_test_graph() -> nx.DiGraph: + """ + Two communities with a bridge and a spam link farm. + + Community A (0-3): densely linked, should have high PageRank + Community B (4-7): densely linked, should have high PageRank + Bridge 3→4: the "bad link" candidate — crosses Sidon bin boundary + Spam farm (8-10): link farm → node 0 (artificial authority injection) + """ + g = nx.DiGraph() + for i in range(3): + for j in range(i + 1, 4): + g.add_edge(i, j) + g.add_edge(j, i) + for i in range(4, 7): + for j in range(i + 1, 8): + g.add_edge(i, j) + g.add_edge(j, i) + g.add_edge(3, 4) # inter-community bridge + for i in [8, 9, 10]: + g.add_edge(i, 0) # spam → target + g.add_edge(i, 8 if i != 8 else 10) # farm cross-links + return g + +# ── Run classical algorithms ────────────────────────────────────────────────── + +def run_pagerank(g: nx.DiGraph, alpha: float = 0.85) -> Dict[int, float]: + return nx.pagerank(g, alpha=alpha, max_iter=300) + +def run_hits(g: nx.DiGraph) -> Tuple[Dict[int, float], Dict[int, float]]: + try: + hubs, auths = nx.hits(g, max_iter=300, normalized=True) + except nx.PowerIterationFailedConvergence: + default = {n: 1.0 / len(g) for n in g} + hubs, auths = default, default + return hubs, auths + +def run_fiedler(g: nx.DiGraph) -> np.ndarray: + ug = g.to_undirected() + nodes = sorted(ug.nodes()) + L = nx.laplacian_matrix(ug, nodelist=nodes).astype(float) + try: + vals, vecs = eigsh(L, k=2, which='SM', tol=1e-8, maxiter=1000) + # Second column = Fiedler vector (index 1 after sorting by eigenvalue) + order = np.argsort(vals) + return vecs[:, order[1]], nodes + except Exception as exc: + print(f" [fiedler] eigsh failed: {exc}") + return np.zeros(len(nodes)), nodes + +# ── Survey ──────────────────────────────────────────────────────────────────── + +def run_survey(g: nx.DiGraph): + nodes = sorted(g.nodes()) + pr = run_pagerank(g) + hubs, auths = run_hits(g) + fiedler_vec, fiedler_nodes = run_fiedler(g) + fiedler_map = {fiedler_nodes[i]: fiedler_vec[i] + for i in range(len(fiedler_nodes))} + + max_pr = max(pr.values()) or 1.0 + max_auth = max(auths.values()) or 1.0 + max_hub = max(hubs.values()) or 1.0 + + print("=" * 72) + print("EIGENVALUE SORT SURVEY — mapped to 8-bin Sidon spectral basis") + print("=" * 72) + + # ── PageRank (bin-0 dominant) ────────────────────────────────────────── + pr_sorted = sorted(pr.items(), key=lambda x: -x[1]) + print("\n[PageRank] bin 0 ← DC / stationary distribution") + print(f" {'node':>5} {'PR':>8} {'bin0':>6} {'gap':>5} community") + for node, score in pr_sorted: + sig = pagerank_sig(score, max_pr) + active = [i for i, v in enumerate(sig) if v > 0.05] + gap_ok = verify_spectral_gap(active) + comm = 'A' if node <= 3 else ('B' if node <= 7 else 'spam') + print(f" {node:>5} {score:>8.5f} {sig[0]:>6.3f} {'OK' if gap_ok else 'BAD':>5} {comm}") + + # ── HITS (bins 0 + 2) ───────────────────────────────────────────────── + auth_sorted = sorted(auths.items(), key=lambda x: -x[1]) + print("\n[HITS Authority] bin 0 (auth) + bin 2 (hub); gap=2 satisfies Sidon") + print(f" {'node':>5} {'auth':>8} {'hub':>8} {'gap':>5} community") + for node, auth in auth_sorted: + hub = hubs.get(node, 0.0) + sig = hits_sig(auth / max_auth, hub / max_hub) + active = [i for i, v in enumerate(sig) if v > 0.05] + gap_ok = verify_spectral_gap(active) + comm = 'A' if node <= 3 else ('B' if node <= 7 else 'spam') + print(f" {node:>5} {auth:>8.5f} {hub:>8.5f} {'OK' if gap_ok else 'BAD':>5} {comm}") + + # ── Fiedler (bin 1) ─────────────────────────────────────────────────── + print("\n[Fiedler] bin 1 ← community boundary (second Laplacian eigenvector)") + fiedler_sorted = sorted(fiedler_map.items(), key=lambda x: x[1]) + neg_side = [n for n, v in fiedler_sorted if v < 0] + pos_side = [n for n, v in fiedler_sorted if v >= 0] + print(f" Community A (Fiedler < 0): {neg_side}") + print(f" Community B (Fiedler ≥ 0): {pos_side}") + print(f" {'node':>5} {'F_val':>9} {'bin1':>6} community") + for node, val in fiedler_sorted: + sig = fiedler_sig(val) + comm = 'A' if val < 0 else 'B' + print(f" {node:>5} {val:>+9.5f} {sig[1]:>6.3f} {comm}") + + # ── Bad-link gate on bridge and spam ────────────────────────────────── + print("\n[Bad-link gate]") + + # Bridge 3→4 + src_sig = pagerank_sig(pr[3], max_pr) + dst_sig = pagerank_sig(pr[4], max_pr) + edge_sig = [0.5] + [0.0] * 7 # bridge carries half bin-0 authority + merged = piecewise_merge(piecewise_merge(src_sig, edge_sig), dst_sig) + active_m = [i for i, v in enumerate(merged) if v > 0.05] + is_bad = not verify_spectral_gap(active_m) + print(f" Bridge 3→4: merged active bins={active_m} " + f"gap={'FAIL (bad link)' if is_bad else 'OK'}") + + for spam_node in [8, 9, 10]: + s_sig = pagerank_sig(pr[spam_node], max_pr) + t_sig = pagerank_sig(pr[0], max_pr) + e_sig = [0.5] + [0.0] * 7 + merged = piecewise_merge(piecewise_merge(s_sig, e_sig), t_sig) + active_m = [i for i, v in enumerate(merged) if v > 0.05] + is_bad = not verify_spectral_gap(active_m) + print(f" Spam {spam_node}→0: merged active bins={active_m} " + f"gap={'FAIL (bad link)' if is_bad else 'OK'}") + + # ── Spectral overlap as PPR dot product ─────────────────────────────── + print("\n[Spectral overlap — PPR dot product with seed = community A sig]") + seed_sig = pagerank_sig(1.0, 1.0) # pure bin-0 seed + scores = [] + for node in nodes: + sig = pagerank_sig(pr[node], max_pr) + score = spectral_overlap(sig, seed_sig) + scores.append((node, score)) + scores.sort(key=lambda x: -x[1]) + print(f" {'node':>5} {'overlap':>9} community") + for node, sc in scores: + comm = 'A' if node <= 3 else ('B' if node <= 7 else 'spam') + print(f" {node:>5} {sc:>9.5f} {comm}") + + # ── Summary ─────────────────────────────────────────────────────────── + print("\n[Summary — eigenvalue sort → spectral bin mapping]") + print(" PageRank → bin 0 (DC authority; bad link = bin-0 bleed)") + print(" HITS auth → bin 0 (top singular vector = pure authority)") + print(" HITS hub → bin 2 (gap=2 from authority; Sidon-safe)") + print(" Fiedler → bin 1 (community cut = spectral gap)") + print(" PPR seed → bin k (propagates along matching-bin paths)") + print(" CheiRank → 1-bin0 (givers have LOW bin-0; receivers HIGH)") + print(" verifySpectralGap ← detects ALL of the above failures at once") + print(" erdosHooleyDelta ← density threshold for Sidon property in 8-bin space") + + +if __name__ == "__main__": + g = build_test_graph() + print(f"Graph: {len(g.nodes())} nodes, {len(g.edges())} edges\n") + run_survey(g) diff --git a/4-Infrastructure/shim/s2_citation_crawler.py b/4-Infrastructure/shim/s2_citation_crawler.py new file mode 100644 index 00000000..6cca2172 --- /dev/null +++ b/4-Infrastructure/shim/s2_citation_crawler.py @@ -0,0 +1,333 @@ +#!/usr/bin/env python3 +""" +s2_citation_crawler.py + +BFS citation graph crawler for Semantic Scholar via the s2-proxy Cloudflare +Worker fleet. Feeds the result into pagerank_eigenvalue_survey.run_survey(). + +Architecture: + coordinator (this script) + └─ round-robins requests across N worker URLs + each worker: KV cache → S2 API (with per-worker key) → KV cache write + +Usage: + # Single worker (no key, slow but free): + python s2_citation_crawler.py --seed "PageRank Brin Page 1998" --hops 2 + + # Worker fleet (fast): + python s2_citation_crawler.py \\ + --seed "PageRank Brin Page 1998" \\ + --hops 2 \\ + --workers https://s2-proxy.YOUR.workers.dev https://s2-proxy-2.YOUR.workers.dev \\ + --out citation_graph.json + +Requires: requests, networkx, scipy, numpy + uv pip install requests networkx scipy numpy +""" + +import argparse +import json +import time +import sys +from collections import deque +from dataclasses import dataclass, asdict, field +from pathlib import Path +from typing import List, Dict, Optional, Set, Iterator +from itertools import cycle + +import requests +import networkx as nx + +# ── Import spectral survey from sibling script ──────────────────────────────── +sys.path.insert(0, str(Path(__file__).parent)) +try: + from pagerank_eigenvalue_survey import run_survey, BIN_COUNT, verify_spectral_gap + HAS_SURVEY = True +except ImportError: + HAS_SURVEY = False + print("[warn] pagerank_eigenvalue_survey not found — graph will be saved but not analysed") + +# ── S2 field set ────────────────────────────────────────────────────────────── + +PAPER_FIELDS = "paperId,title,year,referenceCount,citationCount,influentialCitationCount" +REF_FIELDS = "paperId,title,year,citationCount" +CITE_FIELDS = "paperId,title,year,citationCount" + +# ── Data types ──────────────────────────────────────────────────────────────── + +@dataclass +class S2Paper: + paper_id: str + title: str + year: Optional[int] + ref_count: int + cite_count: int + influential_cite_count: int + references: List[str] = field(default_factory=list) # paper_ids + citations: List[str] = field(default_factory=list) # paper_ids + +# ── Worker fleet client ─────────────────────────────────────────────────────── + +class WorkerFleet: + """Round-robins requests across the s2-proxy worker fleet.""" + + DIRECT_BASE = "https://api.semanticscholar.org" + + def __init__(self, worker_urls: List[str], api_key: Optional[str] = None, + rate_delay: float = 1.1): + self.worker_urls = worker_urls + self.api_key = api_key + self.rate_delay = rate_delay # seconds between requests (per worker slot) + self._fleet = cycle(worker_urls) if worker_urls else None + self._last_req = 0.0 + self._stats = {"hits": 0, "misses": 0, "errors": 0, "requests": 0} + + def _next_url_base(self) -> str: + if self._fleet: + return next(self._fleet) + return self.DIRECT_BASE + + def _throttle(self): + gap = time.monotonic() - self._last_req + if gap < self.rate_delay: + time.sleep(self.rate_delay - gap) + self._last_req = time.monotonic() + + def get(self, s2_path: str, params: Dict = None) -> Optional[dict]: + self._throttle() + base = self._next_url_base() + self._stats["requests"] += 1 + + try: + if self._fleet: + # Via worker + worker_params = {"path": s2_path, **(params or {})} + resp = requests.get(f"{base}/fetch", params=worker_params, timeout=15) + cache_status = resp.headers.get("X-S2-Cache", "?") + if cache_status == "HIT": + self._stats["hits"] += 1 + else: + self._stats["misses"] += 1 + else: + # Direct S2 + headers = {} + if self.api_key: + headers["x-api-key"] = self.api_key + resp = requests.get( + f"{self.DIRECT_BASE}{s2_path}", + params=params or {}, + headers=headers, + timeout=15, + ) + + # Back-off on rate limit + if resp.status_code == 429: + retry_after = int(resp.headers.get("retry-after", 60)) + print(f" [rate-limit] waiting {retry_after}s") + time.sleep(retry_after) + return self.get(s2_path, params) # retry once + + resp.raise_for_status() + return resp.json() + + except Exception as exc: + self._stats["errors"] += 1 + print(f" [fetch error] {s2_path}: {exc}") + return None + + def stats(self) -> str: + s = self._stats + hit_rate = (s["hits"] / max(s["requests"], 1)) * 100 + return (f"requests={s['requests']} cache_hits={s['hits']}({hit_rate:.0f}%)" + f" misses={s['misses']} errors={s['errors']}") + + +# ── S2 API wrappers ─────────────────────────────────────────────────────────── + +def search_paper(fleet: WorkerFleet, query: str) -> Optional[str]: + """Return the top S2 paper_id for a title/author query.""" + data = fleet.get("/graph/v1/paper/search", { + "query": query, "fields": "paperId,title,year", "limit": 1 + }) + if not data or not data.get("data"): + return None + hit = data["data"][0] + print(f" seed paper: [{hit.get('year','')}] {hit.get('title','?')}") + return hit["paperId"] + + +def fetch_paper(fleet: WorkerFleet, paper_id: str) -> Optional[S2Paper]: + data = fleet.get(f"/graph/v1/paper/{paper_id}", {"fields": PAPER_FIELDS}) + if not data or "paperId" not in data: + return None + return S2Paper( + paper_id=data["paperId"], + title=data.get("title") or "", + year=data.get("year"), + ref_count=data.get("referenceCount", 0), + cite_count=data.get("citationCount", 0), + influential_cite_count=data.get("influentialCitationCount", 0), + ) + + +def fetch_neighbors(fleet: WorkerFleet, paper_id: str, + max_refs: int = 50, max_cites: int = 50 + ) -> tuple[List[str], List[str]]: + """Return (reference_ids, citation_ids) for a paper, capped at max_*.""" + refs, cites = [], [] + + ref_data = fleet.get(f"/graph/v1/paper/{paper_id}/references", { + "fields": REF_FIELDS, "limit": max_refs + }) + if ref_data and "data" in ref_data: + refs = [e["citedPaper"]["paperId"] + for e in ref_data["data"] + if e.get("citedPaper", {}).get("paperId")] + + cite_data = fleet.get(f"/graph/v1/paper/{paper_id}/citations", { + "fields": CITE_FIELDS, "limit": max_cites + }) + if cite_data and "data" in cite_data: + cites = [e["citingPaper"]["paperId"] + for e in cite_data["data"] + if e.get("citingPaper", {}).get("paperId")] + + return refs, cites + + +# ── BFS crawler ─────────────────────────────────────────────────────────────── + +def crawl(fleet: WorkerFleet, seed_id: str, hops: int, + max_refs: int = 40, max_cites: int = 40, + max_papers: int = 500) -> Dict[str, S2Paper]: + """BFS from seed_id for `hops` hops. Returns dict of paper_id → S2Paper.""" + visited: Set[str] = set() + papers: Dict[str, S2Paper] = {} + queue: deque = deque([(seed_id, 0)]) + + print(f"\nBFS crawl: seed={seed_id} hops={hops} max_papers={max_papers}") + + while queue and len(papers) < max_papers: + pid, depth = queue.popleft() + if pid in visited: + continue + visited.add(pid) + + paper = fetch_paper(fleet, pid) + if not paper: + continue + papers[pid] = paper + print(f" [{depth}/{hops}] {len(papers):>4} papers " + f"cite={paper.cite_count} {paper.title[:60]}") + + if depth < hops: + refs, cites = fetch_neighbors(fleet, pid, max_refs, max_cites) + paper.references = refs + paper.citations = cites + for nid in refs + cites: + if nid not in visited: + queue.append((nid, depth + 1)) + + print(f"\nCrawl complete: {len(papers)} papers, {fleet.stats()}") + return papers + + +# ── Build NetworkX graph ────────────────────────────────────────────────────── + +def build_graph(papers: Dict[str, S2Paper]) -> nx.DiGraph: + """Citation graph: edge A→B means A cites B.""" + g = nx.DiGraph() + pid_to_int = {pid: i for i, pid in enumerate(papers)} + + for pid, paper in papers.items(): + node_id = pid_to_int[pid] + g.add_node(node_id, + paper_id=pid, + title=paper.title[:80], + year=paper.year or 0, + cite_count=paper.cite_count) + + for pid, paper in papers.items(): + src = pid_to_int[pid] + for ref_id in paper.references: + if ref_id in pid_to_int: + g.add_edge(src, pid_to_int[ref_id]) # src cites ref + + print(f"Graph: {g.number_of_nodes()} nodes, {g.number_of_edges()} edges") + return g + + +# ── Main ────────────────────────────────────────────────────────────────────── + +def main(): + ap = argparse.ArgumentParser(description="S2 citation BFS crawler + spectral survey") + ap.add_argument("--seed", default="PageRank Brin Page 1998 web search engine", + help="Title/author query to find seed paper") + ap.add_argument("--seed-id", default=None, + help="Use a specific S2 paper_id instead of title search") + ap.add_argument("--hops", type=int, default=2, + help="BFS depth (default 2; 3 gets very large)") + ap.add_argument("--max-refs", type=int, default=40, + help="Max references to follow per paper") + ap.add_argument("--max-cites", type=int, default=40, + help="Max citations to follow per paper") + ap.add_argument("--max-papers", type=int, default=300, + help="Hard cap on total papers fetched") + ap.add_argument("--workers", nargs="*", default=[], + help="s2-proxy worker base URLs (omit to hit S2 directly)") + ap.add_argument("--api-key", default=None, + help="S2 API key for direct mode (not used in worker mode)") + ap.add_argument("--rate-delay", type=float, default=1.1, + help="Seconds between requests per worker slot") + ap.add_argument("--out", default="citation_graph.json", + help="Output file for the paper graph JSON") + ap.add_argument("--no-survey", action="store_true", + help="Skip spectral survey (just crawl and save)") + args = ap.parse_args() + + fleet = WorkerFleet(args.workers, api_key=args.api_key, rate_delay=args.rate_delay) + + # Resolve seed + if args.seed_id: + seed_id = args.seed_id + else: + print(f"Searching for seed: {args.seed}") + seed_id = search_paper(fleet, args.seed) + if not seed_id: + print("ERROR: seed paper not found") + sys.exit(1) + print(f"Seed paper_id: {seed_id}") + + # Crawl + papers = crawl(fleet, seed_id, + hops=args.hops, + max_refs=args.max_refs, + max_cites=args.max_cites, + max_papers=args.max_papers) + + # Save + out = Path(args.out) + with out.open("w") as f: + json.dump({pid: asdict(p) for pid, p in papers.items()}, f, indent=2) + print(f"Saved {len(papers)} papers → {out}") + + # Build graph and run spectral survey + g = build_graph(papers) + if HAS_SURVEY and not args.no_survey: + print("\n" + "=" * 72) + print("Running spectral survey on real S2 citation graph...") + print("=" * 72) + run_survey(g) + else: + # Basic stats without survey + pr = nx.pagerank(g, alpha=0.85) + top10 = sorted(pr.items(), key=lambda x: -x[1])[:10] + print("\nTop 10 by PageRank:") + for nid, score in top10: + title = g.nodes[nid].get("title", "?") + year = g.nodes[nid].get("year", "?") + print(f" [{year}] {score:.5f} {title}") + + +if __name__ == "__main__": + main() diff --git a/4-Infrastructure/shim/snap_pist_spectral_crossval.py b/4-Infrastructure/shim/snap_pist_spectral_crossval.py new file mode 100644 index 00000000..46c6bdff --- /dev/null +++ b/4-Infrastructure/shim/snap_pist_spectral_crossval.py @@ -0,0 +1,444 @@ +#!/usr/bin/env python3 +""" +snap_pist_spectral_crossval — Cross-validate SNAP citation graph through PIST +matrix projection vs spectral survey. + +Pipeline: +1. Load SNAP edge list → NetworkX DiGraph +2. (Optionally) sample subgraph for spectral analysis +3. Serialize each paper's citation neighborhood as text +4. Run PIST token→strand 8×8 matrix projection +5. Run spectral survey (PageRank, HITS, Fiedler) +6. Cross-tabulate PIST matrix hashes vs spectral bin assignments +7. Emit comparison receipt JSON + +Claim boundary: cross-validation comparison only; no classifier output. +""" + +import gzip +import hashlib +import json +import math +import os +import random +import re +import sys +from collections import defaultdict, Counter + +import networkx as nx +import numpy as np + +N_STRANDS = 8 +random.seed(42) +np.random.seed(42) + +# ═══════════════════════════════════════════════════════════════════════ +# §1 SNAP PARSING +# ═══════════════════════════════════════════════════════════════════════ + +def load_snap_edgelist(path: str) -> nx.DiGraph: + g = nx.DiGraph() + with gzip.open(path, "rt") if path.endswith(".gz") else open(path) as f: + for line in f: + line = line.strip() + if not line or line.startswith("#"): + continue + parts = line.split() + if len(parts) < 2: + continue + g.add_edge(int(parts[0]), int(parts[1])) + return g + + +# ═══════════════════════════════════════════════════════════════════════ +# §2 CITATION CONTEXT TEXT + PIST TOKENIZER +# ═══════════════════════════════════════════════════════════════════════ + +def tokenize(text: str) -> list[str]: + """PIST-compatible tokenizer: split on space, _, and : (added space for citation IDs).""" + normalized = text.replace("-", "_") + return [p for p in re.split(r"[ _:]", normalized) if p] + + +def node_context_text(g: nx.DiGraph, node: int) -> str: + """Serialize a paper's citation neighborhood as space-separated text. + + Format: paper_id cited_id1 cited_id2 ... cited_by_id1 cited_by_id2 ... + + After space-splitting tokenization, adjacent paper IDs become bigram + pairs → mapped to 8×8 strand adjacency in PIST space. + """ + parts = [str(node)] + parts.extend(str(v) for v in g.successors(node)) + parts.extend(str(v) for v in g.predecessors(node)) + return " ".join(parts) + + +# ═══════════════════════════════════════════════════════════════════════ +# §3 PIST 8×8 ADJACENCY MATRIX +# ═══════════════════════════════════════════════════════════════════════ + +def build_global_vocabulary(all_texts: list[str]) -> dict[str, int]: + all_tokens = set() + for text in all_texts: + for t in tokenize(text): + all_tokens.add(t) + return {t: i for i, t in enumerate(sorted(all_tokens))} + + +def strand_for_token(token: str, vocab: dict[str, int]) -> int: + return vocab[token] % N_STRANDS + + +def build_adjacency_matrix(tokens: list[str], vocab: dict[str, int]) -> list[list[int]]: + M = [[0] * N_STRANDS for _ in range(N_STRANDS)] + strands = [strand_for_token(t, vocab) for t in tokens] + for i in range(len(strands) - 1): + s1, s2 = strands[i], strands[i + 1] + M[s1][s2] += 1 + return M + + +def matrix_hash(M: list[list[int]]) -> str: + return hashlib.sha256( + json.dumps(M, separators=(",", ":")).encode("utf-8") + ).hexdigest() + + +# ═══════════════════════════════════════════════════════════════════════ +# §4 SPECTRAL BIN ASSIGNMENT +# ═══════════════════════════════════════════════════════════════════════ + +def classify_pagerank_bin(score: float, max_score: float) -> str: + if score >= 0.1 * max_score: return "authority" + elif score >= 0.01 * max_score: return "secondary" + return "noise" + + +def classify_hits_bin(auth: float, max_auth: float, hub: float, max_hub: float) -> str: + auth_norm = auth / max_auth if max_auth else 0 + hub_norm = hub / max_hub if max_hub else 0 + if auth_norm > 0.1 and auth_norm > hub_norm: return "authority" + elif hub_norm > 0.1: return "hub" + elif auth_norm > 0.01 or hub_norm > 0.01: return "secondary" + return "noise" + + +def classify_fiedler_bin(val: float) -> str: + return "community_A" if val < 0 else "community_B" + + +def spectral_bin_summary(pr_bin: str, hits_bin: str) -> str: + return f"{pr_bin}|{hits_bin}" + + +# ═══════════════════════════════════════════════════════════════════════ +# §5 FIEDLER ON DEGREE-WEIGHTED SUBSAMPLE + PAGERANK-KNN EXTRAPOLATION +# ═══════════════════════════════════════════════════════════════════════ + +def compute_fiedler_full(g: nx.DiGraph) -> tuple[dict[int, str], dict[int, float], dict]: + """Compute Fiedler vector on a random 2-hop sampled subgraph, + then extend to all nodes via PageRank-kNN extrapolation. + + Strategy: + 1. Pick 500 random seeds, expand 2 hops for connectivity + 2. Compute Fiedler via small eigsh (which='SM') + 3. Extrapolate to all 34K nodes via kNN in PageRank space + (PageRank is an eigenvector centrality that preserves spectral + ordering; random sampling avoids degree bias in the computed set.) + """ + from scipy.sparse.linalg import eigsh + + nodes = sorted(g.nodes()) + ug = g.to_undirected() + n_total = len(nodes) + + # 1. Random 2-hop sampling (representative, no degree bias) + seeds = random.sample(nodes, min(500, n_total)) + sample_set: set[int] = set(seeds) + for s in seeds: + sample_set.update(g.successors(s)) + sample_set.update(g.predecessors(s)) + for s2 in list(g.successors(s)) + list(g.predecessors(s)): + sample_set.update(g.successors(s2)) + sample_set.update(g.predecessors(s2)) + sample = sorted(sample_set & set(nodes)) + print(f" Fiedler: {len(seeds)} seeds, 2-hop → {len(sample)} nodes...", flush=True) + + # Build undirected subgraph, take largest connected component + ug_s = g.subgraph(sample).to_undirected() + cc = max(nx.connected_components(ug_s), key=len) + ug_s = ug_s.subgraph(cc) + cc_nodes = sorted(cc) + print(f" Fiedler component: {len(cc_nodes)} nodes, " + f"{ug_s.number_of_edges()} edges", flush=True) + + # 2. Compute Fiedler on the subgraph + L_s = nx.laplacian_matrix(ug_s, nodelist=cc_nodes).astype(float) + v0_s = np.array([ug_s.degree(n) for n in cc_nodes], dtype=float) + v0_s = v0_s / np.linalg.norm(v0_s) + + vals, vecs = eigsh(L_s, k=2, which="SM", maxiter=10000, tol=1e-4, v0=v0_s) + order = np.argsort(vals) + fiedler_val = vals[order[1]] + fiedler_vec_real = vecs[:, order[1]] + sample_f = {cc_nodes[i]: float(fiedler_vec_real[i]) for i in range(len(cc_nodes))} + print(f" Fiedler eigenvalue: {fiedler_val:.6f} " + f"zero eig: {vals[order[0]]:.6e}", flush=True) + + # 3. Extend to all nodes via PageRank-kNN + pr = nx.pagerank(g, alpha=0.85, max_iter=300) + pr_vals = {n: pr[n] for n in nodes} + sample_keys = set(sample_f.keys()) + + fiedler_bins: dict[int, str] = {} + fiedler_values: dict[int, float] = {} + k = 5 + for n in nodes: + if n in sample_f: + fv = sample_f[n] + else: + nearest = sorted(sample_keys, + key=lambda sn: abs(pr_vals.get(n, 0.0) - pr_vals.get(sn, 0.0)))[:k] + signs = [1 if sample_f[sn] < 0 else 0 for sn in nearest] + majority = sum(signs) / len(signs) if nearest else 0.5 + fv = -1.0 if majority > 0.5 else 1.0 + fiedler_values[n] = fv + fiedler_bins[n] = classify_fiedler_bin(fv) + + n_extrapolated = len(nodes) - len(sample_f) + print(f" Fiedler extension: {len(sample_f)} computed, " + f"{n_extrapolated} PageRank-kNN extrapolated", flush=True) + + n_neg = sum(1 for v in sample_f.values() if v < 0) + n_pos = len(sample_f) - n_neg + collapse_report = { + "n_total": len(sample_f), + "n_neg": n_neg, + "n_pos": n_pos, + "collapse_ratio": max(n_neg, n_pos) / len(sample_f) if len(sample_f) else 0.0, + "var_sign": float(np.var([1 if v < 0 else 0 for v in sample_f.values()])), + "fiedler_eigenvalue": round(fiedler_val, 6), + "component_nodes": len(sample_f), + "total_nodes": len(nodes), + } + + return fiedler_bins, fiedler_values, collapse_report + + +# ═══════════════════════════════════════════════════════════════════════ +# §6 CROSS-TABULATION +# ═══════════════════════════════════════════════════════════════════════ + +def cross_tabulate( + pist_hash_to_nodes: dict[str, list[int]], + node_spectral: dict[int, dict], +) -> dict: + tab = {} + for phash, nodes in pist_hash_to_nodes.items(): + combined_bins = Counter() + pr_bins = Counter() + hits_bins = Counter() + for n in nodes: + s = node_spectral.get(n, {}) + combined_bins[s.get("spectral_summary", "unknown")] += 1 + pr_bins[s.get("pr_bin", "unknown")] += 1 + hits_bins[s.get("hits_bin", "unknown")] += 1 + total = len(nodes) + + def _entropy(cnt: Counter) -> float: + if total == 0: + return 0.0 + return -sum((c / total) * math.log2(c / total) for c in cnt.values() if c > 0) + + tab[phash] = { + "count": total, + "pagerank_mode": pr_bins.most_common(1)[0][0] if pr_bins else None, + "hits_mode": hits_bins.most_common(1)[0][0] if hits_bins else None, + "spectral_entropy": round(_entropy(combined_bins), 4), + } + return tab + + +# ═══════════════════════════════════════════════════════════════════════ +# §7 MAIN +# ═══════════════════════════════════════════════════════════════════════ + +def main(): + path = sys.argv[1] if len(sys.argv) > 1 else "/tmp/snap_data/cit-HepPh.txt.gz" + + print(f"Loading {path}...", flush=True) + g = load_snap_edgelist(path) + print(f"Graph: {g.number_of_nodes()} nodes, {g.number_of_edges()} edges", flush=True) + + nodes = sorted(g.nodes()) + + # ── 1. Serialize each paper as citation context text ─────────────── + node_texts = {n: node_context_text(g, n) for n in nodes} + all_texts = list(node_texts.values()) + + # ── 2. Build PIST global vocabulary ──────────────────────────────── + vocab = build_global_vocabulary(all_texts) + gvh = hashlib.sha256("|".join(vocab.keys()).encode("utf-8")).hexdigest() + print(f"Vocabulary: {len(vocab)} tokens hash={gvh}", flush=True) + + # ── 3. Compute PIST 8×8 matrix per paper ─────────────────────────── + node_pist = {} + for n in nodes: + tokens = tokenize(node_texts[n]) + M = build_adjacency_matrix(tokens, vocab) + node_pist[n] = {"tokens": tokens, "matrix_8x8": M, "matrix_hash": matrix_hash(M)} + + pist_hash_groups = defaultdict(list) + for n, info in node_pist.items(): + pist_hash_groups[info["matrix_hash"]].append(n) + + n_unique = len(pist_hash_groups) + print(f"Unique PIST matrix hashes: {n_unique} from {len(nodes)} papers", flush=True) + hash_sizes = [len(v) for v in pist_hash_groups.values()] + print(f" Hash group sizes: min={min(hash_sizes)} max={max(hash_sizes)} " + f"median={sorted(hash_sizes)[len(hash_sizes)//2]}", flush=True) + + # Show a few sample matrices for verification + for n in nodes[:3]: + info = node_pist[n] + flat = [str(c) for row in info["matrix_8x8"] for c in row] + print(f" node={n} tokens={info['tokens'][:5]}... " + f"matrix_sum={sum(int(c) for c in flat)} " + f"nonzero={sum(1 for c in flat if int(c) > 0)}", flush=True) + + # ── 4. Run spectral survey ───────────────────────────────────────── + print("PageRank...", flush=True) + pr = nx.pagerank(g, alpha=0.85, max_iter=300) + + print("HITS...", flush=True) + try: + hubs, auths = nx.hits(g, max_iter=300, normalized=True) + except nx.PowerIterationFailedConvergence: + hubs = auths = {n: 1.0 / len(g) for n in g} + + # Fiedler: diagnostic only (excluded from spectral scoring). + # Under k-hop projection, the Fiedler sign collapses to the empirical + # majority class (94.5% single-sign in our regime). See + # docs/fiedler_non_identifiability.md for the operator-theoretic + # statement of this non-identifiability condition. + print("Fiedler (diagnostic)...", flush=True) + fiedler_bins, fiedler_values, collapse_report = compute_fiedler_full(g) + collapse_ratio = ( + max(collapse_report["n_neg"], collapse_report["n_pos"]) + / collapse_report["n_total"] + ) if collapse_report["n_total"] else 0.0 + print(f" Fiedler collapse ratio: {collapse_ratio:.4f} " + f"(Var(sign)) = {collapse_report['var_sign']:.4f})", flush=True) + + max_pr = max(pr.values()) or 1.0 + max_auth = max(auths.values()) or 1.0 + max_hub = max(hubs.values()) or 1.0 + + node_spectral = {} + for n in nodes: + pr_b = classify_pagerank_bin(pr[n], max_pr) + hits_b = classify_hits_bin(auths.get(n, 0), max_auth, hubs.get(n, 0), max_hub) + node_spectral[n] = { + "pr": round(pr[n], 6), "pr_bin": pr_b, + "auth": round(auths.get(n, 0), 6), + "hub": round(hubs.get(n, 0), 6), "hits_bin": hits_b, + "fiedler": round(fiedler_values.get(n, 0.0), 6), + "fiedler_bin": fiedler_bins.get(n, "community_unknown"), + "spectral_summary": spectral_bin_summary(pr_b, hits_b), + } + + spectral_bin_counts = Counter(v["spectral_summary"] for v in node_spectral.values()) + print(f"Spectral combos (PR|HITS): {len(spectral_bin_counts)}", flush=True) + for label, count in spectral_bin_counts.most_common(10): + print(f" {label}: {count}", flush=True) + + # ── 5. Cross-tabulate ────────────────────────────────────────────── + tab = cross_tabulate(pist_hash_groups, node_spectral) + pure_groups = sum(1 for v in tab.values() if v["spectral_entropy"] == 0.0) + mixed_groups = len(tab) - pure_groups + print(f"\nPIST hash groups pure spectral: {pure_groups}/{len(tab)} " + f"mixed: {mixed_groups}/{len(tab)}", flush=True) + + # Mean entropy across all groups + mean_entropy = np.mean([v["spectral_entropy"] for v in tab.values()]) + print(f"Mean spectral entropy across PIST groups: {mean_entropy:.4f}", flush=True) + + # Diagnose mixing within multi-paper PIST groups. + pr_pure = 0 + hits_pure = 0 + for phash, nodes_p in pist_hash_groups.items(): + if len(nodes_p) < 2: + continue + pr_bins = set(node_spectral[n]["pr_bin"] for n in nodes_p) + hits_bins = set(node_spectral[n]["hits_bin"] for n in nodes_p) + if len(pr_bins) == 1: pr_pure += 1 + if len(hits_bins) == 1: hits_pure += 1 + + multi_paper_groups = sum(1 for v in pist_hash_groups.values() if len(v) > 1) + print(f"\nMulti-paper PIST groups: {multi_paper_groups}", flush=True) + print(f" PR-bin pure: {pr_pure}/{multi_paper_groups} " + f"({100*pr_pure/multi_paper_groups:.1f}%)", flush=True) + print(f" HITS-bin pure: {hits_pure}/{multi_paper_groups} " + f"({100*hits_pure/multi_paper_groups:.1f}%)", flush=True) + print(f" Fiedler: excluded from scoring " + f"(collapse_ratio={collapse_ratio:.4f}, " + f"Var(sign)={collapse_report['var_sign']:.4f})", flush=True) + + # ── 6. Output ────────────────────────────────────────────────────── + predictions = [] + for n in nodes: + info = node_pist.get(n, {}) + spec = node_spectral.get(n, {}) + predictions.append({ + "node_id": n, + "matrix_hash": info.get("matrix_hash", ""), + "pagerank": spec.get("pr", 0), "pr_bin": spec.get("pr_bin", "unknown"), + "auth": spec.get("auth", 0), "hub": spec.get("hub", 0), + "hits_bin": spec.get("hits_bin", "unknown"), + "fiedler": spec.get("fiedler", 0), "fiedler_bin": spec.get("fiedler_bin", "unknown"), + "spectral_summary": spec.get("spectral_summary", "unknown"), + }) + + artifact = { + "schema": "snap_pist_spectral_crossval_v2", + "claim_boundary": "cross-validation-only;pr+hits-scoring;fiedler-diagnostic-only", + "source": os.path.basename(path), + "graph": {"nodes": g.number_of_nodes(), "edges": g.number_of_edges()}, + "pist": { + "global_vocab_hash": gvh, + "unique_matrix_hashes": n_unique, + "matrix_schema": "token_strand_adjacency_8x8_v1", + }, + "spectral": { + "distinct_combinations": len(spectral_bin_counts), + "fiedler_collapse": collapse_report, + }, + "cross_tabulation": { + "pist_groups_with_pure_spectral_bins": pure_groups, + "pist_groups_with_mixed_spectral_bins": mixed_groups, + "total_pist_groups": len(tab), + "purity_ratio": round(pure_groups / len(tab), 4) if tab else 0.0, + "mean_spectral_entropy": round(mean_entropy, 4), + "scoring_dimensions": "pr_bin|hits_bin (fiedler excluded)", + }, + } + + out_path = f"/tmp/snap_pist_spectral_crossval_{os.path.basename(path)}.json" + with open(out_path, "w") as f: + json.dump(artifact, f, indent=2, sort_keys=True) + print(f"\nWrote {out_path}", flush=True) + + # ── 7. Top PIST groups by size ───────────────────────────────────── + print("\nTop PIST groups by size:", flush=True) + for phash, pnodes in sorted(pist_hash_groups.items(), key=lambda x: -len(x[1]))[:25]: + info = tab[phash] + spec = node_spectral.get(pnodes[0], {}) + print(f" PIST={phash[:16]} count={len(pnodes):>4} " + f"PR={spec.get('pr_bin','?'):>10} HITS={spec.get('hits_bin','?'):>10} " + f"H={info['spectral_entropy']:.3f}", flush=True) + + return 0 + +if __name__ == "__main__": + sys.exit(main()) diff --git a/4-Infrastructure/shim/validate_fiedler.py b/4-Infrastructure/shim/validate_fiedler.py new file mode 100644 index 00000000..354a1290 --- /dev/null +++ b/4-Infrastructure/shim/validate_fiedler.py @@ -0,0 +1,438 @@ +#!/usr/bin/env python3 +""" +validate_fiedler — Sanity checks for Fiedler extrapolation in SNAP crossval. + +Checks: +1. Internal cross-validation: withhold 20% of computed nodes, PR-kNN predict sign +2. No-propagation baseline: global majority Fiedler assignment +3. Shuffle test: degree-preserving randomization, re-run purity + +Usage: python3 validate_fiedler.py /tmp/cit-HepPh.txt.gz +""" + +import gzip +import hashlib +import json +import math +import os +import random +import sys +from collections import defaultdict, Counter + +import networkx as nx +import numpy as np +from scipy.sparse.linalg import eigsh + +random.seed(42) +np.random.seed(42) + + +def load_graph(path: str) -> nx.DiGraph: + g = nx.DiGraph() + with gzip.open(path, "rt") if path.endswith(".gz") else open(path) as f: + for line in f: + line = line.strip() + if not line or line.startswith("#"): + continue + parts = line.split() + if len(parts) < 2: + continue + g.add_edge(int(parts[0]), int(parts[1])) + return g + + +def compute_fiedler_component(g: nx.DiGraph): + """Return (computed_fiedler, all_nodes, ug) — same as crossval's random 2-hop.""" + nodes = sorted(g.nodes()) + ug = g.to_undirected() + n_total = len(nodes) + + seeds = random.sample(nodes, min(500, n_total)) + sample_set: set[int] = set(seeds) + for s in seeds: + sample_set.update(g.successors(s)) + sample_set.update(g.predecessors(s)) + for s2 in list(g.successors(s)) + list(g.predecessors(s)): + sample_set.update(g.successors(s2)) + sample_set.update(g.predecessors(s2)) + sample = sorted(sample_set & set(nodes)) + + ug_s = g.subgraph(sample).to_undirected() + cc = max(nx.connected_components(ug_s), key=len) + ug_s = ug_s.subgraph(cc) + cc_nodes = sorted(cc) + + L_s = nx.laplacian_matrix(ug_s, nodelist=cc_nodes).astype(float) + v0_s = np.array([ug_s.degree(n) for n in cc_nodes], dtype=float) + v0_s = v0_s / np.linalg.norm(v0_s) + + vals, vecs = eigsh(L_s, k=2, which="SM", maxiter=10000, tol=1e-4, v0=v0_s) + order = np.argsort(vals) + fiedler_val = vals[order[1]] + fiedler_vec = vecs[:, order[1]] + sample_f = {cc_nodes[i]: float(fiedler_vec[i]) for i in range(len(cc_nodes))} + + print(f" Fiedler eigenvalue: {fiedler_val:.6f} zero eig: {vals[order[0]]:.6e}") + print(f" Component: {len(cc_nodes)} nodes, {len(nodes) - len(cc_nodes)} uncomputed") + return sample_f, nodes, ug + + +def internal_crossval(sample_f: dict[int, float], nodes: list[int], ug: nx.Graph): + """Withhold 20% of computed nodes, PR-kNN predict Fiedler sign.""" + print("\n═══ CHECK 1: Internal cross-validation ═══") + pr = nx.pagerank(ug.to_directed(), alpha=0.85, max_iter=300) + pr_vals = {n: pr.get(n, 0.0) for n in nodes} + + computed = list(sample_f.keys()) + random.shuffle(computed) + split = int(len(computed) * 0.8) + train_set = set(computed[:split]) + held_out = computed[split:] + + train_f = {n: sample_f[n] for n in train_set} + train_keys = set(train_f.keys()) + + correct = 0 + fv_actual = [] + fv_predicted = [] + k = 5 + for n in held_out: + nearest = sorted(train_keys, + key=lambda sn: abs(pr_vals.get(n, 0.0) - pr_vals.get(sn, 0.0)))[:k] + signs = [1 if train_f[sn] < 0 else 0 for sn in nearest] + majority = sum(signs) / len(signs) if nearest else 0.5 + pred_fv = -1.0 if majority > 0.5 else 1.0 + actual_fv = sample_f[n] + actual_sign = -1.0 if actual_fv < 0 else 1.0 + if (pred_fv < 0) == (actual_fv < 0): + correct += 1 + fv_actual.append(actual_sign) + fv_predicted.append(pred_fv) + + accuracy = correct / len(held_out) + + # Cosine similarity between actual and predicted Fiedler vectors + a = np.array(fv_actual) + p = np.array(fv_predicted) + cos_sim = np.dot(a, p) / (np.linalg.norm(a) * np.linalg.norm(p)) + + # Majority class baseline + n_neg = sum(1 for v in train_f.values() if v < 0) + maj_class = -1 if n_neg > len(train_f) - n_neg else 1 + maj_accuracy = sum(1 for n in held_out if (sample_f[n] < 0) == (maj_class < 0)) / len(held_out) + + # Confusion matrix + tp = sum(1 for i in range(len(held_out)) if fv_actual[i] < 0 and fv_predicted[i] < 0) + fp = sum(1 for i in range(len(held_out)) if fv_actual[i] > 0 and fv_predicted[i] < 0) + tn = sum(1 for i in range(len(held_out)) if fv_actual[i] > 0 and fv_predicted[i] > 0) + fn = sum(1 for i in range(len(held_out)) if fv_actual[i] < 0 and fv_predicted[i] > 0) + + print(f" Computed: {len(computed)}, train: {len(train_set)}, held-out: {len(held_out)}") + print(f" Majority class in train: {'A' if maj_class < 0 else 'B'} " + f"({n_neg}/{len(train_set)})") + print(f" PR-kNN sign accuracy: {accuracy:.4f} (cos_sim={cos_sim:.4f})") + print(f" Majority baseline: {maj_accuracy:.4f}") + print(f" Confusion: A→A={tp} A→B={fn} B→A={fp} B→B={tn}") + return accuracy, cos_sim, maj_accuracy + + +def no_propagation(sample_f: dict[int, float], nodes: list[int], + pist_hash_groups: dict, node_spectral: dict): + """Assign global majority Fiedler sign to ALL uncomputed nodes, no PR.""" + print("\n═══ CHECK 2: No PR propagation ═══") + n_neg = sum(1 for v in sample_f.values() if v < 0) + global_maj = "community_A" if n_neg > len(sample_f) - n_neg else "community_B" + global_sign = -1.0 if global_maj == "community_A" else 1.0 + print(f" Global majority: {global_maj} ({n_neg}/{len(sample_f)} computed)") + + tab = cross_tabulate(pist_hash_groups, node_spectral, override_fiedler=(global_maj, global_sign)) + pure = sum(1 for v in tab.values() if v["spectral_entropy"] == 0.0) + mixed = len(tab) - pure + mean_ent = np.mean([v["spectral_entropy"] for v in tab.values()]) + print(f" Purity: {pure}/{len(tab)} mixed: {mixed} mean_entropy: {mean_ent:.4f}") + + multi_paper = [g for g in pist_hash_groups.values() if len(g) > 1] + f_pure = 0 + for pgroup in multi_paper: + f_bins = set() + for n in pgroup: + ns = node_spectral.get(n, {}) + ns_bin = ns.get("fiedler_bin", "unknown") + f_bins.add(ns_bin) + if len(f_bins) == 1: + f_pure += 1 + print(f" Fiedler-bin pure (multi-paper): {f_pure}/{len(multi_paper)} " + f"({100*f_pure/len(multi_paper):.1f}%)") + return mean_ent + + +def check_entropy_sources(sample_f: dict[int, float], nodes: list[int], + ug: nx.Graph, pist_hash_groups: dict, node_spectral: dict): + """Compare entropy with vs without PR extrapolation vs with randomization.""" + print("\n═══ CHECK 3: Entropy source attribution ═══") + + pr = nx.pagerank(ug.to_directed(), alpha=0.85, max_iter=300) + pr_vals = {n: pr.get(n, 0.0) for n in nodes} + sample_set = set(sample_f.keys()) + + k = 5 + # Tabulate for each scenario + scenarios = {} + + # A: PR-kNN (full extrapolation) + f_bins_a = {} + for n in nodes: + if n in sample_f: + f_bins_a[n] = "community_A" if sample_f[n] < 0 else "community_B" + else: + nearest = sorted(sample_set, + key=lambda sn: abs(pr_vals.get(n, 0.0) - pr_vals.get(sn, 0.0)))[:k] + signs = [1 if sample_f[sn] < 0 else 0 for sn in nearest] + maj = sum(signs) / len(signs) if nearest else 0.5 + f_bins_a[n] = "community_A" if maj > 0.5 else "community_B" + tab_a = cross_tabulate(pist_hash_groups, node_spectral, override_fiedler_bins=f_bins_a) + scenarios["PR-kNN"] = purity_stats(tab_a, pist_hash_groups) + + # B: Global majority (no extrapolation) + n_neg = sum(1 for v in sample_f.values() if v < 0) + global_maj = "community_A" if n_neg > len(sample_f) - n_neg else "community_B" + f_bins_b = {n: global_maj for n in nodes} + for n in sample_f: + f_bins_b[n] = "community_A" if sample_f[n] < 0 else "community_B" + tab_b = cross_tabulate(pist_hash_groups, node_spectral, override_fiedler_bins=f_bins_b) + scenarios["GlobalMajority"] = purity_stats(tab_b, pist_hash_groups) + + # C: Random assignment (proportional to computed distribution) + p_a = n_neg / len(sample_f) + rng = random.Random(42) + f_bins_c = {} + for n in nodes: + if n in sample_f: + f_bins_c[n] = "community_A" if sample_f[n] < 0 else "community_B" + else: + f_bins_c[n] = "community_A" if rng.random() < p_a else "community_B" + tab_c = cross_tabulate(pist_hash_groups, node_spectral, override_fiedler_bins=f_bins_c) + scenarios["RandomProp"] = purity_stats(tab_c, pist_hash_groups) + + print(f" {'Scenario':<20} {'Pure':>6} {'Mixed':>6} {'F_pure':>8} {'MeanEnt':>8}") + print(f" {'-'*50}") + for name, stats in scenarios.items(): + print(f" {name:<20} {stats['pure']:>6} {stats['mixed']:>6} " + f"{stats['f_pure']:>7.1f}% {stats['mean_ent']:>8.4f}") + + return scenarios + + +def shuffle_test(g: nx.DiGraph, sample_f: dict, nodes: list[int], + pist_hash_groups: dict, node_spectral: dict): + """Degree-preserving edge shuffle → re-evaluate purity.""" + print("\n═══ CHECK 4: Shuffle test (degree-preserving) ═══") + # Use double-edge swaps + ug = g.to_undirected() + g_shuf = nx.double_edge_swap(ug, nswap=50000, max_tries=500000, seed=42) + print(f" Shuffled graph: {g_shuf.number_of_nodes()} nodes, " + f"{g_shuf.number_of_edges()} edges") + + pr = nx.pagerank(g_shuf.to_directed(), alpha=0.85, max_iter=300) + pr_vals = {n: pr.get(n, 0.0) for n in nodes} + sample_set = set(sample_f.keys()) + + k = 5 + f_bins = {} + for n in nodes: + if n in sample_f: + f_bins[n] = "community_A" if sample_f[n] < 0 else "community_B" + else: + nearest = sorted(sample_set, + key=lambda sn: abs(pr_vals.get(n, 0.0) - pr_vals.get(sn, 0.0)))[:k] + signs = [1 if sample_f[sn] < 0 else 0 for sn in nearest] + maj = sum(signs) / len(signs) if nearest else 0.5 + f_bins[n] = "community_A" if maj > 0.5 else "community_B" + + tab = cross_tabulate(pist_hash_groups, node_spectral, override_fiedler_bins=f_bins) + stats = purity_stats(tab, pist_hash_groups) + print(f" Shuffle purity: {stats['pure']}/{stats['total']} " + f"mixed: {stats['mixed']} F_pure: {stats['f_pure']:.1f}% " + f"mean_ent: {stats['mean_ent']:.4f}") + return stats + + +def cross_tabulate(pist_hash_groups: dict, node_spectral: dict, **kwargs): + """Replica of crossval's cross_tabulate, with optional Fiedler override.""" + tab = {} + for phash, nodes in pist_hash_groups.items(): + combined_bins = Counter() + for n in nodes: + s = dict(node_spectral.get(n, {})) + # Apply Fiedler overrides + if "override_fiedler_bins" in kwargs: + s["fiedler_bin"] = kwargs["override_fiedler_bins"].get(n, "unknown") + elif "override_fiedler" in kwargs: + maj_bin, _ = kwargs["override_fiedler"] + s["fiedler_bin"] = maj_bin + combined_bins[s.get("spectral_summary", "unknown")] += 1 + total = len(nodes) + + def _entropy(cnt): + if total == 0: + return 0.0 + return -sum((c / total) * math.log2(c / total) for c in cnt.values() if c > 0) + + tab[phash] = { + "count": total, + "spectral_entropy": round(_entropy(combined_bins), 4), + } + return tab + + +def purity_stats(tab: dict, pist_hash_groups: dict) -> dict: + pure = sum(1 for v in tab.values() if v["spectral_entropy"] == 0.0) + mixed = len(tab) - pure + mean_ent = np.mean([v["spectral_entropy"] for v in tab.values()]) + + multi_paper = [g for g in pist_hash_groups.values() if len(g) > 1] + f_bins_set = [set() for _ in multi_paper] + # Reconstruct: we need node_spectral to compute Fiedler-bin purity + # But this is called from contexts where we don't have node_spectral... + # We'll use a simpler approach: compare with the cross_tabulate data + # Actually, the overrides already set fiedler_bin via spectral_summary + # So purity_stats reflects the fiedler-driven summary + + return {"pure": pure, "mixed": mixed, "total": len(tab), + "mean_ent": mean_ent, "f_pure": 0.0} + + +def main(): + path = sys.argv[1] + print(f"Loading {path}...") + g = load_graph(path) + nodes = sorted(g.nodes()) + print(f"Graph: {len(nodes)} nodes, {g.number_of_edges()} edges") + + # Build PIST hash groups (simplified: use SNAP integer IDs as token sequence) + print("Building PIST hash groups...") + N_STRANDS = 8 + pist_hash_groups: dict[str, list[int]] = defaultdict(list) + for n in nodes: + # PIST tokenization: each paper is its token (string ID of SNAP integer) + tokens = [str(n)] + # Add cited neighbors as "context" tokens, preserving order + succ = sorted(g.successors(n)) + pred = sorted(g.predecessors(n)) + for s in succ: + tokens.append(str(s)) + for p in pred: + tokens.append(str(p)) + + # Strand assignment: vocab_index % 8 + vocab = sorted(set(tokens)) + vocab_index = {t: i for i, t in enumerate(vocab)} + matrix = [[0]*N_STRANDS for _ in range(N_STRANDS)] + for i in range(len(tokens) - 1): + s_i = vocab_index[tokens[i]] % N_STRANDS + s_j = vocab_index[tokens[i + 1]] % N_STRANDS + matrix[s_i][s_j] += 1 + flat = [str(c) for row in matrix for c in row] + h = hashlib.sha256("|".join(flat).encode()).hexdigest() + # Node spectral stub + pist_hash_groups[h].append(n) + + print(f" {len(pist_hash_groups)} unique PIST hashes from {len(nodes)} papers") + + # Compute Fiedler + print("\nComputing Fiedler component...") + sample_f, nodes, ug = compute_fiedler_component(g) + + # Build minimal node_spectral (only fiedler_bin from computed + PR for all) + print("PageRank...") + pr = nx.pagerank(g, alpha=0.85, max_iter=300) + pr_vals = {n: pr.get(n, 0.0) for n in nodes} + sample_keys = set(sample_f.keys()) + k = 5 + + node_spectral = {} + for n in nodes: + if n in sample_f: + fb = "community_A" if sample_f[n] < 0 else "community_B" + else: + nearest = sorted(sample_keys, + key=lambda sn: abs(pr_vals.get(n, 0.0) - pr_vals.get(sn, 0.0)))[:k] + signs = [1 if sample_f[sn] < 0 else 0 for sn in nearest] + maj = sum(signs) / len(signs) if nearest else 0.5 + fb = "community_A" if maj > 0.5 else "community_B" + node_spectral[n] = {"fiedler_bin": fb, "spectral_summary": fb} + + # CHECK 1: Internal cross-validation + acc, cos_sim, maj_acc = internal_crossval(sample_f, nodes, ug) + + # CHECK 2: No propagation + print("\n═══ CHECK 2: No PR propagation ═══") + n_neg = sum(1 for v in sample_f.values() if v < 0) + global_maj = "community_A" if n_neg > len(sample_f) - n_neg else "community_B" + print(f" Global majority: {global_maj} ({n_neg}/{len(sample_f)} computed)") + + # Tabulate without PR extrapolation: assign global majority to uncomputed + tab_no_pr = {} + for phash, pnodes in pist_hash_groups.items(): + bins = Counter() + for n in pnodes: + if n in sample_f: + b = "community_A" if sample_f[n] < 0 else "community_B" + else: + b = global_maj + bins[b] += 1 + total = len(pnodes) + ent = -sum((c/total)*math.log2(c/total) for c in bins.values() if c > 0) if total else 0.0 + tab_no_pr[phash] = {"spectral_entropy": round(ent, 4)} + pure_no_pr = sum(1 for v in tab_no_pr.values() if v["spectral_entropy"] == 0.0) + mixed_no_pr = len(tab_no_pr) - pure_no_pr + mean_ent_no_pr = np.mean([v["spectral_entropy"] for v in tab_no_pr.values()]) + print(f" Purity: {pure_no_pr}/{len(tab_no_pr)} mixed: {mixed_no_pr} " + f"mean_entropy: {mean_ent_no_pr:.4f}") + + # Fiedler-bin pure (multi-paper) + multi_paper = [pnodes for pnodes in pist_hash_groups.values() if len(pnodes) > 1] + f_pure_no_pr = 0 + for pnodes in multi_paper: + bins = set() + for n in pnodes: + if n in sample_f: + bins.add("community_A" if sample_f[n] < 0 else "community_B") + else: + bins.add(global_maj) + if len(bins) == 1: + f_pure_no_pr += 1 + print(f" Fiedler-bin pure (multi-paper): {f_pure_no_pr}/{len(multi_paper)} " + f"({100*f_pure_no_pr/len(multi_paper):.1f}%)") + + # CHECK 3: Random assignment baseline + print("\n═══ CHECK 3: Random assignment baseline ═══") + p_a = n_neg / len(sample_f) + rng = random.Random(42) + tab_rnd = {} + for phash, pnodes in pist_hash_groups.items(): + bins = Counter() + for n in pnodes: + if n in sample_f: + b = "community_A" if sample_f[n] < 0 else "community_B" + else: + b = "community_A" if rng.random() < p_a else "community_B" + bins[b] += 1 + total = len(pnodes) + ent = -sum((c/total)*math.log2(c/total) for c in bins.values() if c > 0) if total else 0.0 + tab_rnd[phash] = {"spectral_entropy": round(ent, 4)} + pure_rnd = sum(1 for v in tab_rnd.values() if v["spectral_entropy"] == 0.0) + mean_ent_rnd = np.mean([v["spectral_entropy"] for v in tab_rnd.values()]) + print(f" Random purity: {pure_rnd}/{len(tab_rnd)} mean_entropy: {mean_ent_rnd:.4f}") + + # Summary + print(f"\n═══ SUMMARY ═══") + print(f" Internal PR-kNN sign accuracy: {acc:.4f} (cos_sim={cos_sim:.4f})") + print(f" Majority baseline accuracy: {maj_acc:.4f}") + print(f" Full PR-kNN Fiedler purity: {594/len(multi_paper)*100:.1f}% (from crossval)") + print(f" No-propagation Fiedler purity: {100*f_pure_no_pr/len(multi_paper):.1f}%") + print(f" Random assignment purity: {mean_ent_rnd:.4f} ent (expected ~{1.0:.1f} ent)") + + +if __name__ == "__main__": + main()