#!/usr/bin/env python3
"""
compute-hsi.py -- HSI (Hybrid Similarity Index) Computation Pipeline
=====================================================================
Ported from V4 hsi_semantic_surprise.py and V2 compute_lsa.py.
Reads room/*.md artifacts, computes TF-IDF/SVD structural similarity
and embedding semantic similarity, outputs .hsi-results.json with scored pairs.

v1.6.0 "Powerhouse" upgrade: Spectral OM-HMM scoring based on Markov chain
theory (Seabrook & Wiskott 2022, arxiv 2207.02296). Replaces keyword-density
proxy with transition-matrix spectral analysis of thinking modes. The spectral
gap of the thinking-mode Markov chain measures genuine integrative thinking
quality -- fast mixing across modes = rich cross-domain thinking.

Usage:
    python3 scripts/compute-hsi.py /path/to/room [--tier 1|2] [--threshold 0.30] [--output path]
                                                 [--scope-to-nodes]

Tier system (per HSI-05):
    Tier 0: Handled by analyze-room bash script (keyword-only) -- NOT this script
    Tier 1: sklearn + MiniLM (default, CPU-only, ~80MB model)
    Tier 2: sklearn + Pinecone (uses existing embeddings if configured)
"""

# Reference implementation, the MINDRIAN_RS_BACKEND=python fallback and
# scout's HSI step 1. Its surprise_type strings use the retired Convention B;
# since Phase 355 D-07 no reader trusts them: hsi-to-graph re-derives from
# lsa_sim / semantic_sim through lib/core/direction-convention.cjs.
import argparse
import hashlib
import json
import math
import os
import re
import sqlite3
import sys
from datetime import datetime, timezone
from pathlib import Path

# Auto-install Python ML deps if missing (v1.10.9, plan 85-10, LAWRENCE-001).
# compute-hsi.py is the direct entry point for /mos:act and /mos:mos-reason's
# HSI recomputation cascade; it never had this wiring even though the sibling
# whitespace scripts did (RCA phase-134-python-elimination-false-complete --
# Change 1). sentence-transformers is included even though its own import is
# soft-guarded (Tier 1 embeddings) because an unremediated miss there still
# hard-fails a few lines later once the Tier 2 Pinecone fallback also misses.
sys.path.insert(0, str(Path(__file__).resolve().parent / "lib"))
from ensure_ml_deps import ensure
ensure(["numpy", "scikit-learn", "sentence-transformers"])

# --- Guarded imports ---

try:
    import numpy as np
except ImportError:
    print("HSI requires numpy. Run: pip install -r requirements-hsi.txt", file=sys.stderr)
    sys.exit(1)

try:
    from sklearn.decomposition import TruncatedSVD
    from sklearn.feature_extraction.text import TfidfVectorizer
    from sklearn.metrics.pairwise import cosine_similarity
except ImportError:
    print("HSI requires scikit-learn. Run: pip install -r requirements-hsi.txt", file=sys.stderr)
    sys.exit(1)


# --- Spectral OM-HMM: Thinking Mode Classification & Markov Analysis ---
# Based on Seabrook & Wiskott (2022) "Tutorial on Spectral Theory of Markov Chains"
# arxiv 2207.02296 -- applied to venture artifact integrative thinking detection.
#
# Core insight: Texts that rapidly transition between thinking modes (analytical,
# integrative, descriptive, evaluative, creative) exhibit higher integrative
# thinking quality than texts stuck in a single mode. The spectral gap of the
# thinking-mode transition matrix quantifies this mixing rate.

_INTEGRATIVE_KEYWORDS = frozenset({
    "cross-domain", "synthesis", "combine", "bridge", "transfer",
    "connect", "integrate", "hybrid", "convergence", "interdisciplinary",
    "analogy", "metaphor", "parallel", "intersection", "fusion",
})

# Thinking mode classifiers (sentence-level)
# Each mode represents a hidden state in the Markov chain
_THINKING_MODES = {
    "analytical": re.compile(
        r"\b(because|therefore|consequently|evidence|data|measure|quantif|statistic|analyz|assess|evaluat|compar)\w*\b", re.I
    ),
    "integrative": re.compile(
        r"\b(connect|bridge|synthes|combin|integrat|cross|interdisciplin|convergence|fusion|hybrid|anolog|metaphor|transfer)\w*\b", re.I
    ),
    "descriptive": re.compile(
        r"\b(is|are|was|were|has|have|consist|compris|includ|contain|describ|defin|refer|represent)\w*\b", re.I
    ),
    "evaluative": re.compile(
        r"\b(should|must|better|worse|risk|opportunit|strength|weakness|advantage|disadvantage|critical|important|significant)\w*\b", re.I
    ),
    "creative": re.compile(
        r"\b(novel|innovati|reimagin|redefin|what.if|could|might|envision|transform|disrupt|pioneer|breakthrough|radical)\w*\b", re.I
    ),
}

