section11_self_audit.tex

\documentclass[11pt]{article}
\usepackage{amsmath,amssymb,amsthm}
\usepackage{geometry}
\usepackage{booktabs}
\usepackage{graphicx}
\usepackage{hyperref}
\usepackage{enumerate}
\usepackage{listings}
\usepackage{xcolor}

\geometry{letterpaper, margin=1in}

\lstset{
  basicstyle=\ttfamily\small,
  keywordstyle=\color{blue},
  commentstyle=\color{gray},
  stringstyle=\color{red},
  breaklines=true,
  columns=fullflexible,
  keepspaces=true,
  language=Python
}

\title{\textbf{Formation Engine Methodology} \\
\Large Section 11 — Deterministic Pipeline and Self-Audit}
\author{Formation Engine Development Team}
\date{Version 0.1 — As Implemented in Codebase}

\begin{document}
\maketitle

\section{Introduction: Reproducibility as First Principle}

Most computational chemistry frameworks treat reproducibility as a nice-to-have feature. This framework treats it as a \textbf{non-negotiable requirement}. Every formation produced by this engine must be bit-identical across runs, platforms, and library versions when given the same inputs.

This section documents the \textbf{self-audit infrastructure} that enforces this requirement:

\begin{enumerate}
\item \textbf{Deterministic execution} (seed control, IEEE 754 compliance)
\item \textbf{Failure classification} (autonomous categorization of 10,000+ runs)
\item \textbf{Coverage-driven exploration} (gap targeter filling parameter space)
\item \textbf{Regression detection} (compare parameter changes, flag invariant breaks)
\end{enumerate}

Everything documented here \textbf{exists in the codebase} with Python implementations, not theoretical proposals.

\section{Reproducibility Requirements}

\subsection{The Determinism Contract}

Given:
\begin{equation}
(\text{formula}, \; \text{seed}, \; \text{parameters})
\end{equation}

the engine must produce \textbf{bit-identical output}. Not "approximately the same." Not "statistically equivalent." Bit-identical.

\subsection{The Four Pillars}

\paragraph{1. Seed-Controlled RNG}

Every stochastic operation uses \texttt{std::mt19937} with explicit seed:

\begin{lstlisting}[language=C++]
std::mt19937 rng(seed);
std::normal_distribution<double> dist(0.0, sigma);
for (uint32_t i = 0; i < N; ++i) {
    v[i] = dist(rng);
}
\end{lstlisting}

The same seed always produces the same random sequence.

\paragraph{2. Index-Ordered Force Evaluation}

Force accumulation must not depend on thread execution order:

\begin{lstlisting}[language=C++]
// WRONG: race condition on F[i]
#pragma omp parallel for
for (int i = 0; i < N; ++i) {
    for (int j = i+1; j < N; ++j) {
        Vec3 f = compute_force(i, j);
        F[i] += f;   // RACE
        F[j] -= f;   // RACE
    }
}

// RIGHT: per-thread accumulation
#pragma omp parallel
{
    std::vector<Vec3> F_local(N, {0,0,0});
    #pragma omp for
    for (int i = 0; i < N; ++i) {
        for (int j = i+1; j < N; ++j) {
            Vec3 f = compute_force(i, j);
            F_local[i] += f;
            F_local[j] -= f;
        }
    }
    #pragma omp critical
    for (int i = 0; i < N; ++i) {
        F[i] = F[i] + F_local[i];
    }
}
\end{lstlisting}

\paragraph{3. Deterministic Floating-Point}

\begin{itemize}
\item \textbf{No fast-math flags:} \texttt{-ffast-math} allows reordering of floating-point operations, breaking reproducibility
\item \textbf{IEEE 754 compliance:} \texttt{-fno-finite-math-only}, \texttt{-fno-associative-math}
\item \textbf{Explicit rounding modes:} Use default (round-to-nearest-even)
\item \textbf{Avoid platform-specific intrinsics:} \texttt{rsqrt()} may differ between CPUs
\end{itemize}

\paragraph{4. Fixed Library Versions}

The engine pins exact versions:
\begin{itemize}
\item Eigen 3.4.0 (not "latest")
\item ImGui 1.89.2 (not "master")
\item C++17 standard (not C++20, which changes \texttt{constexpr} rules)
\end{itemize}

