ARCHITECTURE — SpectraFit-Core¶

Overview¶

SpectraFit-Core is a high-performance numerical fitting framework. The computation kernel is written in Rust and exposed to Python via pyo3/maturin. The Python layer provides Pydantic schemas, a DAG composition interface, and an HTML dashboard for result visualisation.

Goals¶

1. Replace lmfit's runtime Python dispatch with compiled Rust model kernels.
2. Replace binary-tree operator-overloaded composition with an explicit DAG IR.
3. Replace per-iteration asteval expression evaluation with pre-compiled expression trees evaluated in Rust.
4. Provide analytical Jacobians for built-in model types where possible (most do; a handful still fall back to finite differences — see the Model::jacobian note below).
5. Support global fits with correct DOF via rayon parallel residual evaluation.

Stack¶

Directory Layout¶

A high-level build/data-flow view — the crate-level dependency DAG lives in Rust crate overview and isn't duplicated here.

flowchart LR
    crates["crates/<br/><small>Rust workspace, 11 crates</small>"] --> python["python/<br/><small>pyo3 binding + oracles</small>"]
    python --> web["web/<br/><small>dashboard + report.html</small>"]
    tests["tests/"] -. validates .-> crates
    tests -. validates .-> python
    benchmark["benchmark/"] --> web
    benchmark --> perf["docs/performance/"]
    docs["docs/"] -. consumes .-> crates
    docs -. consumes .-> python
    docs -. consumes .-> web

The listing below is generated by rrt tree and then filtered by scripts/render_architecture_tree.py to drop every top-level directory listed in scripts/publish_exclusions.py's PRIVATE_ROOTS (working artifacts and per-machine tool state that never reach the public mirror, or are empty scratch dirs with no architectural content). A raw, unfiltered dump would be technically accurate on GitLab but actively misleading here, since the GitHub mirror doesn't have those directories on disk at all. Refresh after a structural change with:

rrt tree --root . --max-depth 2 --snapshot
uv run python scripts/render_architecture_tree.py