# Legacy feature patterns (kept for backward compatibility in scoring)
_FEATURE_PATTERNS = [
    (r"\bsimple\b|\bstraightforward\b|\bbasic\b", "simple"),
    (r"\bcomplex\b|\bcomplicated\b|\bintricate\b", "complex"),
    (r"\blinear\b|\bsequential\b|\bstep.by.step\b", "linear"),
    (r"\bmulti\w*\b|\bparallel\b|\bsimultaneous\b", "multidirectional"),
    (r"\bpart\b|\bcomponent\b|\bpiece\b|\bfragment\b", "part"),
    (r"\bholistic\b|\bwhole\b|\bsystem\b|\bentire\b", "holistic"),
    (r"\btrade.?off\b|\bcompromise\b|\bbalance\b", "tradeoff"),
    (r"\bcreative\b|\bnovel\b|\binnovati\\w+\b|\boriginal\b", "creative"),
]

# --- Shared corpus exclude-list (Phase 200-01 / SEED-018, extended 233-03) ---
#
# SKIP_DIRS / SKIP_FILES / MIN_BODY_CHARS come from the ONE canonical source,
# lib/core/rs_corpus_exclude.py, whose own docstring states the contract: every
# room-artifact walker imports those three names from that module and keeps no
# local literal.
#
# compute-hsi.py VIOLATED that contract until Phase 233-03. It was the one walker
# Phase 200-01 never migrated (rs-engine.py, rs_hybrid.py and rs_rooms.py were),
# and its local copy drifted exactly the way the shared source exists to prevent:
# it never learned .snapshots or sub-rooms, so on motj it scored 207 artifacts
# whose top-20 pairs were ALL .snapshots/sub-rooms near-duplicate state dumps,
# and hsi-to-graph.cjs then wrote ZERO edges because none of those ids are
# Artifact nodes in the parent room.db. See the RCA at
# .planning/debug/graph-derive-silent-clear-dead-api-derivation.md Section 9
# Defect #4. Keep NO local literal here.
#
# (The retired local set also carried ".hsi-cache.json", which was already dead:
# os.walk's `dirs` list holds directory names only, never file names.)
_THIS_FILE = Path(__file__).resolve()
_LIB_CORE = _THIS_FILE.parent.parent / 'lib' / 'core'
if str(_LIB_CORE) not in sys.path:
    sys.path.insert(0, str(_LIB_CORE))

from rs_corpus_exclude import SKIP_DIRS, SKIP_FILES, MIN_BODY_CHARS  # noqa: E402


def parse_frontmatter(content):
    """Extract frontmatter fields using regex (no PyYAML dependency)."""
    fm_match = re.match(r'^---\n([\s\S]*?)\n---', content)
    if not fm_match:
        return {}
    fm_text = fm_match.group(1)
    fields = {}
    for line in fm_text.split('\n'):
        if ':' in line:
            key, _, val = line.partition(':')
            fields[key.strip()] = val.strip().strip('"').strip("'")
    return fields


def extract_title(content, filepath):
    """Extract title from first # heading."""
    match = re.search(r'^# (.+)$', content, re.MULTILINE)
    if match:
        return match.group(1).strip()
    return Path(filepath).stem.replace('-', ' ').title()


def extract_body(content):
    """Extract body text after frontmatter --- block."""
    fm_match = re.match(r'^---\n[\s\S]*?\n---\n?', content)
    if fm_match:
        return content[fm_match.end():]
    return content


def discover_artifacts(room_dir):
    """Walk room_dir for .md files, build artifact list."""
    artifacts = []
    room_path = Path(room_dir).resolve()

    for root, dirs, files in os.walk(room_path):
        # Filter out skip dirs
        dirs[:] = [d for d in dirs if d not in SKIP_DIRS]

        rel_root = Path(root).relative_to(room_path)
        # Skip root-level files (no section)
        if str(rel_root) == '.':
            continue

        section = str(rel_root).split(os.sep)[0]

        for fname in sorted(files):
            if not fname.endswith('.md'):
                continue
            if fname in SKIP_FILES:
                continue

            fpath = Path(root) / fname
            try:
                content = fpath.read_text(encoding='utf-8')
            except (OSError, UnicodeDecodeError):
                continue

            artifact_id = str(Path(rel_root) / Path(fname).stem).replace(os.sep, '/')
            title = extract_title(content, fpath)
            body = extract_body(content)

            if len(body.strip()) < MIN_BODY_CHARS:
                continue  # skip near-empty artifacts

            artifacts.append({
                'id': artifact_id,
                'section': section,
                'title': title,
                'path': str(fpath.relative_to(room_path)),
                'text': body.strip(),
            })

    return artifacts


