Audit any forward model¶

twindex.general.audit_model searches for different parameters that explain the same observations under bounded noise. Supply a Python callable wrapping your ODE solver, finite element model, PINN, neural surrogate, or other forward code.

This produces numerical evidence, not a proof. The exact PDE checker remains available through twindex audit and twindex check. Numerical reports have a different format and cannot be replayed as certificates.

This runs your code

audit_model calls your Python function, and twindex audit-model module:function imports and executes your module. Both are local developer tools: run them only on code you trust, and never expose them as a public service. The result is numerical evidence, not a proof.

First audit¶

Install the optional numerical dependencies:

pip install 'twindex[general]'

Suppose the observation only measures the sum of two parameters:

import json
import numpy as np
from twindex.general import audit_model

def forward(theta):
    return np.array([theta[0] + theta[1]])

report = audit_model(
    forward,
    data=[1.0],
    prior=[[0, 1], [0, 1]],
    epsilon=0.01,
    estimate=[0.5, 0.5],  # optional; defaults to prior midpoint
    starts=4,
    seed=0,
)
for parameter in report.parameters:
    print(parameter['parameter'], parameter['observed_gap'])
    print(parameter['low'], parameter['high'])
with open('audit.json', 'w') as stream:
    json.dump(report.to_dict(), stream, indent=2, allow_nan=False)

Both parameters have observed gaps close to their full prior width. Each endpoint includes its parameter vector, fresh forward prediction, and maximum residual divided by its noise bound. This demonstrates ambiguity of the supplied numerical map at sampled points. It does not bound the true PDE solution error.

What the report means¶

Field	Interpretation
`evidence: numerical_only`, `certified: false`	No exact theorem has been checked.
`local.reference`	Best least-squares reference found by the starts. It need not be feasible.
`local.scaled_jacobian`	Derivative in prior-normalized free coordinates, with each observation divided by its noise bound.
`local.singular_values`, `rank`, `null_directions_scaled`	Local linear sensitivity diagnostics. Null directions are in prior-normalized coordinates, with zeros in fixed coordinates. The numerical rank threshold is recorded.
`local.linearised_ranges`	LP extrema of the linear approximation intersected with the prior. These are neither nonlinear outer bounds nor certificates.
`parameters[k].low`, `.high`	Evaluated feasible samples with smallest/largest sampled coordinate `k`. Both include all coordinates.
`parameters[k].observed_gap`, `.prior_fraction`	Observed sample separation and its fraction of the prior width. Never a global upper bound. The fraction is `null` for fixed coordinates.
`search.runs`, `.forward_calls`	Solver termination messages and actual forward evaluation count.

A small gap is compatible with a search that missed an entire disconnected component. For example, forward(theta) = [theta[0]**2], data [1], and prior [[-2, 2]] has two branches, although its Jacobian at theta = [1] has full rank. A zero derivative at another point may also be inconclusive for nonlinear identifiability. The report deliberately does not classify parameters as "identifiable", "weakly identifiable", or "not identifiable".

no_feasible_sample_found describes the search, not emptiness of the compatible set. Forward errors or malformed/nonfinite outputs raise an exception. They are never converted into evidence of infeasibility.

Forward-map contract¶

The callable accepts one real NumPy vector of length N and returns a finite real vector of length M. Use the same parameter meaning, observation order, units, mesh, and solver settings on every call. It must be deterministic and defined throughout the prior. The current API supports finite box priors, including fixed coordinates. epsilon is positive, either a scalar or one bound per observation: abs(forward(theta)[i] - data[i]) <= epsilon[i].

These are bounded error sets, not Gaussian confidence regions. Account for measurement and numerical/model error in your interpretation; Twindex does not estimate or certify the latter. A Fisher information interpretation would require an additional likelihood model and is not part of this API.

The audit uses multistart bounded least squares, a local Jacobian and scaled LPs, then multistart SLSQP searches for each coordinate's extrema. Feasibility is checked again against the forward map, without adding a tolerance to epsilon. Consequently, a point just outside the bound because of rounding is rejected. Optimizer success alone never establishes feasibility. The method uses SciPy's local least-squares solver, linear programming, and SLSQP.

Pass jacobian=your_jacobian for an (M, N) derivative matrix in physical parameter units for the local diagnostics. Otherwise, bounded two-sided or one-sided finite differences are used with diff_step in prior-normalized coordinates. Finite difference rank can be sensitive to solver accuracy and step size. The default relative rank tolerance is max(1e-7, diff_step) for finite differences, and max(M, Nfree) * machine_epsilon for supplied derivatives. Pass rank_rtol to override it when solver accuracy requires a larger tolerance. singular_value_ratios records each singular value divided by the largest (all zeros for a zero Jacobian). rank_rtol, rank_threshold and derivative_source make the numerical rank decision explicit. Null directions are numerically rank-deficient directions whose ratio is at most rank_rtol, including structural null directions. The supplied Jacobian is currently used for local diagnostics only.

