"""Conditional ideal-output algebra for a supplied three-field Oregonator. The finite diffusion matrix, dimensionless clock and concentration observation are premises. This detached calculation runs no trajectory, installs no TNFR law and does not calibrate a chemical experiment or differentiate noisy data. """ from __future__ import annotations from dataclasses import dataclass from fractions import Fraction from .._exact_time import exact_or_represented_real from ..mathematics._exact_linear_algebra import _ordered, _rational_vector __all__ = [ "ChemicalObservationModel", "ChemicalState", "ChemicalOutputJet", "chemical_output_jet", "reconstruct_chemical_state", ] Vector = tuple[Fraction, ...] Matrix = tuple[Vector, ...] _MAX_SITES = 32 def _vector(values, label: str, *, size: int | None = None, nonnegative=False): result = _rational_vector(values, label) if not 1 <= len(result) <= _MAX_SITES: raise ValueError(f"{label} must contain between 1 and {_MAX_SITES} entries") if size is not None and len(result) != size: raise ValueError(f"{label} must have {size} entries") if nonnegative and any(value < 0 for value in result): raise ValueError(f"{label} must be nonnegative") return result def _diffusion(values) -> Matrix: rows = _ordered(values, "diffusion_matrix") if not 1 <= len(rows) <= _MAX_SITES: raise ValueError(f"diffusion_matrix must have between 1 and {_MAX_SITES} rows") matrix = tuple( _vector(row, f"diffusion_matrix[{i}]", size=len(rows)) for i, row in enumerate(rows) ) if any( matrix[i][j] < 0 for i in range(len(matrix)) for j in range(len(matrix)) if i != j ): raise ValueError("diffusion_matrix off-diagonal entries must be nonnegative") if any(sum(row, Fraction(0)) != 0 for row in matrix): raise ValueError("diffusion_matrix row sums must be exactly zero") return matrix def _positive(value, label: str) -> Fraction: result = exact_or_represented_real(value, label) if result <= 0: raise ValueError(f"{label} must be positive") return result def _nonnegative(value, label: str) -> Fraction: result = exact_or_represented_real(value, label) if result < 0: raise ValueError(f"{label} must be nonnegative") return result @dataclass(frozen=True) class ChemicalObservationModel: """Supplied constant coefficients on at most 32 ordered finite sites. ``diffusion_matrix`` has nonnegative off-diagonal entries and zero row sums. Symmetry, continuum convergence and physical no-flux realization are not inferred from this finite matrix admission. """ diffusion_matrix: Matrix epsilon_u: Fraction epsilon_w: Fraction q: Fraction f: Fraction delta: Fraction def __post_init__(self): object.__setattr__(self, "diffusion_matrix", _diffusion(self.diffusion_matrix)) for name in ("epsilon_u", "epsilon_w", "q"): object.__setattr__(self, name, _positive(getattr(self, name), name)) for name in ("f", "delta"): object.__setattr__(self, name, _nonnegative(getattr(self, name), name)) @dataclass(frozen=True) class ChemicalState: """Nonnegative, ordered concentration coordinates in the supplied chart.""" u: Vector v: Vector w: Vector def __post_init__(self): u = _vector(self.u, "u", nonnegative=True) v = _vector(self.v, "v", size=len(u), nonnegative=True) w = _vector(self.w, "w", size=len(u), nonnegative=True) for name, value in (("u", u), ("v", v), ("w", w)): object.__setattr__(self, name, value) @dataclass(frozen=True) class ChemicalOutputJet: """Ideal ``v`` and its first three derivatives in the declared clock. Direct construction admits the scalar/vector domains only; it does not certify that an arbitrary record is generated by a chemical law. """ y: Vector y1: Vector y2: Vector y3: Vector def __post_init__(self): y = _vector(self.y, "y", nonnegative=True) object.__setattr__(self, "y", y) for name in ("y1", "y2", "y3"): object.__setattr__( self, name, _vector(getattr(self, name), name, size=len(y)) ) def _admit_model(model) -> ChemicalObservationModel: if not isinstance(model, ChemicalObservationModel): raise TypeError("model must be a ChemicalObservationModel") # Frozen records can still be forged or altered through low-level access. # Rebuild consumed primitives rather than trusting an earlier constructor. return ChemicalObservationModel( model.diffusion_matrix, model.epsilon_u, model.epsilon_w, model.q, model.f, model.delta, ) def _product(matrix: Matrix, vector: Vector) -> Vector: return tuple( sum( (coefficient * value for coefficient, value in zip(row, vector)), Fraction(0), ) for row in matrix ) def chemical_output_jet(model, u, v, w, illumination) -> ChemicalOutputJet: """Evaluate ``(v, v', v'', v''')`` from re-admitted primitive state. The supplied field is ``u'=(u-u**2+w*(q-u))/epsilon_u + L*u``, ``v'=u-v``, and ``w'=(illumination+f*v-w*(q+u))/epsilon_w + delta*L*w``. Illumination is an ordered nonnegative vector. Its derivative is not consumed at this jet order; no constant-input trajectory is asserted. """ admitted = _admit_model(model) size = len(admitted.diffusion_matrix) state = ChemicalState(u, v, w) if len(state.u) != size: raise ValueError("concentration vectors must match diffusion_matrix dimension") phi = _vector(illumination, "illumination", size=size, nonnegative=True) u, v, w = state.u, state.v, state.w lu = _product(admitted.diffusion_matrix, u) lw = _product(admitted.diffusion_matrix, w) du = tuple( (a - a * a + c * (admitted.q - a)) / admitted.epsilon_u + diffusion for a, c, diffusion in zip(u, w, lu) ) dw = tuple( (light + admitted.f * b - c * (admitted.q + a)) / admitted.epsilon_w + admitted.delta * diffusion for a, b, c, light, diffusion in zip(u, v, w, phi, lw) ) y1 = tuple(a - b for a, b in zip(u, v)) y2 = tuple(rate - first for rate, first in zip(du, y1)) ldu = _product(admitted.diffusion_matrix, du) ddu = tuple( ((1 - 2 * a - c) * rate + (admitted.q - a) * hidden_rate) / admitted.epsilon_u + diffusion for a, c, rate, hidden_rate, diffusion in zip(u, w, du, dw, ldu) ) return ChemicalOutputJet(v, y1, y2, tuple(a - b for a, b in zip(ddu, y2))) def reconstruct_chemical_state(model, y, y1, y2, *, y3=None) -> ChemicalState: """Recover state algebraically from an ideal full-site output jet. At a regular site ``u != q``, derivatives through order two determine ``w``. At ``u == q`` the order-two numerator must vanish; supplied third derivatives enter this pointwise inverse when ``u' != 0``. An unsupported singular site rejects the entire reconstruction; this does not exclude reconstruction using neighboring observations or additional information. All concentrations must be nonnegative. No numerical differentiation or noise tolerance is inserted. Optional ``y3`` is admitted everywhere but consumed only at singular sites. The inverse does not certify a supplied third jet at regular sites, or recover/admit an illumination history. State reconstruction is conditional on the declared ideal observation and complete law. """ admitted = _admit_model(model) size = len(admitted.diffusion_matrix) y = _vector(y, "y", size=size, nonnegative=True) y1 = _vector(y1, "y1", size=size) y2 = _vector(y2, "y2", size=size) third = None if y3 is None else _vector(y3, "y3", size=size) u = tuple(value + first for value, first in zip(y, y1)) if any(value < 0 for value in u): raise ValueError("reconstructed u must be nonnegative") du = tuple(first + second for first, second in zip(y1, y2)) lu = _product(admitted.diffusion_matrix, u) ldu = _product(admitted.diffusion_matrix, du) if third is not None else None hidden = [] for i, (value, rate, diffusion) in enumerate(zip(u, du, lu)): numerator = admitted.epsilon_u * (rate - diffusion) - value + value * value denominator = admitted.q - value if denominator: hidden.append(numerator / denominator) continue if numerator: raise ValueError( f"inconsistent singular chemical output at site {i}: numerator != 0" ) if third is None or rate == 0: raise ValueError( f"unsupported singular reconstruction at site {i}: " "y3 and nonzero u' required by this inverse" ) derivative = ( admitted.epsilon_u * (y2[i] + third[i] - ldu[i]) + (2 * value - 1) * rate ) hidden.append(-derivative / rate) return ChemicalState(u, y, tuple(hidden))