Version drift across runs violates determinism.

\subsection{Validation Test}

\textbf{Test protocol:}
\begin{enumerate}
\item Run H$_2$O formation with seed 42, parameters $\{T=300\text{K}, \rho=1.0\text{g/cm}^3\}$
\item Record final geometry (64-bit binary dump of positions)
\item Repeat 10 times on different machines (Windows, Linux, macOS)
\item Compute SHA-256 hash of each geometry
\end{enumerate}

\textbf{Pass criterion:} All 10 hashes are identical.

\textbf{Current status:} Pass on same platform, \textit{not yet validated} across OS (floating-point library differences).

\section{Failure Classifier}

\textbf{Implementation:} \texttt{tools/failure\_classifier.py}, 250 lines

The failure classifier categorizes every failed simulation into one of six buckets using regex pattern matching on error messages.

\subsection{The Six Categories}

From \texttt{failure\_classifier.py}, lines 14--19:

\begin{lstlisting}[language=Python]
class FailureCategory(Enum):
    NUMERICAL = "numerical"      # dt too big, overflow, NaN
    PHYSICS = "physics"          # clash explosion, unphysical bonds
    OOD = "out_of_domain"        # composition/scale outside training
    CONVERGENCE = "convergence"  # optimizer failed to converge
    TIMEOUT = "timeout"          # exceeded time limit
    UNKNOWN = "unknown"          # unclassified
\end{lstlisting}

\subsection{Pattern Matching Rules}

Lines 34--58:

\begin{lstlisting}[language=Python]
PATTERNS = {
    FailureCategory.NUMERICAL: [
        r"(?i)nan|inf|overflow|underflow",
        r"(?i)timestep.*too.*large",
        r"(?i)numerical.*unstable",
        r"(?i)division.*by.*zero",
        r"(?i)sqrt.*negative",
    ],
    FailureCategory.PHYSICS: [
        r"(?i)clash|explosion|overlap",
        r"(?i)unphysical.*bond",
        r"(?i)energy.*diverged",
        r"(?i)force.*too.*large",
        r"(?i)invalid.*geometry",
    ],
    FailureCategory.OOD: [
        r"(?i)unknown.*element",
        r"(?i)unsupported.*composition",
        r"(?i)atom.*count.*exceeded",
    ],
    FailureCategory.CONVERGENCE: [
        r"(?i)failed.*converge",
        r"(?i)max.*iterations",
        r"(?i)optimizer.*failed",
    ],
    FailureCategory.TIMEOUT: [
        r"(?i)timeout|time.*limit|killed",
    ],
}
\end{lstlisting}

\subsection{Classification Algorithm}

Lines 71--85:

\begin{lstlisting}[language=Python]
def classify(self, error_message: str, stderr: str, 
             returncode: int) -> FailureCategory:
    """Classify failure based on error patterns."""
    
    # Timeout is special (returncode)
    if returncode in [-9, 124, 137]:  # SIGKILL, timeout
        return FailureCategory.TIMEOUT
    
    # Combine error sources
    full_error = f"{error_message}\n{stderr}"
    
    # Match patterns
    for category, patterns in self.PATTERNS.items():
        for pattern in patterns:
            if re.search(pattern, full_error):
                return category
    
    return FailureCategory.UNKNOWN
\end{lstlisting}

\subsection{Minimal Reproduction Command}

Every failure stores a minimal reproduction command (lines 102--112):

\begin{lstlisting}[language=Python]
def generate_minimal_repro(self, seed: int, formula: str, 
                          params: Dict) -> str:
    """Generate minimal reproduction command."""
    
    cmd_parts = [
        f"./build/meso-sim",
        f"--formula {formula}",
        f"--seed {seed}",
    ]
    
    for key, val in params.items():
        cmd_parts.append(f"--{key} {val}")
    
    return " ".join(cmd_parts)
\end{lstlisting}

\paragraph{Example output:}
\begin{verbatim}
./build/meso-sim --formula NaCl --seed 12345 --temp 5000 --dt 10
\end{verbatim}

This allows instant reproduction of any failure without searching through logs.

\subsection{The UNKNOWN Budget}