def _file_uri_path(p):
    """Percent-encode the URI-significant bytes of a path for a SQLite file: URI.

    SQLite's URI parser treats ? as the query separator, # as the fragment start,
    and % as the escape byte. A room path containing any of them, spliced raw into
    'file:%s?mode=ro', silently breaks the open: the parser swallows part of the
    path as a query string, and the mode=ro that follows is then never recognized
    as the read-only mode flag. That is a broken read-only GUARANTEE (T-233-09
    claims read-only is enforced mechanically at the SQLite layer, not by
    convention), not merely a broken path.

    This is byte-for-byte the rule lib/core/graph-derivation.cjs::_fileUriPath
    already applies to the SAME room.db path for its sub-room ATTACH: the same
    three bytes, encoded in the same order (% first, so it cannot double-encode
    the escapes it just introduced). Kept identical on purpose. Everything else
    (spaces, quotes, unicode) passes through untouched, because SQLite decodes
    %XX and takes the rest literally. Two encoders that disagree about what a
    room path means would be a worse bug than either one alone.
    """
    return str(p).replace('%', '%25').replace('?', '%3F').replace('#', '%23')


def load_graph_artifact_ids(room_dir):
    """Return the set of Artifact node ids in <room>/.mindrian/room.db, or None.

    Phase 233-03 (RCA Section 9 Defect #4/#5). hsi-to-graph.cjs already refuses to
    write an edge whose endpoint is not an Artifact NODE (its findArtifact guard,
    correct behavior). Without this intersection, compute-hsi.py spends a full
    LSA + embedding pass on artifacts that could never become an edge, and the
    top-20 cut is then consumed by pairs that get silently discarded downstream.
    Scoping to the node set is both cheaper AND more correct.

    READ-ONLY by construction (T-233-09): the connection is opened through the
    `file:...?mode=ro` URI, so a write attempt fails at the SQLite layer rather
    than by convention. This script must never mutate the room's authoritative db.
    The path is run through _file_uri_path first, so a room directory containing
    ?, # or % cannot break the URI parse and take the mode=ro flag down with it.

    Returns None (never raises, never returns an empty set on error) when the db
    is missing, unreadable, or has no `nodes` table yet -- the Tier 0 /
    pre-structural-index case. The caller then scores everything, exactly like
    today. Degrading to "no filter" is right; degrading to "zero artifacts" would
    turn an immature room into a silent no-op, which is the false-success class
    this whole phase exists to close.
    """
    db_path = Path(room_dir) / '.mindrian' / 'room.db'
    if not db_path.exists():
        print(
            "HSI: --scope-to-nodes: no room.db at %s, scoring every discovered "
            "artifact" % db_path,
            file=sys.stderr,
        )
        return None

    conn = None
    try:
        conn = sqlite3.connect('file:%s?mode=ro' % _file_uri_path(db_path), uri=True)
        rows = conn.execute(
            "SELECT id FROM nodes WHERE type = 'Artifact'"
        ).fetchall()
    except sqlite3.Error as exc:
        print(
            "HSI: --scope-to-nodes: room.db unreadable (%s), scoring every "
            "discovered artifact" % exc,
            file=sys.stderr,
        )
        return None
    finally:
        if conn is not None:
            try:
                conn.close()
            except sqlite3.Error:
                pass

    return {r[0] for r in rows}


def compute_content_hashes(artifacts):
    """Compute MD5 hash of each artifact's text content."""
    hashes = {}
    for art in artifacts:
        h = hashlib.md5(art['text'].encode('utf-8')).hexdigest()[:12]
        hashes[art['id']] = h
    return hashes


def check_cache(room_dir, current_hashes):
    """Check if all artifact hashes match cached hashes. Return True if unchanged."""
    cache_path = Path(room_dir) / '.hsi-cache.json'
    if not cache_path.exists():
        return False
    try:
        cached = json.loads(cache_path.read_text(encoding='utf-8'))
        cached_hashes = cached.get('hashes', {})
        if set(cached_hashes.keys()) != set(current_hashes.keys()):
            return False
        return all(cached_hashes.get(k) == v for k, v in current_hashes.items())
    except (json.JSONDecodeError, OSError):
        return False


def write_cache(room_dir, hashes):
    """Write content hashes to cache file."""
    cache_path = Path(room_dir) / '.hsi-cache.json'
    cache_path.write_text(json.dumps({
        'timestamp': datetime.now(timezone.utc).isoformat(),
        'hashes': hashes,
    }, indent=2), encoding='utf-8')


