Skills › Scientific Agent Skills › 13c-metabolic-flux › _mfa_fit.py
_mfa_fit.py
254 lines · skills/13c-metabolic-flux/scripts/_mfa_fit.py at 5f2bc03
1"""Constrained multistart fitting and nuisance-reoptimized flux profiles."""
2from __future__ import annotations
4import numpy as np
5from scipy.linalg import null_space
6from scipy.optimize import Bounds, LinearConstraint, linprog, minimize
7from scipy.stats import chi2
9from _mfa_model import InputError, require
12def fixed_range(low, high):
13 """Compare a feasible range to its own magnitude, independent of flux units."""
14 return high - low <= 1e-10 * max(abs(low), abs(high), np.finfo(float).tiny)
17class FluxSpace:
18 def __init__(self, network, fixed=None):
19 self.network = network
20 self.a = network.equalities.copy()
21 self.b = network.rhs.copy()
22 if fixed:
23 for rid, value in fixed.items():
24 row = np.eye(len(network.ids))[network.ids.index(rid)]
25 self.a = np.vstack((self.a, row))
26 self.b = np.append(self.b, value)
27 anchors = np.abs(self.b[self.b != 0])
28 self.scale = float(np.median(anchors)) if len(anchors) else float(np.median(network.upper))
29 self.bounds = list(zip(network.lower / self.scale, network.upper / self.scale))
30 feasible = self.linear(np.zeros(len(network.ids)))
31 require(feasible.success, "Flux constraints are infeasible")
32 # A profile endpoint can collapse a face of the feasible polytope.
33 # Promote bounds that fix a reaction throughout this face to equalities.
34 initial_basis = null_space(self.a)
35 for i in range(len(network.ids)):
36 if np.linalg.norm(initial_basis[i]) < 1e-10:
37 continue
38 axis = np.eye(len(network.ids))[i]
39 low, high = self.linear(axis), self.linear(-axis)
40 require(low.success and high.success, "Cannot determine feasible flux ranges")
41 if fixed_range(low.x[i], high.x[i]):
42 self.a = np.vstack((self.a, axis))
43 self.b = np.append(self.b, (high.x[i] + low.x[i]) / 2)
44 feasible = self.linear(np.zeros(len(network.ids)))
45 require(feasible.success, "Flux constraints are infeasible after reducing fixed ranges")
46 self.base = feasible.x
47 self.basis = null_space(self.a) * self.scale
48 self.dimension = self.basis.shape[1]
50 def linear(self, objective):
51 result = linprog(objective, A_eq=self.a, b_eq=self.b / self.scale, bounds=self.bounds,
52 method="highs", options={"time_limit": 10., "primal_feasibility_tolerance": 1e-9,
53 "dual_feasibility_tolerance": 1e-9})
54 if result.success:
55 result.x = result.x * self.scale
56 if not self.feasible(result.x):
57 result.success = False
58 result.message = "LP returned a vector that fails explicit feasibility checks"
59 return result
61 def coordinates(self, flux):
62 return self.basis.T @ (flux - self.base) / self.scale ** 2
64 def flux(self, coordinates):
65 return self.base + self.basis @ coordinates
67 def feasible(self, flux):
68 scale = max(float(np.max(np.abs(flux))), float(np.max(np.abs(self.b))), 1e-300)
69 return self.network.feasible(flux) and np.max(np.abs(self.a @ flux - self.b)) < 1e-8 * scale
71 def starting_points(self, count, rng, warm=None):
72 if not self.dimension:
73 return [self.base]
74 vertices = [self.base]
75 for _ in range(max(12, 3 * self.dimension)):
76 result = self.linear(rng.normal(size=len(self.base)))
77 if result.success:
78 vertices.append(result.x)
79 vertices = np.array(vertices)
80 center = vertices.mean(axis=0)
81 starts = [center]
82 if warm is not None:
83 projected = self.flux(self.coordinates(warm))
84 if self.feasible(projected):
85 starts[0] = projected
86 for _ in range(count - 1):
87 weights = rng.dirichlet(np.ones(len(vertices)) * .3)
88 starts.append(.05 * center + .95 * weights @ vertices)
89 return starts
92def fit(observations, starts=8, seed=2026, maxiter=1000, fixed=None, warm=None):
93 require(type(starts) is int and 1 <= starts <= 100, "starts must be 1..100")
94 require(type(maxiter) is int and 1 <= maxiter <= 10000, "maxiter must be 1..10000")
95 require(observations.size > 0, "Fitting needs measurements")
96 network = observations.network
97 space = FluxSpace(network, fixed)
98 rng = np.random.default_rng(seed)
99 records, solutions = [], []
100 numerical_errors = (InputError, np.linalg.LinAlgError, FloatingPointError, ZeroDivisionError)
102 def rss(flux):
103 residual = observations.residuals(flux)
104 value = float(residual @ residual)
105 require(np.isfinite(value), "Nonfinite objective")
106 return value
108 if not space.dimension:
109 value = rss(space.base)
110 return dict(vector=space.base, rss=value, space=space,
111 starts=[dict(success=True, rss=value, message="All fluxes fixed by constraints")])
112 varying = np.linalg.norm(space.basis / space.scale, axis=1) > 1e-10
113 linear = LinearConstraint(space.basis[varying] / space.scale,
114 (network.lower - space.base)[varying] / space.scale,
115 (network.upper - space.base)[varying] / space.scale)
116 # Loose finite coordinate bounds limit excursions while the linear constraint
117 # enforces all actual reaction bounds. They cannot exclude a feasible vector.
118 radius = np.linalg.norm(network.upper - network.lower) / space.scale + 1
119 for initial in space.starting_points(starts, rng, warm):
120 try:
121 scale = max(1., rss(initial))
123 def objective(z):
124 try:
125 return rss(space.flux(z)) / scale
126 except numerical_errors:
127 return 1e20
129 result = minimize(objective, space.coordinates(initial), method="SLSQP",
130 bounds=Bounds(-radius, radius), constraints=[linear],
131 options={"maxiter": maxiter, "ftol": 1e-12})
132 vector = space.flux(result.x)
133 value = rss(vector)
134 success = bool(result.success and space.feasible(vector) and np.isfinite(value))
135 records.append(dict(success=success, rss=value, message=str(result.message),
136 iterations=int(result.nit)))
137 if success:
138 solutions.append((value, vector))
139 except numerical_errors as exc:
140 records.append(dict(success=False, rss=None, message=str(exc)))
141 require(bool(solutions), "No converged feasible fit; check flux bounds, zero-throughput pools, and starting conditions")
142 value, vector = min(solutions, key=lambda item: item[0])
143 return dict(vector=vector, rss=value, space=space, starts=records)
146def diagnostics(observations, fitted):
147 network, vector, space = observations.network, fitted["vector"], fitted["space"]
148 residual, predictions = observations.residuals(vector, details=True)
149 columns = []
150 z = space.coordinates(vector)
151 for i in range(space.dimension):
152 step = np.eye(space.dimension)[i] * 1e-5
153 plus, minus = space.flux(z + step), space.flux(z - step)
154 if space.feasible(plus) and space.feasible(minus):
155 columns.append((observations.residuals(plus) - observations.residuals(minus)) / 2e-5)
156 elif space.feasible(plus):
157 columns.append((observations.residuals(plus) - residual) / 1e-5)
158 elif space.feasible(minus):
159 columns.append((residual - observations.residuals(minus)) / 1e-5)
160 else:
161 columns = None
162 break
163 rank, singular, null_directions = None, [], []
164 if columns is not None:
165 jacobian = np.column_stack(columns) if columns else np.empty((observations.size, 0))
166 _, values, vh = np.linalg.svd(jacobian, full_matrices=True)
167 threshold = max(1e-4, (float(values[0]) if len(values) else 0.) * 1e-6)
168 rank = int(np.sum(values > threshold))
169 singular = values.tolist()
170 for direction in vh[rank:]:
171 flux_direction = space.basis @ direction
172 flux_direction /= np.max(np.abs(flux_direction))
173 null_directions.append({rid: float(v) for rid, v in zip(network.ids, flux_direction)
174 if abs(v) > 1e-5})
175 active = [rid for i, rid in enumerate(network.ids) if rid not in network.fixed and
176 (abs(vector[i] - network.lower[i]) < 1e-6 * space.scale or
177 abs(vector[i] - network.upper[i]) < 1e-6 * space.scale)]
178 dof = observations.size - space.dimension
179 regular = rank == space.dimension and dof > 0 and not active
180 warnings = list(observations.warnings)
181 if rank is None:
182 warnings.append("Sensitivity rank could not be evaluated along feasible directions; inspect profiles")
183 elif rank < space.dimension:
184 warnings.append("Locally rank-deficient fit: some flux combinations are not constrained by these measurements")
185 if active:
186 warnings.append("Active flux bounds invalidate ordinary interior asymptotic confidence claims")
187 if any(not record["success"] for record in fitted["starts"]):
188 warnings.append("Some optimizer starts failed; inspect start records and repeat with more starts")
189 return dict(
190 fluxes={rid: float(v) for rid, v in zip(network.ids, vector)}, flux_unit=network.spec["flux_unit"],
191 rss=fitted["rss"], independent_measurements=observations.size,
192 free_flux_dimensions=space.dimension, local_sensitivity_rank=rank,
193 sensitivity_singular_values=singular, weak_flux_directions=null_directions,
194 active_bounds=active, max_mass_balance_error=float(np.max(np.abs(network.balance @ vector))),
195 goodness_of_fit=dict(degrees_of_freedom=dof,
196 approximate_p_value=float(chi2.sf(fitted["rss"], dof)) if regular else None,
197 interpretation="Asymptotic Gaussian diagnostic; not evidence that the biological model is correct"
198 if regular else "Not reported: rank, bounds, or degrees of freedom violate regular assumptions"),
199 predictions=predictions, optimizer_starts=fitted["starts"], warnings=warnings)
202def profile(observations, fitted, reaction, points=21, starts=4, seed=2026,
203 maxiter=1000, confidence=.95):
204 network, space = observations.network, fitted["space"]
205 require(reaction in network.ids, f"Unknown profile reaction {reaction}")
206 require(type(points) is int and 5 <= points <= 201, "profile points must be 5..201")
207 require(0 < confidence < 1, "confidence must be between zero and one")
208 index = network.ids.index(reaction)
209 objective = np.eye(len(network.ids))[index]
210 lower = space.linear(objective)
211 upper = space.linear(-objective)
212 require(lower.success and upper.success, "Cannot compute feasible profile range")
213 lo, hi = float(lower.x[index]), float(upper.x[index])
214 if fixed_range(lo, hi):
215 return dict(reaction=reaction, status="fixed_by_constraints", feasible_range=[lo, hi], points=[])
216 grid = np.unique(np.append(np.linspace(lo, hi, points), fitted["vector"][index]))
217 records = []
218 cutoff = float(chi2.ppf(confidence, 1))
219 for i, value in enumerate(grid):
220 try:
221 result = fit(observations, starts=starts, seed=seed + i, maxiter=maxiter,
222 fixed={reaction: float(value)}, warm=fitted["vector"])
223 delta = result["rss"] - fitted["rss"]
224 records.append(dict(flux=float(value), rss=result["rss"], delta_rss=float(delta),
225 accepted=bool(delta <= cutoff), converged=True,
226 failed_starts=sum(not start["success"] for start in result["starts"])))
227 except (InputError, np.linalg.LinAlgError, FloatingPointError) as exc:
228 records.append(dict(flux=float(value), rss=None, delta_rss=None,
229 accepted=None, converged=False, error=str(exc)))
230 valid = [r for r in records if r["converged"]]
231 better = any(r["delta_rss"] < -max(1e-5, abs(fitted["rss"]) * 1e-5) for r in valid)
232 failed = any(not r["converged"] for r in records)
233 all_accepted = len(valid) == len(records) and all(r["accepted"] for r in valid)
234 crossings = []
235 for left, right in zip(records, records[1:]):
236 if left["converged"] and right["converged"] and left["accepted"] != right["accepted"]:
237 crossings.append([left["flux"], right["flux"]])
238 if better:
239 status = "baseline_not_optimal_refit_required"
240 elif failed:
241 status = "incomplete_profile"
242 elif all_accepted:
243 status = "unresolved_within_bounds"
244 elif records[0]["accepted"] or records[-1]["accepted"]:
245 status = "bound_limited"
246 else:
247 status = "threshold_crossings_bracketed"
248 return dict(reaction=reaction, status=status, feasible_range=[lo, hi],
249 confidence=confidence, delta_rss_cutoff=cutoff,
250 threshold_crossing_brackets=crossings,
251 lower_bound_accepted=records[0]["accepted"], upper_bound_accepted=records[-1]["accepted"],
252 note="Grid profile with nuisance fluxes reoptimized. Crossings are brackets, not exact confidence limits. "
253 "Chi-square threshold is asymptotic and requires a suitable error model and regularity; "
254 "a flat finite profile is not proof of global structural non-identifiability.", points=records)