The UNKNOWN category must remain below 5\% of total failures. If it exceeds this threshold, the pattern matchers need updating.

\textbf{Rationale:} A large UNKNOWN fraction means the classifier is missing common failure modes, reducing its diagnostic value.

\subsection{Output Format}

The classifier produces a JSON report:

\begin{lstlisting}[language=Python]
{
  "total_runs": 10000,
  "total_failures": 1247,
  "failures_by_category": {
    "numerical": 342,
    "physics": 198,
    "out_of_domain": 67,
    "convergence": 521,
    "timeout": 89,
    "unknown": 30
  },
  "unknown_fraction": 0.024,
  "failures": [
    {
      "category": "numerical",
      "reason": "NaN in force evaluation at step 127",
      "seed": 12345,
      "formula": "H2O",
      "params": {"temp": 300, "dt": 0.001},
      "minimal_repro": "./build/meso-sim --formula H2O --seed 12345 ..."
    },
    ...
  ]
}
\end{lstlisting}

\section{Gap Targeter (Coverage-Driven Exploration)}

\textbf{Implementation:} \texttt{tools/gap\_targeter.py}, 200 lines

The gap targeter partitions parameter space into a 3D grid and schedules runs to fill sparse regions.

\subsection{The Parameter Grid}

From \texttt{gap\_targeter.py}, lines 38--40:

\begin{lstlisting}[language=Python]
# Define grid
self.scale_bins = [2, 5, 10, 20, 50, 100]  # atom count
self.temp_bins = np.linspace(50, 500, 10)   # K
self.density_bins = np.linspace(0.001, 0.1, 10)  # g/cm^3
\end{lstlisting}

\begin{itemize}
\item \textbf{Scale:} 6 bins (2, 5, 10, 20, 50, 100 atoms)
\item \textbf{Temperature:} 10 bins (50 K to 500 K)
\item \textbf{Density:} 10 bins (0.001 to 0.1 g/cm$^3$)
\end{itemize}

Total: $6 \times 10 \times 10 = 600$ cells.

\subsection{Cell Representation}

Lines 10--30:

\begin{lstlisting}[language=Python]
@dataclass
class Cell:
    """Cell in parameter space grid."""
    scale: int
    temp: float
    density: float
    
    # Coverage metrics
    n_runs: int = 0
    n_success: int = 0
    n_failures: int = 0
    
    # Mismatch metrics
    mean_energy: float = 0.0
    std_energy: float = 0.0
    max_mismatch: float = 0.0  # vs neighbors
    
    def is_sparse(self, threshold: int = 10) -> bool:
        return self.n_runs < threshold
    
    def is_high_mismatch(self, threshold: float = 0.5) -> bool:
        return self.max_mismatch > threshold
\end{lstlisting}

\subsection{Gap Identification}

A cell is considered a \textbf{gap} if:

\paragraph{Condition 1: Sparse coverage}
\begin{equation}
n_{\text{runs}} < 10
\end{equation}

\paragraph{Condition 2: High mismatch with neighbors}
\begin{equation}
\max_{j \in \text{neighbors}} \frac{|\bar{E}_i - \bar{E}_j|}{\sqrt{\sigma_i^2 + \sigma_j^2}} > 0.5
\end{equation}

The mismatch metric detects discontinuities in the energy landscape (phase transitions, force field breakdowns).

\subsection{Mismatch Computation}

Lines 79--100:

\begin{lstlisting}[language=Python]
def compute_mismatches(self):
    """Compute mismatch between neighboring cells."""
    
    for (i, j, k), cell in self.cells.items():
        if cell.n_success < 2:
            continue
        
        # Check 6 neighbors (±1 in each dimension)
        neighbors = [
            (i-1, j, k), (i+1, j, k),
            (i, j-1, k), (i, j+1, k),
            (i, j, k-1), (i, j, k+1),
        ]
        
        mismatches = []
        for neighbor_key in neighbors:
            if neighbor_key in self.cells:
                neighbor = self.cells[neighbor_key]
                if neighbor.n_success >= 2:
                    # Energy mismatch
                    delta = abs(cell.mean_energy - neighbor.mean_energy)
                    combined_std = np.sqrt(cell.std_energy**2 
                                         + neighbor.std_energy**2)
                    if combined_std > 0:
                        mismatch = delta / combined_std
                        mismatches.append(mismatch)
        
        if mismatches:
            cell.max_mismatch = max(mismatches)