def compute_lsa_similarity(texts):
    """Compute LSA structural similarity matrix using TF-IDF + SVD.

    Ported from V4 compute_lsa_similarity and V2 compute_lsa.py.
    Uses 500 max_features (lighter than V2's 2000 -- room artifacts are smaller).
    """
    n = len(texts)
    vectorizer = TfidfVectorizer(stop_words='english', max_features=500)

    try:
        tfidf_matrix = vectorizer.fit_transform(texts)
    except ValueError:
        # Empty vocabulary
        return np.eye(n)

    n_components = min(80, n - 1, tfidf_matrix.shape[1])
    if n_components < 1:
        return np.eye(n)

    svd = TruncatedSVD(n_components=n_components)
    reduced = svd.fit_transform(tfidf_matrix)

    sim_matrix = cosine_similarity(reduced)
    return np.clip(sim_matrix, 0.0, 1.0)


def compute_semantic_similarity_tier1(texts):
    """Tier 1: Local MiniLM embeddings (CPU-only, ~80MB model)."""
    try:
        from sentence_transformers import SentenceTransformer
    except ImportError:
        return None

    model = SentenceTransformer('all-MiniLM-L6-v2')
    embeddings = model.encode(texts, show_progress_bar=False)
    sim_matrix = cosine_similarity(embeddings)
    return np.clip(sim_matrix, 0.0, 1.0)


def compute_semantic_similarity_tier2(artifact_ids, texts):
    """Tier 2: Pinecone embeddings from existing index.

    Falls back to Tier 1 if Pinecone not configured or query fails.
    """
    api_key = os.environ.get('PINECONE_API_KEY')
    index_name = os.environ.get('PINECONE_INDEX')

    if not api_key or not index_name:
        return None  # Fall back to Tier 1

    try:
        from pinecone import Pinecone
        pc = Pinecone(api_key=api_key)
        index = pc.Index(index_name)

        # Fetch existing embeddings for artifact IDs
        fetch_result = index.fetch(ids=artifact_ids)
        vectors = fetch_result.get('vectors', {})

        if len(vectors) < 2:
            return None  # Not enough vectors, fall back

        # Build embedding matrix in artifact order
        dim = None
        embeddings = []
        for aid in artifact_ids:
            if aid in vectors:
                vec = vectors[aid]['values']
                if dim is None:
                    dim = len(vec)
                embeddings.append(vec)
            else:
                return None  # Missing vector, fall back to Tier 1

        embedding_matrix = np.array(embeddings)
        sim_matrix = cosine_similarity(embedding_matrix)
        return np.clip(sim_matrix, 0.0, 1.0)

    except Exception:
        return None  # Any error, fall back to Tier 1


def classify_sentence_mode(sentence):
    """Classify a sentence into its dominant thinking mode.

    Returns the mode with the highest keyword match count.
    Ties broken by priority: integrative > creative > evaluative > analytical > descriptive.
    Falls back to 'descriptive' if no matches (most common baseline mode).
    """
    scores = {}
    for mode, pattern in _THINKING_MODES.items():
        scores[mode] = len(pattern.findall(sentence))

    max_score = max(scores.values())
    if max_score == 0:
        return "descriptive"

    priority = ["integrative", "creative", "evaluative", "analytical", "descriptive"]
    for mode in priority:
        if scores[mode] == max_score:
            return mode
    return "descriptive"


def build_transition_matrix(mode_sequence):
    """Build a Markov transition matrix from a sequence of thinking modes.

    Returns (matrix, mode_names) where matrix[i][j] = P(mode_j | mode_i).
    Based on Seabrook & Wiskott (2022) Section 2: transition matrices encode
    state-to-state probabilities in a Markov chain.
    """
    modes = list(_THINKING_MODES.keys())
    n = len(modes)
    mode_idx = {m: i for i, m in enumerate(modes)}

    # Count transitions
    counts = np.zeros((n, n), dtype=float)
    for k in range(len(mode_sequence) - 1):
        i = mode_idx[mode_sequence[k]]
        j = mode_idx[mode_sequence[k + 1]]
        counts[i][j] += 1.0

    # Normalize rows to get transition probabilities
    # Add Laplace smoothing (alpha=0.1) to avoid zero rows
    alpha = 0.1
    counts += alpha

    row_sums = counts.sum(axis=1, keepdims=True)
    row_sums[row_sums == 0] = 1.0  # safety
    matrix = counts / row_sums

    return matrix, modes