rrt tree --root . --max-depth 2 --check still detects drift against the real repo structure (enforced in CI, see ci.yml's lint job) — it compares against the full, unfiltered snapshot in .rrt/tree.lock.toml, independent of this page's filtering.

|-- benchmark/
|   `-- scenarios/
|-- crates/
|   |-- spectrafit-builder/
|   |-- spectrafit-core/
|   |-- spectrafit-dogleg/
|   |-- spectrafit-graph/
|   |-- spectrafit-levenberg-marquardt/
|   |-- spectrafit-models/
|   |-- spectrafit-newton-cg/
|   |-- spectrafit-solver/
|   |-- spectrafit-trust-region/
|   |-- spectrafit-types/
|   |-- spectrafit-varpro/
|   `-- README.md
|-- docs/
|   |-- contributor-guide/
|   |-- explanation/
|   |-- fonts/
|   |-- getting-started/
|   |-- how-to/
|   |-- images/
|   |-- includes/
|   |-- javascripts/
|   |-- reference/
|   |-- release-notes/
|   |-- stylesheets/
|   |-- tutorials/
|   |-- _render_abbr_map.py
|   |-- _render_benchmark_summary.py
|   |-- _render_model_formulas.py
|   |-- _render_nist_tables.py
|   |-- _render_references.py
|   |-- _render_tags.py
|   |-- glossary.md
|   |-- index.md
|   |-- limitations.md
|   |-- references.ris
|   |-- security.md
|   |-- support.md
|   `-- why-spectrafit-core.md
|-- overrides/
|   |-- partials/
|   |-- 404.html
|   |-- home.html
|   |-- main.html
|   `-- section-index.html
|-- python/
|   |-- oracles/
|   `-- spectrafit_core/
|-- reproducibility/
|   |-- figures/
|   |-- ladder/
|   |-- nist_unimplemented/
|   |-- seed-sweep/
|   |-- spectra/
|   |-- assets.toml
|   |-- checksums.sha256
|   `-- ro-crate-metadata.json
|-- scripts/
|   |-- audit/
|   |-- audit_bindings.py
|   |-- audit_latex.py
|   |-- backport_from_github.py
|   |-- bench_background.py
|   |-- bench_ladder.py
|   |-- bg.sh
|   |-- binding_audit_notes.toml
|   |-- check_crate_conventions.py
|   |-- check_pytest_bg.sh
|   |-- check_stale_github_branches.py
|   |-- check_vendored_assets.py
|   |-- coverage_atlas.py
|   |-- fair.py
|   |-- fast_lane_gate.py
|   |-- generate_rustdoc_index.py
|   |-- guard_public_push.py
|   |-- mcp_spectrafit_reports.py
|   |-- measure_significant_digits.py
|   |-- publish_exclusions.py
|   |-- publish_remove_excluded.py
|   |-- publish_snapshot.sh
|   |-- publish_sync.py
|   |-- purge_github_actions_runs.py
|   |-- remote_bench.sh
|   |-- remote_bench_worker.sh
|   |-- render_architecture_tree.py
|   |-- run_pytest_bg.sh
|   |-- self_heal_automation.py
|   |-- shellcheck_ci.py
|   `-- vendored_assets.toml
|-- tests/
|   |-- audit/
|   |-- inference/
|   |-- integration/
|   |-- meta/
|   |-- parity/
|   |-- scenario/
|   |-- unit/
|   `-- conftest.py
|-- web/
|   |-- public/
|   |-- scripts/
|   |-- src/
|   |-- tests/
|   |-- biome.jsonc
|   |-- index.html
|   |-- openapi.snapshot.json
|   |-- package-lock.json
|   |-- package.json
|   |-- playwright.config.ts
|   |-- README.md
|   |-- tsconfig.json
|   |-- vite.config.ts
|   `-- vitest.config.ts
|-- ARCHITECTURE.md
|-- Cargo.lock
|-- Cargo.toml
|-- CHANGELOG.md
|-- CITATION.cff
|-- CLAUDE.md
|-- CODE_OF_CONDUCT.md
|-- codemeta.json
|-- CODEOWNERS
|-- CONTRIBUTING.md
|-- LICENSE
|-- LIMITATIONS.md
|-- package-lock.json
|-- package.json
|-- pyproject.toml
|-- README.md
|-- SECURITY.md
|-- tsconfig.json
|-- uv.lock
`-- zensical.toml

Data Flow¶

sequenceDiagram
    participant Py as Python
    participant Rs as Rust

    Note over Py: FitGraph (Pydantic) + MeasurementData + FitOptions
    Py->>Rs: .model_dump_json() -> JSON strings (pyo3 call)
    Note over Rs: serde_json::from_str -> FitGraphSpec / MeasurementSpec / FitOptionsSpec
    Note over Rs: CompiledGraph::compile(&spec) -- topological sort,<br/>param binding (free_mask: vary=true AND no expr)
    Note over Rs: LmProblem { compiled: &CompiledGraph,<br/>datasets: &[MeasurementSpec], params: DVector, ... }
    Note over Rs: solver dispatch (lm/trf/geodesic/dogleg/newton-cg/varpro)<br/>iterate: residuals(), jacobian() (or FD when expr_edges are present)
    Note over Rs: cov = (JᵀJ)⁻¹ · (chi2 / DOF)<br/>chi2, reduced_chi2, DOF, AIC, BIC
    Rs->>Py: FitResultSpec -> serde_json::to_string() (pyo3 return)
    Note over Py: FitResult.model_validate_json(result_json)

Model Composition — DAG IR¶

Models are defined as a directed acyclic graph at the Python level, serialised to JSON, and evaluated entirely in Rust. Concretely — two peaks tied to a shared width, summed with a background:

flowchart LR
    n1["peak1: Gaussian<br/><small>amplitude, center, sigma</small>"]
    n2["peak2: Gaussian<br/><small>amplitude, center, sigma</small>"]
    n3["bg: Constant<br/><small>c</small>"]
    sum(("Σ"))
    out["model output y(x)"]

    n2 -. "ExprEdge: peak2.sigma = peak1.sigma" .-> n1
    n1 --> sum
    n2 --> sum
    n3 --> sum
    sum --> out

Nodes¶

Each node is a ModelNodeSpec: a typed model instance with a unique id, one ModelType kernel (Gaussian, Lorentzian, Voigt, …), and its parameters keyed by name. An optional dataset_index scopes the node to one dataset in a multi-dataset ("global analysis") fit — None (the default) makes it a global node contributing to every dataset's points; i restricts it to dataset i's residuals and Jacobian columns.

Edges¶

Edges encode parameter constraints (ties) across nodes, and are evaluated in Rust — expr_edges are parsed into an Expr/TiedPlan AST (spectrafit-graph::expr) at compile time, then re-applied every solver iteration by LmProblem::set_free_and_tied (crates/spectrafit-solver/src/lm_problem.rs, shared by both the nalgebra-LM and faer trust-region front-ends) so each tied target is recomputed from its expression before the model is evaluated:

Each edge is an ExprEdge: a target_node / target_param pair naming the parameter to constrain, plus an expression string that references other nodes' parameters in node_id.param form (e.g. "0.5 * peak1.amplitude").

Aggregation¶

Default: sum of all node outputs at each x point.

Why not operator overloading (lmfit-style)?¶

lmfit's model1 + model2 creates a binary tree evaluated recursively at Python speed, allocating N temporary NumPy arrays per iteration. Our DAG is compiled once to a Rust struct; evaluation is a single O(N_nodes * N_x) loop with no Python round-trips.

Parameter Model¶

class Parameter(BaseModel):
    value: float                # initial value
    min: float = -inf
    max: float = inf
    vary: bool = True           # False → fixed constant; ignored when expr is set
    expr: str | None = None     # constraint expression; evaluated every solver iteration
    scale: float | None = None  # solver step-size hint; None → 1.0 (identity, no-op)

name is the dict key in ModelNodeSpec.parameters — not duplicated as a field.

vary is ignored whenever expr is set — the engine always derives the value from the expression and excludes the parameter from the free set regardless of vary's value. There is no validator requiring vary=True when expr is set; vary simply has no effect in that case.

Three binding kinds resolved at compile time (free_mask = vary AND expr is None, spectrafit-graph::compiler):

Bounds (min, max) are enforced by reflective projection, not clamping: a step that overshoots a bound is mirrored back into range (p < lo → 2*lo - p, and symmetrically at hi), with an extreme overshoot parked at the violated bound instead of reflecting past the opposite one. This runs in LmProblem::apply_free_params (crates/spectrafit-solver/src/lm_problem.rs), called from both solver front-ends' set_params, not inside residuals(). scale is applied as an internal change-of-variables preconditioning: the solver works on theta' = theta / scale (LmProblem::scales, apply_free_params/scale_columns_rowmajor in crates/spectrafit-solver/src/lm_problem.rs), not forwarded to any external x_scale field — the levenberg-marquardt crate this workspace vendors has no such field.

Rust Model Kernels¶

pub trait Model: Send + Sync {
    /// x is a coordinate slice: len=1 for 1-D models, len>=n_dims() for nD models.
    fn eval(&self, x: &[f64], params: &[f64]) -> f64;
    /// Default: forward-difference finite differences. Most built-in kernels
    /// override with an analytical formula; 11 currently don't (asym_ir,
    /// breit_wigner, harmonic_ir, kww, log_normal, moffat, pearson7,
    /// split_gaussian, split_pearson7, students_t, tauc — each overrides
    /// with a hand-written central-difference FD loop instead, tracked as a
    /// follow-up to give them analytical Jacobians too).
    fn jacobian(&self, x: &[f64], params: &[f64]) -> Vec<f64> { /* FD fallback */ }
    /// Owned Vec<Cow<'static, str>> — not a static slice, so runtime-generated
    /// models (e.g. GaussianND{d}'s indexed center_0..center_{d-1}) can name
    /// their own params without a compile-time-static list.
    fn param_names(&self) -> Vec<std::borrow::Cow<'static, str>>;
    fn n_dims(&self) -> usize { 1 }  // override for nD models
}

Built-in models¶

The catalog has grown well past an initial handful — 37 wire variants as of this writing. See the Model Reference for the full, authoritative formula table (canonical parameter names, real LaTeX formulas, one row per model, kept in sync with the Rust kernels and the Python parity oracles by convention) — not duplicated here: an earlier hand-copied subset of that table lived on this page, in plain-text notation, and had already drifted (still A * exp(...) after the model reference itself moved to real LaTeX). One table, one place, linked from here instead.

Every model's mixing/weight parameter is named for what it is — e.g. the pseudo-Voigt Lorentzian fraction is always fraction, never eta/frac (see the Model Reference's naming-history section).