starts, seed, and max_iter control search effort; max_iter applies to each fit/optimization, not a total evaluation budget. Reported evaluations include fresh endpoint checks. This initial implementation uses dense SVD and LPs and searches every free coordinate; expensive or high-dimensional models can require substantial work. No speed or global-coverage guarantee is made.

ODE and 2-D finite-difference examples¶

From a repository checkout, examples/general_models.py supplies:

decay: initial amplitude and decay rate in u' = -rate*u, evaluated with solve_ivp at five times. It checks forward solver convergence.
poisson_2d: diffusivity and source for -a*Laplacian(u) = s on a square, zero boundary values, using a fixed 9 by 9 interior grid. Scaling both parameters equally produces identical measurements.

from examples.general_models import poisson_2d
from twindex.general import audit_model

report = audit_model(poisson_2d, poisson_2d([1, 1]),
                     [[0.5, 2], [0.5, 2]], 1e-4, estimate=[1, 1])

This is an audit of a finite-difference model, not a certified 2-D PDE backend.

PyTorch and PINN wrappers¶

A surrogate taking parameters as inputs can be wrapped directly:

import torch

model.eval()
def forward(theta):
    # Match device and dtype to the trained model; select observations here.
    parameter = next(model.parameters())
    t = torch.as_tensor(theta, dtype=parameter.dtype, device=parameter.device)
    with torch.no_grad():
        return model(t[None, :])[0].detach().cpu().numpy()

For a PINN that solves one PDE at one parameter, the wrapper must instead solve or retrain at each supplied theta and sample the same observation points. A fixed trained x -> u(x) network cannot stand in for theta -> measurements. For DeepXDE-style workflows, export the parameter-dependent solve and sampling as this callable. Training randomness must be controlled. These integration patterns carry the accuracy of your numerical model; no framework-specific certification is claimed.

CLI¶

From the repository root (or with your own module on Python's import path):

python -c 'import json; json.dump([1.0], open("measurements.json", "w"))'
twindex audit-model examples.general_models:ridge \
    --data measurements.json --prior 0,1 --prior 0,1 \
    --eps 0.01 --estimate 0.5,0.5 --starts 4 --seed 0 --output audit.json

--data is a JSON vector; --prior occurs once per parameter. Use --epsilon 0.01,0.02 for two unequal observation bounds. A negative lower prior bound can be passed as --prior=-2,2. The command imports and executes the named Python function. Exit 0 means feasible numerical samples were found, 2 means none were found or the command failed (errors go to stderr). Neither exit code is an identifiability verdict. Use the JSON status to distinguish a completed search without feasible samples from an error.

Certification is a separate replay¶

The existing 1-D measured-data audit uses the exact model and problem definition to generate and replay certificates. A generic callable cannot identify its own governing equation, parameter map, or discretization error, so Twindex does not silently turn this numerical report into a proof. Use the explicit bridge after constructing an exact 1-D Problem:

from twindex.general import certify_audit

result = certify_audit(report, problem)
for parameter in result.certified:
    print(parameter['status'], parameter['attempts'])

The caller asserts that this exact problem describes the callable's physics, sensor order, units, data, noise bound and prior. The bridge cannot prove that correspondence. It generates exact polynomials with the declared boundary values, then replays pair certificates independently for each coordinate. Each attempt retains its certificate and checker report, or a generation error. Default shrink factors 1, 9/10, 1/2 move both endpoints toward the numerical reference after rationalisation; the reference itself need not be compatible, so every attempt still requires replay. certified_pair means that replay was VALID; inspect its exact claimed gap, which may be zero. Failed attempts prove no uniqueness. No paving upper bound is generated by this bridge.

result.numerical is a snapshot; the original report stays numerical-only and is never mutated. result.to_dict() contains separate numerical and certified sections and can be saved as JSON. Individual certificates in the attempts can be saved and passed to twindex check.

Coordinates correspond in order by default. For a custom parameterisation, provide both theta_map(theta) and target_map, one exact target index per audit coordinate. Each certificate uses its mapped exact target and recomputes the actual rational gap; mapped separation is never assumed. By default mapped values are rounded to 12 decimal places (digits controls this). The exact checker decides prior membership and compatibility after rounding.