def compute_spectral_gap(transition_matrix):
    """Compute spectral gap of a Markov chain transition matrix.

    Spectral gap = 1 - |lambda_2| where lambda_2 is the second-largest
    eigenvalue by magnitude. Larger gap = faster mixing = more diverse
    thinking mode transitions.

    From Seabrook & Wiskott (2022) Section 4: the spectral gap determines
    the rate of convergence to the stationary distribution. A chain that
    mixes quickly visits all states rapidly -- in our context, this means
    the text transitions through many thinking modes rather than getting
    stuck in one.
    """
    try:
        eigenvalues = np.linalg.eigvals(transition_matrix)
        # Sort by magnitude descending
        magnitudes = sorted(np.abs(eigenvalues), reverse=True)

        if len(magnitudes) < 2:
            return 0.0

        # lambda_1 should be ~1.0 (Perron-Frobenius)
        # Spectral gap = 1 - |lambda_2|
        spectral_gap = 1.0 - magnitudes[1]
        return max(0.0, min(1.0, spectral_gap))
    except (np.linalg.LinAlgError, ValueError):
        return 0.0


def compute_stationary_distribution(transition_matrix):
    """Compute stationary distribution of the Markov chain.

    The stationary distribution pi satisfies pi * P = pi.
    It reveals the long-run proportion of time spent in each thinking mode.
    From Seabrook & Wiskott (2022) Section 3: for ergodic chains, the
    stationary distribution is unique and equals the left eigenvector
    corresponding to eigenvalue 1.
    """
    try:
        # Left eigenvector: solve pi * P = pi, equivalently P^T * pi^T = pi^T
        eigenvalues, eigenvectors = np.linalg.eig(transition_matrix.T)

        # Find eigenvector for eigenvalue closest to 1
        idx = np.argmin(np.abs(eigenvalues - 1.0))
        stationary = np.real(eigenvectors[:, idx])

        # Normalize to probability distribution
        total = np.sum(np.abs(stationary))
        if total > 0:
            stationary = np.abs(stationary) / total
        else:
            stationary = np.ones(len(stationary)) / len(stationary)

        return stationary
    except (np.linalg.LinAlgError, ValueError):
        n = transition_matrix.shape[0]
        return np.ones(n) / n


def detect_absorbing_tendency(transition_matrix, modes):
    """Detect if the chain has absorbing tendencies (gets stuck in modes).

    A self-loop probability > 0.6 indicates the text tends to stay in
    that mode rather than transitioning -- a sign of shallow or
    single-mode analysis.

    Returns: absorbing_score (0-1) where 0 = no absorption, 1 = fully stuck.
    """
    diag = np.diag(transition_matrix)
    # Average self-loop probability above baseline (1/n would be uniform)
    n = len(modes)
    baseline = 1.0 / n
    excess = np.maximum(diag - baseline, 0)
    absorbing_score = float(np.mean(excess) / (1.0 - baseline)) if baseline < 1.0 else 0.0
    return max(0.0, min(1.0, absorbing_score))


def compute_omhmm_score(text):
    """Spectral OM-HMM integrative thinking score (0-100).

    v1.6.0 upgrade: Uses Markov chain spectral analysis instead of keyword
    density proxy. Based on Seabrook & Wiskott (2022) spectral theory of
    Markov chains applied to thinking-mode transition detection.

    Score components:
      - spectral_gap (40%): Fast mixing across thinking modes = rich thinking
      - integrative_weight (25%): Stationary distribution weight on integrative mode
      - mode_diversity (20%): Entropy of stationary distribution (Ashby's variety)
      - anti_absorption (15%): Penalty for getting stuck in single modes

    Falls back to legacy keyword scoring if text has < 5 sentences
    (insufficient data for meaningful transition matrix).
    """
    # Split into sentences for mode classification
    sentences = re.split(r'[.!?]+', text)
    sentences = [s.strip() for s in sentences if len(s.strip()) > 10]

    # Fallback to legacy scoring for very short texts
    if len(sentences) < 5:
        return _compute_omhmm_legacy(text)

    # Step 1: Classify each sentence into a thinking mode
    mode_sequence = [classify_sentence_mode(s) for s in sentences]

    # Step 2: Build transition matrix (Seabrook & Wiskott Section 2)
    transition_matrix, modes = build_transition_matrix(mode_sequence)

    # Step 3: Compute spectral gap (Section 4)
    # Larger gap = faster mixing = more diverse thinking
    spectral_gap = compute_spectral_gap(transition_matrix)

    # Step 4: Compute stationary distribution (Section 3)
    stationary = compute_stationary_distribution(transition_matrix)

    # Step 5: Extract integrative mode weight from stationary distribution
    mode_idx = {m: i for i, m in enumerate(modes)}
    integrative_weight = float(stationary[mode_idx["integrative"]])

    # Step 6: Compute mode diversity via Shannon entropy
    # High entropy = visits many modes = Ashby's requisite variety
    entropy = 0.0
    for p in stationary:
        if p > 1e-10:
            entropy -= p * math.log2(p)
    max_entropy = math.log2(len(modes))  # Uniform distribution
    mode_diversity = entropy / max_entropy if max_entropy > 0 else 0.0

    # Step 7: Detect absorbing tendency (penalty for shallow thinking)
    absorbing = detect_absorbing_tendency(transition_matrix, modes)
    anti_absorption = 1.0 - absorbing

    # Composite score (0-100)
    score = (
        0.40 * spectral_gap * 100
        + 0.25 * integrative_weight * 100 * len(modes)  # Scale by mode count
        + 0.20 * mode_diversity * 100
        + 0.15 * anti_absorption * 100
    )

    return max(0.0, min(100.0, score))