DAG Graph Engine (Rust)¶

struct CompiledGraph {
    nodes:          Vec<NodeEntry>,
    free_keys:      Vec<String>,               // "node_id.param_name", sorted
    node_free_cols: Vec<Vec<(usize, usize)>>,  // per-node (local_param_idx, jac_col)
    tied_plan:      TiedPlan,                  // dependency-ordered expr_edge plan
    dataset_offsets: Vec<usize>,                // per-dataset point boundaries (global fits)
}

struct NodeEntry {
    id:            String,
    model:         Box<dyn Model>,
    param_names:   Vec<String>,
    free_mask:     Vec<bool>,           // vary=true AND no expr, per param
    dataset_index: Option<usize>,       // None = global node; Some(i) = dataset-local
}

(ParamBinding/CompiledNode names some earlier drafts of this document used do not exist in the current compiler — see crates/spectrafit-graph/src/compiler.rs.)

evaluate(&graph, x_vals, free_params)¶

output[i] = 0.0
for each node n in compiled.nodes:
    p = resolve params via n.free_mask / n.param_names
    for each i, coord in x_vals:               // coord = &x_vals[i]  (len D)
        output[i] += n.model.eval(coord, &p)

jacobian(&graph, x_vals, free_params)¶

jac[i][j] = 0.0
for each node n in compiled.nodes:
    p = resolve params
    node_jac = n.model.jacobian(&x_vals[i], &p)  // analytical (or FD fallback), coord slice
    for each (local_idx, jac_col) in node_free_cols[node_idx]:
        jac[i][jac_col] += node_jac[local_idx]