\end{lstlisting}

\subsection{Run Scheduling}

Lines 102--125:

\begin{lstlisting}[language=Python]
def schedule_runs(self, n_runs: int) -> List[Dict]:
    """Generate run schedule to fill gaps."""
    
    gaps = self.identify_gaps()
    
    if not gaps:
        print("No gaps detected - coverage is complete!")
        return []
    
    print(f"Identified {len(gaps)} gap cells")
    print(f"Scheduling {n_runs} runs to improve coverage...")
    
    schedule = []
    for i in range(n_runs):
        # Round-robin through gaps
        cell = gaps[i % len(gaps)]
        
        run_config = {
            "scale": cell.scale,
            "temp": cell.temp,
            "density": cell.density,
            "formula": self._generate_formula(cell.scale),
            "seed": i + 1000000,  # Unique seed
        }
        
        schedule.append(run_config)
    
    return schedule
\end{lstlisting}

\subsection{Convergence Criterion}

Coverage is considered \textbf{complete} when:
\begin{equation}
\frac{n_{\text{cells with } n_{\text{runs}} \geq 10}}{600} > 0.95
\end{equation}

95\% of cells have at least 10 successful runs.

\subsection{Usage Example}

\begin{lstlisting}[language=bash]
# Initialize gap targeter
python tools/gap_targeter.py --output results/coverage.json

# Run 1000 iterations
for i in {1..1000}; do
    # Get next batch of runs
    python tools/gap_targeter.py --schedule 10 > batch.json
    
    # Execute batch
    while read run; do
        eval $run
    done < batch.json
    
    # Update coverage
    python tools/gap_targeter.py --update results/run_$i.json
done

# Generate coverage report
python tools/gap_targeter.py --report > coverage_report.txt
\end{lstlisting}

\section{Regression Detector}

\textbf{Implementation:} \texttt{tools/regression\_detector.py}, 200 lines

The regression detector reruns a benchmark suite when parameters change and reports:
\begin{itemize}
\item Which molecules changed classification (stable $\to$ unstable)
\item Score deltas exceeding threshold
\item Invariant violations
\end{itemize}

\subsection{Benchmark Suite Format}

Each line is a JSON object:

\begin{lstlisting}[language=Python]
{"formula": "H2O", "seed": 42, "expected_class": "stable"}
{"formula": "CH4", "seed": 43, "expected_class": "stable"}
{"formula": "NaCl", "seed": 44, "expected_class": "ionic"}
...
\end{lstlisting}

\subsection{Configuration Hashing}

Every configuration gets a deterministic hash (lines 48--51):

\begin{lstlisting}[language=Python]
def hash_config(self, config: Dict) -> str:
    """Hash config for change detection."""
    config_str = json.dumps(config, sort_keys=True)
    return hashlib.sha256(config_str.encode()).hexdigest()[:12]
\end{lstlisting}

\paragraph{Example:}
\begin{lstlisting}[language=Python]
config = {"bond_cutoff": 1.2, "angle_threshold": 0.1}
hash = "a3f4b7c8d1e2"
\end{lstlisting}

\subsection{The Four Invariants}

Lines 41--46:

\begin{lstlisting}[language=Python]
INVARIANTS = {
    "energy_monotonic": "Energy should decrease during relaxation",
    "bonds_stable": "Bond count should stabilize after equilibration",
    "score_bounded": "Score should be in [0, 1]",
    "classification_consistent": "Same input -> same classification",
}
\end{lstlisting}

\subsection{Comparison Algorithm}

Lines 109--144:

\begin{lstlisting}[language=Python]
def compare_runs(self, baseline: List[Result], 
                current: List[Result]) -> ChangeReport:
    """Compare two benchmark runs."""
    
    # Index by (formula, seed)
    baseline_map = {(r.formula, r.seed): r for r in baseline}
    current_map = {(r.formula, r.seed): r for r in current}
    
    # Find differences
    reclassified = []
    score_deltas = {}
    
    for key, curr in current_map.items():
        if key not in baseline_map:
            continue
        
        base = baseline_map[key]
        
        # Classification changed?
        if base.classification != curr.classification:
            reclassified.append({
                "formula": curr.formula,
                "seed": curr.seed,
                "old_class": base.classification,
                "new_class": curr.classification,
                "score_delta": curr.score - base.score,
            })
        
        # Score changed significantly?
        delta = abs(curr.score - base.score)
        if delta > 0.01:  # 1% threshold
            score_deltas[curr.formula] = curr.score - base.score
    
    # Check invariants
    invariant_breaks = self._check_invariants(current)
    
    return ChangeReport(
        config_hash_old="baseline",
        config_hash_new="current",
        config_diff="",
        total_runs=len(current),
        reclassified=len(reclassified),
        score_changed=len(score_deltas),
        invariant_breaks=invariant_breaks,
        reclassified_details=reclassified,
        score_deltas=score_deltas
    )
\end{lstlisting}

\subsection{Invariant Checking}

Lines 146--157:

\begin{lstlisting}[language=Python]
def _check_invariants(self, results: List[Result]) -> List[str]:
    """Check if any invariants are broken."""
    
    breaks = []
    
    # score_bounded
    for r in results:
        if not (0.0 <= r.score <= 1.0):
            breaks.append(
                f"score_bounded: {r.formula} has score "
                f"{r.score:.3f} (out of [0,1])"
            )
    
    # TODO: Check other invariants (needs trajectory data)
    
    return breaks
\end{lstlisting}

\subsection{Output Report}

\begin{lstlisting}[language=Python]
{
  "config_hash_old": "a3f4b7c8d1e2",
  "config_hash_new": "e9d2c1b7f4a3",
  "config_diff": "bond_cutoff: 1.2 -> 1.3",
  "total_runs": 100,
  "reclassified": 3,
  "score_changed": 12,
  "invariant_breaks": [
    "score_bounded: C6H6 has score 1.23 (out of [0,1])"
  ],
  "reclassified_details": [
    {
      "formula": "H2O",
      "seed": 42,
      "old_class": "stable",
      "new_class": "unstable",
      "score_delta": -0.23
    }
  ],
  "score_deltas": {
    "CH4": +0.015,
    "NH3": -0.022,
    ...
  }
}
\end{lstlisting}

\subsection{Binary Verdict}

The detector produces:
\begin{itemize}
\item ✅ \textbf{NO REGRESSIONS} — if reclassified $= 0$ and no invariant breaks
\item ⚠️ \textbf{REVIEW REQUIRED} — otherwise
\end{itemize}

Parameter changes causing regressions are not accepted until:
\begin{enumerate}
\item The regressions are understood (root cause identified)
\item Either fixed or documented as acceptable trade-offs
\end{enumerate}

\section{The Self-Audit Loop}

The three tools work together in a continuous improvement cycle:

\begin{center}
\begin{tabular}{ccc}
\toprule
Tool & Input & Output \\
\midrule
Failure Classifier & Error logs & Failure buckets + minimal repro \\
Gap Targeter & Coverage grid & Sparse/mismatch cells \\
Regression Detector & Config changes & Invariant breaks + verdict \\
\bottomrule
\end{tabular}
\end{center}

\subsection{The Weekly Cycle}

\paragraph{Monday:} Run 10,000 simulations
\begin{itemize}
\item Failure classifier categorizes all failures
\item UNKNOWN category analyzed: new patterns added
\end{itemize}

\paragraph{Tuesday:} Gap analysis
\begin{itemize}
\item Gap targeter identifies sparse cells
\item Schedule 1,000 runs to fill gaps
\end{itemize}

\paragraph{Wednesday:} Execute gap-filling runs
\begin{itemize}
\item Update coverage metrics
\item Check if 95\% threshold reached
\end{itemize}

\paragraph{Thursday:} Parameter tuning
\begin{itemize}
\item Propose scoring weight changes
\item Run regression detector on benchmark suite
\end{itemize}

\paragraph{Friday:} Review regressions
\begin{itemize}
\item If NO REGRESSIONS: accept changes
\item If REVIEW REQUIRED: investigate, fix, or revert
\end{itemize}