def _compute_omhmm_legacy(text):
    """Legacy OM-HMM scoring for short texts (< 5 sentences).

    Original V4 keyword-density approach. Used as fallback when there
    is insufficient data to build a meaningful transition matrix.

    Score = 0.5 * likelihood_ratio * 100
          + 0.3 * (feature_diversity / 8) * 100
          + 0.2 * min(complexity_ratio * 50, 100)
    """
    text_lower = text.lower()
    words = text_lower.split()
    total_words = max(len(words), 1)

    # Feature diversity
    features_found = set()
    for pattern, label in _FEATURE_PATTERNS:
        if re.search(pattern, text_lower):
            features_found.add(label)
    feature_diversity = len(features_found)

    # Complexity ratio: sentence length variance / mean
    sentences = re.split(r'[.!?]+', text)
    sentences = [s.strip() for s in sentences if s.strip()]
    if len(sentences) > 1:
        lengths = [len(s.split()) for s in sentences]
        mean_len = sum(lengths) / len(lengths)
        if mean_len > 0:
            variance = sum((ln - mean_len) ** 2 for ln in lengths) / len(lengths)
            complexity_ratio = (variance ** 0.5) / mean_len
        else:
            complexity_ratio = 0.0
    else:
        complexity_ratio = 0.0

    # Likelihood ratio: integrative keyword density
    integrative_count = sum(
        1 for w in words if w.strip('.,;:!?') in _INTEGRATIVE_KEYWORDS
    )
    likelihood_ratio = integrative_count / total_words

    score = (
        0.5 * likelihood_ratio * 100
        + 0.3 * (feature_diversity / 8) * 100
        + 0.2 * min(complexity_ratio * 50, 100)
    )

    return max(0.0, min(100.0, score))


def compute_artifact_spectral_profile(text):
    """Compute full spectral profile for an artifact.

    Returns dict with spectral metrics for downstream graph analysis:
    - omhmm_score: composite integrative thinking score (0-100)
    - spectral_gap: mixing rate of thinking-mode Markov chain (0-1)
    - dominant_mode: mode with highest stationary probability
    - mode_entropy: Shannon entropy of stationary distribution (0-1 normalized)
    - absorbing_score: tendency to get stuck in single mode (0-1)
    - mode_distribution: full stationary distribution dict
    """
    sentences = re.split(r'[.!?]+', text)
    sentences = [s.strip() for s in sentences if len(s.strip()) > 10]

    if len(sentences) < 5:
        # Legacy fallback -- limited profile
        score = _compute_omhmm_legacy(text)
        return {
            'omhmm_score': score,
            'spectral_gap': 0.0,
            'dominant_mode': 'unknown',
            'mode_entropy': 0.0,
            'absorbing_score': 0.0,
            'mode_distribution': {},
            'spectral_method': 'legacy',
        }

    mode_sequence = [classify_sentence_mode(s) for s in sentences]
    transition_matrix, modes = build_transition_matrix(mode_sequence)
    spectral_gap = compute_spectral_gap(transition_matrix)
    stationary = compute_stationary_distribution(transition_matrix)

    mode_idx = {m: i for i, m in enumerate(modes)}
    dominant_idx = int(np.argmax(stationary))
    dominant_mode = modes[dominant_idx]

    entropy = 0.0
    for p in stationary:
        if p > 1e-10:
            entropy -= p * math.log2(p)
    max_entropy = math.log2(len(modes))
    mode_entropy = entropy / max_entropy if max_entropy > 0 else 0.0

    absorbing = detect_absorbing_tendency(transition_matrix, modes)

    mode_dist = {m: round(float(stationary[mode_idx[m]]), 4) for m in modes}

    return {
        'omhmm_score': compute_omhmm_score(text),
        'spectral_gap': round(spectral_gap, 4),
        'dominant_mode': dominant_mode,
        'mode_entropy': round(mode_entropy, 4),
        'absorbing_score': round(absorbing, 4),
        'mode_distribution': mode_dist,
        'spectral_method': 'markov',
    }