When tied_plan is non-empty, ties are re-applied before each evaluation and the analytical Jacobian is swapped for a finite-difference one so the tied-parameter chain-rule terms are captured correctly.

Solver¶

LmProblem<'a> (crates/spectrafit-solver/src/lm_problem.rs — a crate-private module, so this type is not importable from outside spectrafit-solver) is one of ten solver strategies spectrafit-solver dispatches to (Solver in crates/spectrafit-solver/src/dispatch.rs): lm/lm-legacy/trf/geodesic/ dogleg/newton-cg/irls/global/varpro/auto — LM is the default and auto picks VarPro vs. the LM family from the graph shape. How it minimizes weighted residuals, and the chi2/DOF/AIC/BIC/covariance statistics it reports, are covered in Solver — not duplicated here, for the same reason the model-formula table above lives in one place. See Multi-Dataset & Multi-Dimensional Fitting below for the shared-parameter DOF variant.

Standalone Evaluation (No Fitting)¶

The framework is usable without the solver — a first-class path in both Python and Rust.

Python¶

# Evaluate compiled graph at given parameters, no fitting
y_model: np.ndarray = graph.eval(params, data)

# Per-node component outputs
components: dict[str, np.ndarray] = graph.eval_components(params, data)

FitGraph.eval() serialises graph + params + data to JSON, calls the Rust evaluate pyfunction, and returns a numpy array. No LM solver is invoked.

Rust public API¶

The Python-independent, typed API lives in spectrafit-graph/src/lib.rs — pub fn taking real Rust values (not JSON strings), usable as an rlib dependency without pyo3:

pub fn evaluate(
    graph: &FitGraphSpec,
    params_flat: &HashMap<String, f64>,
    x: &[f64],
) -> Result<Vec<f64>, CoreError>;

pub fn evaluate_components(
    graph: &FitGraphSpec,
    params_flat: &HashMap<String, f64>,
    x: &[f64],
) -> Result<HashMap<String, Vec<f64>>, CoreError>;

pub fn jacobian(
    graph: &FitGraphSpec,
    params_flat: &HashMap<String, f64>,
    x: &[f64],
) -> Result<DMatrix<f64>, CoreError>;

Separately, crates/spectrafit-core/src/lib.rs (the pyo3 cdylib) has its own private, JSON-string-based evaluate/evaluate_components functions decorated #[pyfunction] — these wrap the typed API above for the Python FFI boundary and are not themselves pub/rlib-consumable.

Multi-Dataset & Multi-Dimensional Fitting¶