\section{Validation and Testing}

\subsection{Determinism Validation}

\textbf{Test 1: Bit-identical runs}

\begin{lstlisting}[language=bash]
# Run 10 times
for i in {1..10}; do
    ./build/meso-sim --formula H2O --seed 42 --output run_$i.xyz
    sha256sum run_$i.xyz >> hashes.txt
done

# Check uniqueness
uniq hashes.txt | wc -l  # Should be 1
\end{lstlisting}

\textbf{Pass criterion:} All 10 hashes identical.

\textbf{Test 2: Platform independence}

Run same command on Windows, Linux, macOS. Compare hashes.

\textbf{Current status:} Pass within platform, \textit{not yet validated} across platforms.

\subsection{Failure Classifier Validation}

\textbf{Test:} Hand-label 100 failures, compare to classifier output.

\textbf{Metric:} Cohen's kappa $\kappa > 0.8$ (strong agreement).

\textbf{Current status:} Not yet performed (requires labeled dataset).

\subsection{Gap Targeter Validation}

\textbf{Test:} After 1000 iterations, check coverage:
\begin{equation}
\frac{n_{\text{cells with } n_{\text{runs}} \geq 10}}{600} > 0.95
\end{equation}

\textbf{Current status:} Framework implemented, not yet executed at scale.

\subsection{Regression Detector Validation}

\textbf{Test:} Inject known regression (break energy monotonicity), verify detector flags it.

\textbf{Current status:} Manual testing only.

\section{Known Limitations and Future Work}

\subsection{Determinism: Cross-Platform}

\begin{itemize}
\item Windows and Linux floating-point libraries differ (glibc vs. MSVCRT)
\item \texttt{sin()}, \texttt{cos()}, \texttt{exp()} may vary in last bits
\item Solution: Use fixed-precision math library (e.g., MPFR) or accept $\varepsilon = 10^{-12}$ tolerance
\end{itemize}

\subsection{Failure Classifier: Pattern Coverage}

\begin{itemize}
\item Current patterns based on manual inspection of $\sim$100 failures
\item Needs validation on larger corpus (10,000+ failures)
\item UNKNOWN fraction currently not monitored in production
\end{itemize}

\subsection{Gap Targeter: Formula Generation}

\begin{itemize}
\item Current formula selection is heuristic (lines 144--153)
\item Does not adapt to chemistry (e.g., metals vs. organics)
\item Better approach: Use composition library with verified structures
\end{itemize}

\subsection{Regression Detector: Trajectory Analysis}

\begin{itemize}
\item Energy monotonicity check requires trajectory data
\item Current implementation only checks final results
\item Extension: Store intermediate frames, check invariants at each step
\end{itemize}

\section{Conclusion: The Self-Correcting Engine}

This section has documented the self-audit infrastructure \textbf{as it exists in the codebase}:

\begin{itemize}
\item \textbf{Failure classifier:} 6 categories, regex pattern matching, minimal repro generation (250 lines Python)
\item \textbf{Gap targeter:} 3D parameter grid (600 cells), mismatch detection, coverage-driven scheduling (200 lines)
\item \textbf{Regression detector:} Config hashing, 4 invariants, binary verdict (200 lines)
\end{itemize}

The framework operates in a continuous loop:
\begin{enumerate}
\item Run simulations at scale (10,000+ runs)
\item Classify failures autonomously
\item Identify parameter space gaps
\item Schedule targeted runs to fill gaps
\item Detect regressions when parameters change
\item Accept changes only if no invariants break
\end{enumerate}

This is \textbf{not a theoretical proposal}. The tools exist, have clear interfaces, and integrate into the formation pipeline. What remains is \textit{execution at scale}---running the full loop for weeks and validating the 95\% coverage criterion.

The principle: \textbf{trust but verify}. The formation engine does not assume its results are correct. It continuously audits itself, flags anomalies, and demands human review when invariants break. This is the difference between a research prototype and a production scientific instrument.

\paragraph{Transition to \S12:} With the self-audit infrastructure established (failure classification, gap targeting, regression detection), Section~12 documents the \textit{validation doctrine}: the hierarchical test suite that certifies the formation engine for production use.

\end{document}