def compute_hsi_matrix(artifacts, lsa_matrix, semantic_matrix, threshold=0.30):
    """Compute HSI innovation differential matrix and extract top pairs.

    v1.6.0 upgrade: Uses spectral OM-HMM scores and includes spectral
    metadata (spectral_gap, dominant_mode) in pair output for downstream
    graph analysis and random walk innovation pathway discovery.

    Formula:
    - semantic_surprise = abs(semantic_sim - lsa_sim)
    - integrative_factor = sqrt(omhmm_i * omhmm_j) / 100
    - innovation_differential = 0.6 * semantic_surprise + 0.4 * integrative_factor
    - breakthrough_potential = 0.7 * differential + 0.3 * min(lsa, semantic)
    """
    n = len(artifacts)
    texts = [a['text'] for a in artifacts]

    # Compute spectral profiles for all artifacts
    print("HSI: computing spectral OM-HMM profiles...", file=sys.stderr)
    spectral_profiles = [compute_artifact_spectral_profile(t) for t in texts]
    omhmm_scores = [p['omhmm_score'] for p in spectral_profiles]

    # Count spectral vs legacy
    n_spectral = sum(1 for p in spectral_profiles if p['spectral_method'] == 'markov')
    n_legacy = n - n_spectral
    print(f"HSI: {n_spectral} spectral, {n_legacy} legacy (< 5 sentences)", file=sys.stderr)

    pairs = []
    for i in range(n):
        for j in range(i + 1, n):
            lsa_sim = float(lsa_matrix[i, j])
            sem_sim = float(semantic_matrix[i, j])

            semantic_surprise = abs(sem_sim - lsa_sim)
            integrative_factor = math.sqrt(omhmm_scores[i] * omhmm_scores[j]) / 100.0
            innovation_diff = 0.6 * semantic_surprise + 0.4 * integrative_factor

            if innovation_diff < threshold:
                continue

            # Classify
            if lsa_sim > sem_sim:
                surprise_type = 'structural_transfer'
            else:
                surprise_type = 'semantic_implementation'

            # Breakthrough potential (from V2)
            breakthrough = 0.7 * innovation_diff + 0.3 * min(lsa_sim, sem_sim)

            # Spectral metadata for pair (average of both artifacts)
            avg_spectral_gap = (
                spectral_profiles[i]['spectral_gap']
                + spectral_profiles[j]['spectral_gap']
            ) / 2.0

            pairs.append({
                'left_id': artifacts[i]['id'],
                'right_id': artifacts[j]['id'],
                'lsa_sim': round(lsa_sim, 4),
                'semantic_sim': round(sem_sim, 4),
                'hsi_score': round(innovation_diff, 4),
                'surprise_type': surprise_type,
                'breakthrough_potential': round(breakthrough, 4),
                'spectral_gap_avg': round(avg_spectral_gap, 4),
                'left_dominant_mode': spectral_profiles[i]['dominant_mode'],
                'right_dominant_mode': spectral_profiles[j]['dominant_mode'],
            })

    # Sort by hsi_score descending, take top 20
    pairs.sort(key=lambda p: p['hsi_score'], reverse=True)
    return pairs[:20], spectral_profiles