Multi-dimensional independent variables¶

The Model trait accepts x: &[f64] — a coordinate slice of length D. For standard 1-D models D = 1; callers pass &[x_i]. nD models override n_dims() and expect e.g. &[x_i, t_i] for a 2-D energy/temperature surface.

On the Python side MeasurementData.x is a numpy array of shape (N,) for 1-D or (N, D) for nD. Serialised to JSON as a flat list[list[float]] (each inner list is one coordinate vector); 1-D is [[x0], [x1], ...].

class MeasurementData(BaseModel):
    schema_version: str = "0.1"
    x: list[list[float]] | list[float]  # (N, D) matrix, OR flat (N,) — a
                                         # before-validator promotes flat
                                         # input to (N, 1); D=1 for 1-D data
    y: list[float]              # shape (N,)
    sigma: list[float] | None   # shape (N,) or None → uniform weight
    label: str | None = None    # optional dataset identifier

Multi-dataset global fitting¶

fit() accepts a single dataset or a list of datasets. All datasets share the same FitGraph (and therefore the same free parameters). Residuals from all datasets are concatenated before the LM solver sees them.

# single dataset
result = fit(graph, data, options)

# global fit over multiple datasets
result = fit(graph, [data_A, data_B, data_C], options)

On the Rust side LmProblem borrows &[MeasurementSpec] (concatenated into x_concat/y_concat once at construction). Rayon parallelism (spectrafit-graph::executor) is point-wise, not a per-dataset split — a size-based auto-switch between sequential and par_iter_mut/par_chunks_exact inside evaluate_compiled_indexed and the residual/Jacobian kernels, keyed off rayon::current_num_threads() and problem size:

evaluate_compiled_indexed(compiled, node_param_bufs, x_concat, out):
    if x_concat.len() is large enough (rayon::current_num_threads()-aware):
        out.par_iter_mut().zip(x_concat.par_chunks_exact(stride)) ...      // rayon
    else:
        out.iter_mut().zip(x_concat.chunks_exact(stride)) ...              // sequential

DOF = sum_d(N_d) - N_free_shared

The FitResult is the same schema regardless of whether one or many datasets were used; best_fit and residuals are the concatenated arrays in dataset order.

FitResult¶

class ParameterResult(BaseModel):
    value: float
    min: float
    max: float
    vary: bool
    expr: str | None
    scale: float | None         # inherited from Parameter; solver step-size hint
    name: str | None = None     # dotted "node_id.param" name, if known
    stderr: float | None        # None if fit did not converge

class DatasetSlice(BaseModel):
    label: str | None           # from MeasurementData.label
    n_points: int
    best_fit: list[float]       # model output for this dataset
    residuals: list[float]      # (y - f) for this dataset, unweighted
    chi2: float                 # partial chi2 contribution

class FitResult(BaseModel):
    schema_version: str = "0.1"
    parameters: dict[str, ParameterResult]
    covariance: list[list[float | None]] | None
    chi2: float
    reduced_chi2: float
    r_squared: float
    dof: int
    aic: float
    bic: float
    n_iter: int
    n_func_evals: int | None    # None if unavailable
    n_jac_evals: int | None     # None if unavailable
    success: bool
    message: str
    best_fit: list[float]       # concatenated across all datasets
    residuals: list[float]      # concatenated across all datasets
    init_fit: list[float]       # model at initial parameter values
    components: dict[str, list[float]]   # node_id → best_fit per DAG node
    dataset_slices: list[DatasetSlice] | None  # None for single-dataset
    condition_number: float | None       # cond(J^T J) at the solution
    n_de_generations: int | None         # differential-evolution generations, if used
    cost_history: list[float]            # per-iteration cost ½‖r‖² trajectory
    gradient_norm_history: list[float]   # per-iteration ‖J^T r‖_∞
    params_history: list[list[float]]    # per-iteration parameter vector
    covariance_param_order: list[str] | None  # free-param names indexing `covariance`

Python / Rust Boundary¶

All data crosses the boundary as JSON strings, with exactly one allowlisted exemption (fit_arrays_numpy, which returns (json, ndarray) so the fitted curve on the zero-copy path never round-trips through JSON — see CLAUDE.md §4):