def main():
    parser = argparse.ArgumentParser(
        description='HSI computation pipeline for room artifacts'
    )
    parser.add_argument('room_dir', help='Path to room directory')
    parser.add_argument('--tier', type=int, choices=[1, 2], default=1,
                        help='Embedding tier: 1=MiniLM (default), 2=Pinecone')
    parser.add_argument('--threshold', type=float, default=0.30,
                        help='Minimum HSI score threshold (default: 0.30)')
    parser.add_argument('--output', default=None,
                        help='Output JSON path (default: {room_dir}/.hsi-results.json)')
    parser.add_argument('--scope-to-nodes', action='store_true',
                        help='Score only artifacts that already exist as Artifact '
                             'nodes in <room>/.mindrian/room.db (read-only). '
                             'Degrades to scoring everything if the db is absent.')

    args = parser.parse_args()
    room_dir = Path(args.room_dir).resolve()

    if not room_dir.is_dir():
        print(f"Error: {room_dir} is not a directory", file=sys.stderr)
        sys.exit(1)

    output_path = args.output or str(room_dir / '.hsi-results.json')

    # Step 1: Discover artifacts
    artifacts = discover_artifacts(room_dir)

    # Step 1b: optional intersection with the graph's own Artifact node set, so
    # scoring covers exactly the pairs hsi-to-graph.cjs can actually write.
    if args.scope_to_nodes:
        node_ids = load_graph_artifact_ids(room_dir)
        if node_ids is not None:
            discovered = len(artifacts)
            artifacts = [a for a in artifacts if a['id'] in node_ids]
            print(
                "HSI: --scope-to-nodes: %d/%d discovered artifacts exist as "
                "Artifact nodes" % (len(artifacts), discovered),
                file=sys.stderr,
            )

    if len(artifacts) < 2:
        # Write empty results and exit
        result = {
            'metadata': {
                'timestamp': datetime.now(timezone.utc).isoformat(),
                'room_dir': str(room_dir),
                'tier': args.tier,
                'artifact_count': len(artifacts),
                'pair_count': 0,
            },
            'artifacts': [],
            'hsi_pairs': [],
            'reverse_salients': [],
        }
        Path(output_path).write_text(json.dumps(result, indent=2), encoding='utf-8')
        print(f"HSI: {len(artifacts)} artifacts found (minimum 2 required), wrote empty results",
              file=sys.stderr)
        sys.exit(0)

    # Step 2: Check cache
    content_hashes = compute_content_hashes(artifacts)
    if check_cache(room_dir, content_hashes):
        print("HSI: all artifacts unchanged (cache hit), skipping computation",
              file=sys.stderr)
        sys.exit(0)

    # Step 3: Compute LSA structural similarity
    texts = [a['text'] for a in artifacts]
    print(f"HSI: computing LSA similarity for {len(artifacts)} artifacts...",
          file=sys.stderr)
    lsa_matrix = compute_lsa_similarity(texts)

    # Step 4: Compute semantic similarity
    tier_used = args.tier
    semantic_matrix = None

    if args.tier == 2:
        artifact_ids = [a['id'] for a in artifacts]
        semantic_matrix = compute_semantic_similarity_tier2(artifact_ids, texts)
        if semantic_matrix is not None:
            tier_used = 2
            print("HSI: using Tier 2 (Pinecone) embeddings", file=sys.stderr)

    if semantic_matrix is None:
        # Tier 1 fallback
        print("HSI: using Tier 1 (MiniLM) embeddings...", file=sys.stderr)
        semantic_matrix = compute_semantic_similarity_tier1(texts)
        tier_used = 1

        if semantic_matrix is None:
            print(
                "HSI: sentence-transformers not installed. "
                "Run: pip install -r requirements-hsi.txt",
                file=sys.stderr
            )
            sys.exit(1)

    # Step 5: Compute HSI matrix and extract top pairs (with spectral profiles)
    pairs, spectral_profiles = compute_hsi_matrix(
        artifacts, lsa_matrix, semantic_matrix, threshold=args.threshold
    )

    # Step 6: Build output with spectral metadata
    artifact_list = [
        {
            'id': a['id'],
            'section': a['section'],
            'title': a['title'],
            'path': a['path'],
            'spectral': {
                'omhmm_score': round(spectral_profiles[i]['omhmm_score'], 2),
                'spectral_gap': spectral_profiles[i]['spectral_gap'],
                'dominant_mode': spectral_profiles[i]['dominant_mode'],
                'mode_entropy': spectral_profiles[i]['mode_entropy'],
                'absorbing_score': spectral_profiles[i]['absorbing_score'],
                'mode_distribution': spectral_profiles[i]['mode_distribution'],
                'method': spectral_profiles[i]['spectral_method'],
            },
        }
        for i, a in enumerate(artifacts)
    ]

    # Compute room-level spectral summary
    spectral_scores = [p['omhmm_score'] for p in spectral_profiles]
    spectral_gaps = [p['spectral_gap'] for p in spectral_profiles if p['spectral_method'] == 'markov']
    mode_counts = {}
    for p in spectral_profiles:
        dm = p['dominant_mode']
        mode_counts[dm] = mode_counts.get(dm, 0) + 1

    result = {
        'metadata': {
            'timestamp': datetime.now(timezone.utc).isoformat(),
            'room_dir': str(room_dir),
            'tier': tier_used,
            'artifact_count': len(artifacts),
            'pair_count': len(pairs),
            'spectral_version': '1.6.0',
            'spectral_summary': {
                'mean_omhmm': round(sum(spectral_scores) / len(spectral_scores), 2) if spectral_scores else 0,
                'mean_spectral_gap': round(sum(spectral_gaps) / len(spectral_gaps), 4) if spectral_gaps else 0,
                'dominant_mode_distribution': mode_counts,
                'spectral_artifacts': len(spectral_gaps),
                'legacy_artifacts': len(spectral_scores) - len(spectral_gaps),
            },
        },
        'artifacts': artifact_list,
        'hsi_pairs': pairs,
        'reverse_salients': [],  # Populated by detect-reverse-salients.py
    }

    Path(output_path).write_text(json.dumps(result, indent=2), encoding='utf-8')

    # Update cache
    write_cache(room_dir, content_hashes)

    # Summary to stderr
    top_score = pairs[0]['hsi_score'] if pairs else 0.0
    print(
        f"HSI: {len(pairs)} pairs found (top score: {top_score:.3f}), "
        f"tier {tier_used}, {len(artifacts)} artifacts",
        file=sys.stderr
    )


if __name__ == '__main__':
    main()