• Python → Rust: pydantic_obj.model_dump_json() passed as &str via pyo3
• Rust → Python: serde_json::to_string(&result) returned as String

This keeps the pyo3 FFI layer trivial (no custom type conversions) and makes the boundary independently testable with plain string I/O.

Future: replace JSON with MessagePack for large-dataset performance.

Versioning¶

FitGraph/GlobalFitGraph/FitResult/ExprEdge carry schema_version: str = "0.1". Breaking schema changes bump the minor version. Unknown fields raise a validation error (ConfigDict(extra="forbid")) — not silently ignored.

Out of Scope¶

(ExprEdge Rust evaluation and the benchmark/verification engine — both formerly listed here as out of scope — are implemented; see "Edges" above and "Benchmark Engine + Report" below.)

Benchmark Engine + Report (python/oracles)¶

The benchmark/verification suite is registry-driven and pydantic-first, emitting a frozen JSON contract consumed by a Vite + React UI. There is no Jinja2/HTML artifact — one data flow: benchmark run → results.json → FastAPI → React.

Module structure¶

python/oracles/           # was python/benchmark/, then python/extras/bench/ — merged F13
├── cases.py       # CategoryDef registry + CaseSpec/CaseFamily + build_catalog/materialize
├── models.py      # MODEL_REGISTRY — numpy formulas (the parity oracle for the Rust kernels)
├── backends/      # spectrafit (the SUBJECT) + lmfit, jax/optimistix (cross-check oracles)
├── engine.py      # build_report: run_suite (all 160) + run_featured (deep-dive every case);
│                  #   re-imports _multidim/_global_fit from _engine_multidim (see below)
├── _engine_base.py, _engine_multidim.py, _engine_nested.py, _engine_profile.py
│                  #   engine.py facade split (G27) — _engine_multidim.py's _multidim() fits
│                  #   a genuine ≥3-D gaussian_nd recovery (real subject, not a 2-D map);
│                  #   _global_fit() runs a GlobalFitGraph shared-model multi-spectrum joint fit
├── metrics.py     # timing / accuracy / ECDF / spread / pull statistics
├── synth.py       # deterministic synthetic BenchReport (test fixture; never served)
├── bench_contract.py  # the FROZEN BenchReport contract (Pydantic) — single source of truth
├── contract.py    # small shared-leaf module (SolverMeta, _Contract base) — NOT the BenchReport home
├── api.py         # FastAPI: GET /api/report (latest), /api/runs, /api/report/{run_id}
├── reports.py     # run-centric output: .spectrafit_reports/<category>/<date>_run_NNN/
└── cli.py         # `run` (write results.json + manifest.json) · `gate` (regression gate)

Contract → UI¶

bench_contract.py (not contract.py — that's a small shared-leaf module) defines the frozen BenchReport contract and is the single source of truth. The FastAPI app publishes its OpenAPI schema from it; web/src/openapi.gen.ts is generated from that live schema (npm run contract) and web/src/contract/index.ts re-exports the view types — so the React views never drift from the Python models. The web app (web/, Vite + React) fetches /api/report at boot and renders 2 destinations — Standing (#standing, default: facts masthead + per-backend results table) and Evidence (#evidence: all cases side by side, including the N-D and global-fit "Native showcases" section; #audit redirects here) — with no silent fallback and no hardcoded backend ids (enforced by a vitest source-scan test). A vitest suite (web/src/__tests__/*) renders every panel from fixtures without a browser.

Run & gate¶

uv run poe benchmark writes results.json + manifest.json into a fresh run dir; uv run poe serve serves the latest over FastAPI; python -m oracles.cli gate fails if the geomean speedup vs the pinned baseline_solver_id (default lmfit) drops below 1× or max |Δr²| (LM-family cases) exceeds 1e-3. The test suite (tests/unit/benchmark/, tests/audit/) proves the JSON is real (every category deep-dived, analyzed set multiple + unique, per-case plots distinct, all-finite) before the UI consumes it.

Next steps¶

• Learn the two extension points — Extending SpectraFit-Core lists every touchpoint adding a model or a new solver family requires, building on the crate layout and DAG IR covered above.
• See how this whole workspace is gated — CI pipeline & cache architecture covers why the GitLab/GitHub pipelines are built the way they are.