diff --git a/pybnf/algorithms/optimizers/gradient_base.py b/pybnf/algorithms/optimizers/gradient_base.py index 6912857c..4f999562 100644 --- a/pybnf/algorithms/optimizers/gradient_base.py +++ b/pybnf/algorithms/optimizers/gradient_base.py @@ -47,7 +47,7 @@ from .concurrent_multistart import DONE, ConcurrentMultiStartOptimizer from ...gradient import ( GradientNotSupported, - apply_routing, + apply_routings, assemble_constraint_gradient, assemble_gaussian_gradient, route_for_model, @@ -379,31 +379,41 @@ def _setup_gradient_path(self): routings -- idempotent, called once from the leaf's ``start_run`` (before the model scatter, so the request rides the pickle to the workers). - For each model: ``apply_routing`` the **wildtype** routing, whose request - lists (``sensitivity_params`` / ``sensitivity_ic``) are the union over - conditions -- the wildtype pins nothing, so every parameter-/IC-bound free - parameter is requested, a superset of any single condition's (a ``=``-pinned - condition only ever *drops* columns). Then build one - :class:`ExperimentRouting` per scored ``(model, suffix)``, carrying that - condition's chain-rule factors for the assembly. Raises (capability gate) if - the bngsim build lacks ``output_sensitivities``.""" + For each model: build one :class:`ExperimentRouting` per scored + ``(model, suffix)`` (carrying that condition's chain-rule factors for the + assembly), then ``apply_routings`` the **union** of their sensitivity + requests -- plus the wildtype's. The wildtype alone is not a superset: a + condition can route a free parameter to a column no other experiment binds (a + per-condition estimated initial condition, ADR-0076), so that column is + reached only through the condition and must be unioned in (an extra requested + column is harmless; a missing one aborts the assembly). Raises (capability + gate) if the bngsim build lacks ``output_sensitivities``.""" if self._routings is not None: return names = [v.name for v in self.variables] routings = {} for model in self.model_list: - # Union sensitivity request (wildtype) -> sets _sensitivity_request, - # which survives the scatter and is applied at every simulate(). - apply_routing(model, route_for_model(model, names, condition=None)) # Declare which of this model's outputs are scored gradient targets so # an incidental/unscored action (a stochastic diagnostic, a # carried-state pre-equilibration scan) runs sensitivity-free instead # of aborting the whole fit at a differentiability guard (#475). Rides # the scatter alongside the sensitivity request. model.set_scored_suffixes(self.exp_data.get(model.name, {})) + # Apply the UNION sensitivity request over the wildtype and every scored + # condition -> sets _sensitivity_request, which survives the scatter and is + # applied at every simulate(). The wildtype alone is NOT a superset once a + # condition routes a free parameter to a column no other experiment binds -- + # a per-condition estimated initial condition (ADR-0076): its species-IC / + # multiplier column is reached only through that condition, so the union must + # include every scored routing (an extra column is harmless; a missing one aborts). + wildtype = route_for_model(model, names, condition=None) + model_routings = {} for suffix in self.exp_data.get(model.name, {}): condition = self._condition_for_suffix(model, suffix) - routings[(model.name, suffix)] = route_for_model(model, names, condition) + model_routings[suffix] = route_for_model(model, names, condition) + apply_routings(model, [wildtype, *model_routings.values()]) + for suffix, routing in model_routings.items(): + routings[(model.name, suffix)] = routing self._routings = routings def _condition_for_suffix(self, model, suffix): diff --git a/pybnf/bngsim_model/expressions.py b/pybnf/bngsim_model/expressions.py index f466dbf7..9998437b 100644 --- a/pybnf/bngsim_model/expressions.py +++ b/pybnf/bngsim_model/expressions.py @@ -172,6 +172,38 @@ def _parse_net_species_initializers(net_lines): return initializers +def _net_species_ic_seed_map(species_initializers, param_names): + """Map a model parameter to the species whose initial value it *bares* in a .net species + block, for the gradient router's per-condition estimated initial condition composition + (ADR-0076, #511). + + A **bare** initializer ``species <- `` (the whole expression is a single parameter + id) has ``d(IC)/d(param) = 1``, so a free parameter a condition assigns to that parameter + routes straight onto the species' initial-condition sensitivity axis. Such a parameter maps + to its species. A parameter referenced by a **non-bare** initializer expression (``2*p``, + ``a+b``), or that bares more than one species, maps to ``None`` -- present but non-routable + (its ``d(IC)/d(param)`` is not a plain ``1``), so the router refuses rather than emitting a + wrong column. + """ + param_set = set(param_names) + bare = {} # param -> set(species) it bares + nonbare = set() # params a non-bare initializer references (force a refuse) + ident = re.compile(r'^[A-Za-z_]\w*$') + token = re.compile(r'[A-Za-z_]\w*') + for species, expr in species_initializers: + e = expr.strip() + if ident.match(e) and e in param_set: + bare.setdefault(e, set()).add(species) + else: + nonbare.update(t for t in token.findall(e) if t in param_set) + seed_map = {} + for param, species in bare.items(): + seed_map[param] = next(iter(species)) if len(species) == 1 else None + for param in nonbare: + seed_map[param] = None # a non-bare use forces a refuse, even if also bare elsewhere + return seed_map + + def _parse_bngl_param_block(model_lines): """Extract BNGL parameter definitions as ordered (name, expression) pairs.""" params = [] diff --git a/pybnf/bngsim_model/net_model.py b/pybnf/bngsim_model/net_model.py index 005bd16c..6f8d7db6 100644 --- a/pybnf/bngsim_model/net_model.py +++ b/pybnf/bngsim_model/net_model.py @@ -40,6 +40,7 @@ _eval_model_expression, _build_mutant_param_set, _parse_net_species_initializers, + _net_species_ic_seed_map, ) from .scan import ( _with_sim_timeout, @@ -343,7 +344,7 @@ def has_discrete_events(self): def sensitivity_entity_namespace(self): """The bind-by-id namespaces the gradient router classifies free parameters against (#448). - Returns ``(param_ids, species_initializers)``: + Returns ``(param_ids, species_initializers, ic_seed_map)``: * ``param_ids`` -- the model's ``begin parameters`` namespace (the engine's ``param_names``), the kinetic/global ids a free parameter binds to via ``set_param`` @@ -351,12 +352,19 @@ def sensitivity_entity_namespace(self): * ``species_initializers`` -- the ``(species, initial-expr)`` pairs (``_parse_net_species_initializers``); a free parameter that is a species' bare initial-value expression binds to ``Simulator(sensitivity_ic=)`` keyed by the - species (an IC parameter is absent from the ODE RHS, so the parameter axis is zero). + species (an IC parameter is absent from the ODE RHS, so the parameter axis is zero); + * ``ic_seed_map`` -- ``{model parameter -> species}`` for a bare species initializer + (``_net_species_ic_seed_map``), so the router can route a free parameter a condition + assigns to that parameter onto the species' IC axis (a per-condition estimated initial + condition, ADR-0076, #511); a non-routable seed maps to ``None``. This is the only model coupling :mod:`pybnf.gradient.routing` needs, so the routing - core stays backend-agnostic. No simulation -- both namespaces are known at build time. + core stays backend-agnostic. No simulation -- all three are known at build time. """ - return list(self._engine_model.param_names), list(self._net_species_initializers) + param_ids = list(self._engine_model.param_names) + species_initializers = list(self._net_species_initializers) + return (param_ids, species_initializers, + _net_species_ic_seed_map(species_initializers, param_ids)) def _sensitivity_request_kwargs(self, method): """Simulator kwargs requesting forward sensitivities on the gradient path. diff --git a/pybnf/bngsim_sbml_model.py b/pybnf/bngsim_sbml_model.py index f8db418e..d4c42717 100644 --- a/pybnf/bngsim_sbml_model.py +++ b/pybnf/bngsim_sbml_model.py @@ -212,6 +212,7 @@ def _extract_sbml_structure(self): self._initial_dep_names = self._compute_initial_dependency_names(doc.getModel()) self._species_unit_factor, self._unsafe_volume = self._compute_species_unit_factors( doc.getModel()) + self._ic_seed_map = self._compute_ic_seed_map(doc.getModel()) @staticmethod def _collect_ast_names(node, out): @@ -289,6 +290,50 @@ def _compute_initial_dependency_names(self, sbml_model): worklist.append(ref) return dep + def _compute_ic_seed_map(self, sbml_model): + """Map a model parameter to the species whose initial value it *bares*. + + A **bare** ``initialAssignment`` ``species = `` (a single ``ci`` naming a global + parameter, on a concentration-based / unit-factor-1 species) has ``d(IC)/d(param) = 1``, + so the gradient router can route a free parameter a condition assigns to that parameter + (a per-condition estimated initial condition, ADR-0076, #511) straight onto the species' + initial-condition sensitivity axis. Such a parameter maps to its species. + + A parameter that seeds a species initial value through a **non-bare** expression + (``2*init``, ``a+b``), or on an amount species whose value needs a non-unit + concentration factor, or that bares more than one species, maps to ``None`` -- present + but non-routable: its ``d(IC)/d(param)`` is not a plain ``1``, so the router refuses + (``classify_condition_target``) rather than emitting a wrong/silently-zero column. + """ + param_set = set(self._global_param_names) + bare = {} # param -> set(species) it bares + nonbare = set() # params a species IA references non-baruely (force a refuse) + for i in range(sbml_model.getNumInitialAssignments()): + ia = sbml_model.getInitialAssignment(i) + symbol = ia.getSymbol() + if symbol not in self._species_name_set: + continue + math = ia.getMath() + is_bare = ( + math is not None + and math.getType() == libsbml.AST_NAME + and math.getNumChildren() == 0 + and math.getName() in param_set + and self._species_unit_factor.get(symbol, 1.0) == 1.0 + ) + if is_bare: + bare.setdefault(math.getName(), set()).add(symbol) + else: + refs = set() + self._collect_ast_names(math, refs) + nonbare.update(r for r in refs if r in param_set) + seed_map = {} + for param, species in bare.items(): + seed_map[param] = next(iter(species)) if len(species) == 1 else None + for param in nonbare: + seed_map[param] = None # a non-bare use forces a refuse, even if also bare elsewhere + return seed_map + def _compute_species_unit_factors(self, sbml_model): """Per-species factor converting a PyBNF species value to a bngsim concentration, plus a flag for volumes we cannot safely convert. @@ -909,7 +954,7 @@ def has_discrete_events(self): def sensitivity_entity_namespace(self): """The bind-by-id namespaces the gradient router classifies free parameters against (#448). - Returns ``(param_ids, species_initializers)``: + Returns ``(param_ids, species_initializers, ic_seed_map)``: * ``param_ids`` -- the model's global ``parameter`` ids, the kinetic ids a free parameter binds to via ``set_param`` and thus routes to @@ -919,14 +964,19 @@ def sensitivity_entity_namespace(self): for a species sets that species' initial value (via ``set_concentration``, the bind-by-id convention, ADR-0034), so each species' bare initializer expression *is* its own name; such a free parameter routes to the initial-condition axis keyed by the - species (an IC parameter is absent from the ODE RHS, so its parameter axis is zero). + species (an IC parameter is absent from the ODE RHS, so its parameter axis is zero); + * ``ic_seed_map`` -- ``{model parameter -> species}`` for a bare ``initialAssignment`` + ``species = `` (:meth:`_compute_ic_seed_map`), so the router can route a free + parameter a condition assigns to that parameter onto the species' IC axis (a + per-condition estimated initial condition, ADR-0076, #511); a non-routable seed maps + to ``None``. This is the only model coupling :mod:`pybnf.gradient.routing` needs, so the routing core - stays backend-agnostic. No simulation -- both namespaces are known at load time. (A - species whose initial is a non-trivial *expression* of a parameter is the deferred - non-bare-initializer case the router documents, exactly as for the net backend.) + stays backend-agnostic. No simulation -- all three are known at load time. """ - return list(self._global_param_names), [(s, s) for s in self._species_names] + return (list(self._global_param_names), + [(s, s) for s in self._species_names], + dict(self._ic_seed_map)) def _sensitivity_request_kwargs(self, method): """Simulator kwargs requesting forward sensitivities on the gradient path. diff --git a/pybnf/gradient/__init__.py b/pybnf/gradient/__init__.py index 78b57977..3f600679 100644 --- a/pybnf/gradient/__init__.py +++ b/pybnf/gradient/__init__.py @@ -28,13 +28,16 @@ PARAM, IC, NONE, + RouteContribution, ParamRoute, ExperimentRouting, classify_free_param, + classify_condition_target, condition_factor, route_experiment, route_for_model, apply_routing, + apply_routings, ) from .assembly import ( GradientResult, @@ -50,13 +53,16 @@ 'PARAM', 'IC', 'NONE', + 'RouteContribution', 'ParamRoute', 'ExperimentRouting', 'classify_free_param', + 'classify_condition_target', 'condition_factor', 'route_experiment', 'route_for_model', 'apply_routing', + 'apply_routings', 'GradientResult', 'assemble_gaussian_gradient', 'assemble_gradient_and_fisher_hessian', diff --git a/pybnf/gradient/assembly.py b/pybnf/gradient/assembly.py index 86f8c15c..e35b0f39 100644 --- a/pybnf/gradient/assembly.py +++ b/pybnf/gradient/assembly.py @@ -469,24 +469,25 @@ def _accumulate_experiment_fisher(objective, sim_data, exp_data, routing, index, _accumulate_fisher_point(objective, sim_data, exp_data, index, point, hessian) -def _sensitivity(sens, selector, route, sim_row): +def _sensitivity(sens, selector, contribution, sim_row, free_param): """The native forward sensitivity ``d(observable)/d(routed entity)`` at one time row. - Reads the parameter axis (``sensitivity_params``) for a PARAM route and the - initial-condition axis (``sensitivity_ic``) for an IC route, addressing the - entity by the route's ``key`` (parameter id, or species for an IC).""" - if route.target == PARAM: + Reads the parameter axis (``sensitivity_params``) for a PARAM contribution and the + initial-condition axis (``sensitivity_ic``) for an IC contribution, addressing the entity by + the contribution's ``key`` (parameter id, or species for an IC). ``free_param`` names the + routed free parameter for the diagnostic when a requested column is absent.""" + if contribution.target == PARAM: axis, labels = 'parameter', sens.param_names else: # IC axis, labels = 'ic', sens.ic_species - if route.key not in labels: + if contribution.key not in labels: raise GradientNotSupported( "Free parameter '%s' routes to %s '%s', but the simulation's " "sensitivity tensor has no such column (axis labels: %s). Apply the same " "routing to the model before running it (apply_routing)." - % (route.free_param, axis, route.key, ', '.join(map(str, labels)) or '(none)')) + % (free_param, axis, contribution.key, ', '.join(map(str, labels)) or '(none)')) column = sens.slice_for(selector, axis=axis) # (n_times, n_axis) - return column[sim_row, labels.index(route.key)] + return column[sim_row, labels.index(contribution.key)] def _raw_sensitivity_accessor(objective, sim_data, sens, routing, index, n_param, indvar): @@ -522,9 +523,12 @@ def tensor_sens(col_name, row): vec = np.zeros(n_param) selector = _selector_for(sens, col_name) for name, route in routing.routes.items(): - if route.factor == 0.0 or route.target == NONE: - continue - vec[index[name]] = route.factor * _sensitivity(sens, selector, route, row) + total = 0.0 + for c in route.contributions: + if c.target == NONE or c.factor == 0.0: + continue + total += c.factor * _sensitivity(sens, selector, c, row, name) + vec[index[name]] = total return vec def raw_sens(col_name, row): @@ -743,8 +747,11 @@ def raw_sens(model, suffix, observable, row): selector = _selector_for(sens, observable) vec = np.zeros(n_param) for name, route in routing.routes.items(): - if route.factor == 0.0 or route.target == NONE: - continue - vec[index[name]] = route.factor * _sensitivity(sens, selector, route, row) + total = 0.0 + for c in route.contributions: + if c.target == NONE or c.factor == 0.0: + continue + total += c.factor * _sensitivity(sens, selector, c, row, name) + vec[index[name]] = total return vec return raw_sens diff --git a/pybnf/gradient/routing.py b/pybnf/gradient/routing.py index 96a8dec5..6af1b9f1 100644 --- a/pybnf/gradient/routing.py +++ b/pybnf/gradient/routing.py @@ -46,12 +46,24 @@ Multiple mutations on the same id compose: affine maps compose, so the multiplicative parts multiply; any ``=`` pins the parameter, driving the composed factor to ``0``. +Reached through a condition (a per-condition estimated initial condition, ADR-0076) +----------------------------------------------------------------------------------- +A free parameter that binds no model id of its own can still reach the model through a +``condition`` that sets a model entity to *its* value -- ``target = free_param`` +(``is_param_ref``). The referenced free parameter then gains a :class:`RouteContribution` on the +target's own sensitivity column (:func:`classify_condition_target`): :data:`IC` when ``target`` +seeds a species initial value (chain-rule factor ``d(IC)/d(target) = 1`` for a bare +``initialAssignment``), :data:`PARAM` for an ordinary global (e.g. a shared rate multiplier). +Because one free parameter may be assigned to *several* targets in one condition, a route is a +**sum** over its contributions (:class:`ParamRoute`). + Cut-1 scope ----------- -IC routing matches a *bare* initializer (``species <- p``). A parameter that both appears in -the ODE RHS and seeds an initial value, or a species seeded by a non-trivial expression of -the free parameter (``2*p``), is a later layer -- the clean param/ic split is what #449's -first FD-oracle case needs. +IC routing matches a *bare* initializer (``species <- p`` directly, or ``target = p`` where +``target`` bares a species IC). A parameter that both appears in the ODE RHS and seeds an +initial value, or a species seeded by a non-trivial expression (``2*p``), is a later layer -- +:func:`classify_condition_target` raises :class:`GradientNotSupported` for the non-bare seed +rather than emitting a parameter-dependent factor, keeping the clean param/ic split honest. """ from dataclasses import dataclass @@ -68,23 +80,64 @@ @dataclass(frozen=True) -class ParamRoute: - """How one free parameter maps onto an experiment's forward-sensitivity request. +class RouteContribution: + """One ``(native sensitivity column -> free parameter)`` term of a route. ``target`` is :data:`PARAM` (kinetic/global -> ``sensitivity_params``), :data:`IC` - (species initial value -> ``sensitivity_ic``), or :data:`NONE` (bound to no model id -- - a free sigma; no model column). ``key`` is the request key the tensor is read by: the - parameter id for :data:`PARAM`, the *species* for :data:`IC`, ``None`` for :data:`NONE`. - ``factor`` is the chain-rule derivative ``d(perturbed param)/d(free param)`` for this - experiment's condition (#449 multiplies it into the Jacobian column); a pinned (``=``) - parameter has factor ``0`` and is dropped from the request list. + (species initial value -> ``sensitivity_ic``), or :data:`NONE` (no model column). ``key`` + is the request key the tensor is read by: the parameter id for :data:`PARAM`, the *species* + for :data:`IC`, ``None`` for :data:`NONE`. ``factor`` is the chain-rule derivative folded + into this term (#449 multiplies it into the Jacobian column); a zero factor drops it. """ - free_param: str target: str key: object # str (param id / species) for param/ic; None for none factor: float +@dataclass(frozen=True) +class ParamRoute: + """How one free parameter maps onto an experiment's forward-sensitivity request. + + A free parameter's derivative is the **sum** over its ``contributions`` -- one + :class:`RouteContribution` per native sensitivity column it reaches. The common case is a + single contribution: a ``parameter:`` free parameter bound by id (ADR-0034) reaches exactly + one model column. A free parameter routed *only* through a condition (a per-condition + estimated initial condition, ADR-0076) reaches its column through the condition target + instead; and a free parameter a condition assigns to *several* model entities at once (a + shared rate multiplier) reaches several columns, so its derivative is their sum. + + ``.target`` / ``.key`` / ``.factor`` read the sole contribution of a single-column route + (every bind-by-id route); reading them on a multi-column route raises -- use + ``.contributions``. + """ + free_param: str + contributions: tuple # of RouteContribution, in composition order + + @classmethod + def single(cls, free_param, target, key, factor): + """A route with a single :class:`RouteContribution` -- the common bind-by-id case.""" + return cls(free_param, (RouteContribution(target, key, factor),)) + + @property + def target(self): + return self._sole().target + + @property + def key(self): + return self._sole().key + + @property + def factor(self): + return self._sole().factor + + def _sole(self): + if len(self.contributions) != 1: + raise ValueError( + f"ParamRoute for '{self.free_param}' has {len(self.contributions)} " + f"contributions; read .contributions, not .target/.key/.factor.") + return self.contributions[0] + + @dataclass class ExperimentRouting: """The per-experiment routing object: ``{free_param -> ParamRoute}`` plus the derived @@ -94,23 +147,24 @@ class ExperimentRouting: @property def sensitivity_params(self): - """``sensitivity_params=`` for the experiment's Simulator: every parameter-bound free - parameter with a non-zero factor (a pinned ``=`` column is dropped), de-duplicated in - declared free-parameter order.""" + """``sensitivity_params=`` for the experiment's Simulator: every parameter-axis column + any free parameter reaches with a non-zero factor (a pinned ``=`` column is dropped), + de-duplicated in declared free-parameter order.""" return self._request_keys(PARAM) @property def sensitivity_ic(self): - """``sensitivity_ic=`` for the experiment's Simulator: the species of every IC-bound - free parameter with a non-zero factor, de-duplicated in declared free-parameter - order.""" + """``sensitivity_ic=`` for the experiment's Simulator: every species initial-condition + column any free parameter reaches with a non-zero factor, de-duplicated in declared + free-parameter order.""" return self._request_keys(IC) def _request_keys(self, target): keys = [] for route in self.routes.values(): - if route.target == target and route.factor != 0.0 and route.key not in keys: - keys.append(route.key) + for c in route.contributions: + if c.target == target and c.factor != 0.0 and c.key not in keys: + keys.append(c.key) return keys @@ -166,52 +220,123 @@ def classify_free_param(free_param, param_ids, species_initializers): return (NONE, None) -def route_experiment(free_params, param_ids, species_initializers, condition=None): +def classify_condition_target(target, param_ids, species_names, ic_seed_map): + """Classify the model entity a param-ref condition sets: return ``(axis, key, factor)``. + + ``target`` is the model id a ``condition`` assignment ``target = free_param`` sets (a + per-condition estimated initial condition, ADR-0076). A model parameter that seeds a + species' initial value routes to the :data:`IC` axis keyed by that species (checked first -- + an IC seed is *also* a ``begin parameters`` id, but only its IC axis is non-zero); a species + set directly routes to :data:`IC` on itself; any other model parameter routes to + :data:`PARAM`. ``factor`` is ``d(model entity)/d(model parameter)`` -- ``1`` for a bare seed + or a direct set. + + ``ic_seed_map`` maps a model parameter to the species whose initial value it *bares* (a bare + ``initialAssignment`` ``species = ``); the value ``None`` marks a parameter that + seeds a species IC through a **non-bare** expression (a parameter-dependent + ``d(IC)/d(param)`` this bind-by-id routing does not model). Raises + :class:`GradientNotSupported` for that non-bare seed and for a target that binds no + sensitivity entity at all -- keeping a gradient/EFIM fit honest rather than emitting a + silently-zero column. + """ + if target in ic_seed_map: + species = ic_seed_map[target] + if species is None: + raise GradientNotSupported( + f"Condition sets '{target}', which seeds a species initial value whose " + f"d(IC)/d('{target}') is not a plain 1 -- a non-bare initialAssignment " + f"expression, an amount species needing a non-unit concentration factor, or a " + f"parameter seeding several species. The gradient path does not yet route this " + f"per-condition estimated initial condition (ADR-0076). Use a gradient-free " + f"optimizer or sampler for this fit.") + return (IC, species, 1.0) + if target in species_names: + return (IC, target, 1.0) + if target in param_ids: + return (PARAM, target, 1.0) + raise GradientNotSupported( + f"Condition sets '{target}' to the value of a free parameter, but '{target}' is " + f"neither a model parameter nor a species initial value the sensitivity request can " + f"bind; the gradient path cannot route it. Use a gradient-free optimizer or sampler.") + + +def route_experiment(free_params, param_ids, species_initializers, condition=None, + ic_seed_map=None): """Build the :class:`ExperimentRouting` for one experiment (pure -- no model, no sim). ``free_params`` is the ordered free-parameter id list (the config's declared variables); ``param_ids`` the model's ``begin parameters`` namespace; ``species_initializers`` the ``(species, initial-expr)`` pairs; ``condition`` the experiment's - :class:`pybnf.pset.MutationSet` (``None`` for the wildtype experiment). - - A parameter-reference perturbation (a per-condition estimated initial condition, - ADR-0076) raises :class:`GradientNotSupported`: the referenced free parameter reaches the - trajectory *through* the model entity the condition sets (a different id), a coupling this - bind-by-id routing does not model, so its sensitivity column would be silently zero. The - gate keeps a gradient/EFIM fit honest; a gradient-free optimizer/sampler is unaffected. + :class:`pybnf.pset.MutationSet` (``None`` for the wildtype experiment); ``ic_seed_map`` the + ``{model parameter -> species}`` bare-``initialAssignment`` map (:func:`classify_condition_target`). + + A parameter-reference perturbation (a per-condition estimated initial condition, ADR-0076) + ``target = free_param`` **composes** the chain rule: the referenced free parameter reaches + the trajectory through the condition target's own sensitivity column, so it gains a + :class:`RouteContribution` on that column (:data:`IC` when the target seeds a species IC, + :data:`PARAM` for an ordinary global). One free parameter a condition assigns to several + targets at once (a shared rate multiplier) accumulates one contribution per target -- its + derivative is their sum. A target whose ``d(IC)/d(param)`` is not a bare ``1`` (a non-bare + initialAssignment), or that binds no sensitivity entity, raises + :class:`GradientNotSupported` rather than emitting a silently-zero column. """ + param_ids = set(param_ids) + ic_seed_map = ic_seed_map or {} + free_param_set = set(free_params) + species_names = {species for species, _ in species_initializers} + + # Contributions a condition's parameter-reference perturbations add to the *referenced* free + # parameter: it reaches the model through the condition target's own sensitivity column. + ref_contribs = {} # free_param -> [RouteContribution] if condition is not None: for mut in condition: - if getattr(mut, 'is_param_ref', False): + if not getattr(mut, 'is_param_ref', False): + continue + free_param = mut.value + if free_param not in free_param_set: + # References a non-variable; the config layer validates this (Mutation.amount). + continue + if mut.operation != '=': raise GradientNotSupported( - f"Condition sets '{mut.name}' to the value of free parameter '{mut.value}' " - f"(a per-condition estimated initial condition, ADR-0076). The gradient path " - f"cannot yet route a free parameter that reaches the model through a " - f"condition target of a different id; use a gradient-free optimizer or " - f"sampler for this fit.") - param_ids = set(param_ids) + f"Condition perturbs '{mut.name}' {mut.operation} '{mut.value}' by a " + f"non-'=' parameter reference; the gradient path routes only an '=' " + f"per-condition estimated initial condition (ADR-0076). Use a gradient-free " + f"optimizer or sampler for this fit.") + axis, key, factor = classify_condition_target( + mut.name, param_ids, species_names, ic_seed_map) + ref_contribs.setdefault(free_param, []).append( + RouteContribution(axis, key, factor)) + routes = {} for name in free_params: - target, key = classify_free_param(name, param_ids, species_initializers) - routes[name] = ParamRoute( - free_param=name, target=target, key=key, - factor=condition_factor(name, condition), - ) + contribs = [] + base_target, base_key = classify_free_param(name, param_ids, species_initializers) + if base_target != NONE: + contribs.append( + RouteContribution(base_target, base_key, condition_factor(name, condition))) + contribs.extend(ref_contribs.get(name, [])) + if not contribs: + # No model column at all (a free sigma, or a free parameter pinned out of every + # experiment): a single NONE contribution, dropped by the request lists and assembly. + contribs.append(RouteContribution(NONE, None, condition_factor(name, condition))) + routes[name] = ParamRoute(free_param=name, contributions=tuple(contribs)) return ExperimentRouting(routes=routes) def route_for_model(model, free_params, condition=None): """:func:`route_experiment` against a live model's bind-by-id namespaces. - Reads the model's ``begin parameters`` ids and ``(species, initial-expr)`` pairs through + Reads the model's ``begin parameters`` ids, ``(species, initial-expr)`` pairs, and the + bare-``initialAssignment`` seed map through :meth:`BngsimModel.sensitivity_entity_namespace` (the only model coupling), so the routing core stays backend-agnostic. ``condition`` may be a :class:`pybnf.pset.MutationSet`, a condition *name* resolved against ``model.mutants``, or ``None`` for the wildtype experiment. """ - param_ids, species_initializers = model.sensitivity_entity_namespace() + param_ids, species_initializers, ic_seed_map = model.sensitivity_entity_namespace() condition = _resolve_condition(model, condition) - return route_experiment(free_params, param_ids, species_initializers, condition) + return route_experiment(free_params, param_ids, species_initializers, condition, + ic_seed_map=ic_seed_map) def apply_routing(model, routing): @@ -227,6 +352,28 @@ def apply_routing(model, routing): return routing +def apply_routings(model, routings): + """Hand the **union** request over several routings to #447's gradient path on ``model``. + + The model's forward-sensitivity request rides the scatter and is applied at every + simulate(), so it must cover every column any scored experiment reads -- the union of the + per-condition ``sensitivity_params`` / ``sensitivity_ic``. (The wildtype request is *not* a + superset once a condition routes a free parameter to a column no other experiment binds -- + a per-condition estimated initial condition, ADR-0076.) Capability-gated exactly as + :func:`apply_routing`. Returns the applied ``(params, ic)`` lists. + """ + params, ic = [], [] + for routing in routings: + for key in routing.sensitivity_params: + if key not in params: + params.append(key) + for key in routing.sensitivity_ic: + if key not in ic: + ic.append(key) + model.enable_output_sensitivities(params=params, ic=ic) + return params, ic + + def _resolve_condition(model, condition): """Resolve a condition *name* to its :class:`MutationSet` on ``model``; pass anything else (a MutationSet or ``None``) through unchanged.""" diff --git a/tests/test_bngsim_expressions.py b/tests/test_bngsim_expressions.py index 133ac943..5ea226dc 100644 --- a/tests/test_bngsim_expressions.py +++ b/tests/test_bngsim_expressions.py @@ -144,6 +144,31 @@ def test_parse_net_species_initializers(): assert expressions._parse_net_species_initializers(lines) == [('A()', '100'), ('B()', '3.0*Vo')] +def test_net_species_ic_seed_map_bare_param_seed(): + """A bare initializer ``species <- `` maps that parameter to its species (routable, + d(IC)/d(param) = 1), so a condition assigning the parameter to a free parameter can route it + onto the species IC axis (ADR-0076, #511).""" + inits = [('A()', 'initA'), ('B()', 'initB')] + assert expressions._net_species_ic_seed_map(inits, ['initA', 'initB', 'k']) == { + 'initA': 'A()', 'initB': 'B()'} + + +def test_net_species_ic_seed_map_non_bare_and_numeric_are_not_routable(): + """A numeric initializer seeds nothing; a non-bare expression maps its referenced parameters + to None (present-but-non-routable), so the router refuses rather than emitting a + parameter-dependent factor.""" + inits = [('A()', '100'), ('B()', '3.0*Vo'), ('C()', 'k+kf')] + seed = expressions._net_species_ic_seed_map(inits, ['Vo', 'k', 'kf']) + assert seed == {'Vo': None, 'k': None, 'kf': None} + + +def test_net_species_ic_seed_map_multi_species_seed_is_not_routable(): + """A parameter that bares more than one species maps to None: routing it to a single IC axis + would drop the other species' contribution.""" + inits = [('A()', 'init'), ('B()', 'init')] + assert expressions._net_species_ic_seed_map(inits, ['init']) == {'init': None} + + # ---------------------------------------------------------- model expression eval class _FakeModel: def __init__(self, params): diff --git a/tests/test_bngsim_sbml_bridge.py b/tests/test_bngsim_sbml_bridge.py index b5b583a0..c406ae3a 100644 --- a/tests/test_bngsim_sbml_bridge.py +++ b/tests/test_bngsim_sbml_bridge.py @@ -278,15 +278,19 @@ def test_sensitivity_entity_namespace_global_params_and_species_ic(): """The gradient router's bind-by-id namespace (#448/#455): the SBML backend reports its global model parameters as the parameter axis and each species as its own bare initializer, so a free parameter named for a global param routes to ``sensitivity_params`` and one named - for a species routes to ``sensitivity_ic`` keyed by that species.""" + for a species routes to ``sensitivity_ic`` keyed by that species. The third element is the + bare-``initialAssignment`` seed map (ADR-0076, #511) -- empty for this model, which has no + parameter-seeded species initials.""" model = object.__new__(bngsim_sbml_model.BngsimSbmlModelNoTimeout) model._global_param_names = ('kAB', 'kBA') model._species_names = ('A', 'B', 'C') + model._ic_seed_map = {} - param_ids, species_initializers = model.sensitivity_entity_namespace() + param_ids, species_initializers, ic_seed_map = model.sensitivity_entity_namespace() assert param_ids == ['kAB', 'kBA'] assert species_initializers == [('A', 'A'), ('B', 'B'), ('C', 'C')] + assert ic_seed_map == {} def test_make_simulator_threads_sensitivity_request_for_ode(monkeypatch): diff --git a/tests/test_gradient_assembly.py b/tests/test_gradient_assembly.py index 652d77a1..857d133c 100644 --- a/tests/test_gradient_assembly.py +++ b/tests/test_gradient_assembly.py @@ -29,7 +29,7 @@ from pybnf.gradient import ( assemble_constraint_gradient, assemble_constraint_hessian, assemble_fisher_hessian, assemble_gaussian_gradient, GradientNotSupported, - PARAM, IC, NONE, ExperimentRouting, ParamRoute, route_experiment, + PARAM, IC, NONE, ExperimentRouting, ParamRoute, RouteContribution, route_experiment, ) from pybnf.gradient.assembly import _sampling_scale_factors from pybnf.constraint import AlwaysConstraint, AtConstraint, ConstraintSet @@ -160,7 +160,7 @@ def test_residual_jacobian_and_scalar_gradient_param_axis(): obj = ChiSquareObjective() sim = _sim_with_sensitivities(pred, d_param=dk) exp = _exp(obs, sigma) - routing = ExperimentRouting(routes={'k': ParamRoute('k', PARAM, 'k', 1.0)}) + routing = ExperimentRouting(routes={'k': ParamRoute.single('k', PARAM, 'k', 1.0)}) free = _free(('k', 'uniform_var', 0.0, 10.0, 0.3)) res = assemble_gaussian_gradient(obj, [(sim, exp, routing)], free) @@ -188,8 +188,8 @@ def test_initial_condition_axis_and_factor(): exp = _exp(obs, sigma) # k scaled by 4 in this condition (factor 4); S0 unperturbed IC (factor 1). routing = ExperimentRouting(routes={ - 'k': ParamRoute('k', PARAM, 'k', 4.0), - 'S0': ParamRoute('S0', IC, 'S()', 1.0), + 'k': ParamRoute.single('k', PARAM, 'k', 4.0), + 'S0': ParamRoute.single('S0', IC, 'S()', 1.0), }) free = _free(('k', 'uniform_var', 0.0, 10.0, 0.3), ('S0', 'uniform_var', 0.0, 500.0, 120.0)) @@ -203,6 +203,33 @@ def test_initial_condition_axis_and_factor(): np.testing.assert_allclose(res.gradient, J.T @ rho) +def test_multi_contribution_route_sums_columns(): + """A free parameter a condition assigns to several targets at once (a per-condition + estimated initial condition / shared multiplier, ADR-0076, #511) sums the native + sensitivity columns it reaches: its Jacobian column is the sum over its RouteContributions. + Here one free param 'm' is routed to both the parameter axis (k) and the IC axis (S()), + each with chain-rule factor 1, so its column is dk + ds0.""" + pred = np.array([120.0, 80.0, 53.0, 36.0]) + obs = np.array([118.0, 82.0, 50.0, 38.0]) + sigma = 4.0 + dk = np.array([0.0, -80.0, -106.0, -108.0]) # d Stot/d k (parameter axis) + ds0 = np.array([1.0, 0.66, 0.44, 0.30]) # d Stot/d S(0) (ic axis) + + obj = ChiSquareObjective() + sim = _sim_with_sensitivities(pred, d_param=dk, d_ic=ds0) + exp = _exp(obs, sigma) + routing = ExperimentRouting(routes={'m': ParamRoute('m', ( + RouteContribution(PARAM, 'k', 1.0), RouteContribution(IC, 'S()', 1.0)))}) + free = _free(('m', 'uniform_var', 0.0, 10.0, 0.3)) + + res = assemble_gaussian_gradient(obj, [(sim, exp, routing)], free) + rho = (pred - obs) / sigma + J = ((1.0 / sigma) * (dk + ds0)).reshape(-1, 1) # the two columns SUM into m's column + assert res.param_names == ['m'] + np.testing.assert_allclose(res.jacobian, J) + np.testing.assert_allclose(res.gradient, J.T @ rho) + + def test_summation_across_experiments_stacks_residuals_and_sums_gradient(): """Two experiments contribute stacked residual/Jacobian rows and a summed gradient -- equal to assembling each alone and adding the scalar gradients.""" @@ -211,11 +238,11 @@ def test_summation_across_experiments_stacks_residuals_and_sums_gradient(): sim1 = _sim_with_sensitivities([100, 74, 55, 41], d_param=[0, -74, -110, -123]) exp1 = _exp([100, 70, 60, 40], 5.0) - r1 = ExperimentRouting(routes={'k': ParamRoute('k', PARAM, 'k', 1.0)}) + r1 = ExperimentRouting(routes={'k': ParamRoute.single('k', PARAM, 'k', 1.0)}) sim2 = _sim_with_sensitivities([100, 55, 30, 17], d_param=[0, -110, -120, -100]) exp2 = _exp([100, 58, 28, 20], 3.0) - r2 = ExperimentRouting(routes={'k': ParamRoute('k', PARAM, 'k', 4.0)}) + r2 = ExperimentRouting(routes={'k': ParamRoute.single('k', PARAM, 'k', 4.0)}) both = assemble_gaussian_gradient(obj, [(sim1, exp1, r1), (sim2, exp2, r2)], free) g1 = assemble_gaussian_gradient(obj, [(sim1, exp1, r1)], free).gradient @@ -236,9 +263,9 @@ def test_pinned_and_unbound_parameters_have_zero_gradient_columns(): sim = _sim_with_sensitivities([100, 74, 55, 41], d_param=[0, -74, -110, -123]) exp = _exp([100, 70, 60, 40], 5.0) routing = ExperimentRouting(routes={ - 'k': ParamRoute('k', PARAM, 'k', 1.0), - 'kpin': ParamRoute('kpin', PARAM, 'kpin', 0.0), # pinned -> dropped - 'sigma': ParamRoute('sigma', NONE, None, 1.0), # unbound nuisance + 'k': ParamRoute.single('k', PARAM, 'k', 1.0), + 'kpin': ParamRoute.single('kpin', PARAM, 'kpin', 0.0), # pinned -> dropped + 'sigma': ParamRoute.single('sigma', NONE, None, 1.0), # unbound nuisance }) free = _free(('k', 'uniform_var', 0.0, 10.0, 0.3), ('kpin', 'uniform_var', 0.0, 10.0, 0.5), @@ -259,7 +286,7 @@ def test_bootstrap_weight_folds_in_as_sqrt_w(): sim = _sim_with_sensitivities([100, 74, 55, 41], d_param=[0, -74, -110, -123]) exp = _exp([100, 70, 60, 40], 5.0) exp.weights[:, exp.cols['Stot']] = np.array([1.0, 4.0, 0.0, 2.0]) # bootstrap counts - routing = ExperimentRouting(routes={'k': ParamRoute('k', PARAM, 'k', 1.0)}) + routing = ExperimentRouting(routes={'k': ParamRoute.single('k', PARAM, 'k', 1.0)}) free = _free(('k', 'uniform_var', 0.0, 10.0, 0.3)) res = assemble_gaussian_gradient(obj, [(sim, exp, routing)], free) @@ -288,8 +315,8 @@ def test_estimated_sigma_adds_a_scalar_noise_gradient_column(): exp = _exp_dyn(obs) # k -> the parameter axis; noise_sd -> NONE (a free sigma, bound to no model id). routing = ExperimentRouting(routes={ - 'k': ParamRoute('k', PARAM, 'k', 1.0), - 'noise_sd': ParamRoute('noise_sd', NONE, None, 1.0), + 'k': ParamRoute.single('k', PARAM, 'k', 1.0), + 'noise_sd': ParamRoute.single('noise_sd', NONE, None, 1.0), }) free = _free(('k', 'uniform_var', 0.0, 10.0, 0.3), ('noise_sd', 'uniform_var', 0.01, 100.0, sigma)) @@ -321,7 +348,7 @@ def test_fixed_sigma_objective_is_least_squares_exact(): obj = ChiSquareObjective() sim = _sim_with_sensitivities([100, 74, 55, 41], d_param=[0, -74, -110, -123]) exp = _exp([100, 70, 60, 40], 5.0) - routing = ExperimentRouting(routes={'k': ParamRoute('k', PARAM, 'k', 1.0)}) + routing = ExperimentRouting(routes={'k': ParamRoute.single('k', PARAM, 'k', 1.0)}) free = _free(('k', 'uniform_var', 0.0, 10.0, 0.3)) res = assemble_gaussian_gradient(obj, [(sim, exp, routing)], free) @@ -341,10 +368,10 @@ def test_estimated_sigma_gradient_sums_across_experiments(): exp1 = _exp_dyn([100, 70, 60, 40]) sim2 = _sim_with_sensitivities([100, 55, 30, 17], d_param=[0, -110, -120, -100]) exp2 = _exp_dyn([100, 58, 28, 20]) - r1 = ExperimentRouting(routes={'k': ParamRoute('k', PARAM, 'k', 1.0), - 'noise_sd': ParamRoute('noise_sd', NONE, None, 1.0)}) - r2 = ExperimentRouting(routes={'k': ParamRoute('k', PARAM, 'k', 4.0), - 'noise_sd': ParamRoute('noise_sd', NONE, None, 1.0)}) + r1 = ExperimentRouting(routes={'k': ParamRoute.single('k', PARAM, 'k', 1.0), + 'noise_sd': ParamRoute.single('noise_sd', NONE, None, 1.0)}) + r2 = ExperimentRouting(routes={'k': ParamRoute.single('k', PARAM, 'k', 4.0), + 'noise_sd': ParamRoute.single('noise_sd', NONE, None, 1.0)}) both = assemble_gaussian_gradient(obj, [(sim1, exp1, r1), (sim2, exp2, r2)], free) expected = 0.0 @@ -375,7 +402,7 @@ def test_logscale_residual_is_in_additive_space(scale, forward, dforward): obj = _logscale_objective(scale) sim = _sim_with_sensitivities(pred, d_param=dk) exp = _exp(obs, sigma) - routing = ExperimentRouting(routes={'k': ParamRoute('k', PARAM, 'k', 1.0)}) + routing = ExperimentRouting(routes={'k': ParamRoute.single('k', PARAM, 'k', 1.0)}) free = _free(('k', 'uniform_var', 0.0, 10.0, 0.3)) res = assemble_gaussian_gradient(obj, [(sim, exp, routing)], free) @@ -401,7 +428,7 @@ def test_log_scale_collapses_to_linear_when_linear(): obj = ChiSquareObjective() sim = _sim_with_sensitivities(pred, d_param=dk) exp = _exp(obs, sigma) - routing = ExperimentRouting(routes={'k': ParamRoute('k', PARAM, 'k', 1.0)}) + routing = ExperimentRouting(routes={'k': ParamRoute.single('k', PARAM, 'k', 1.0)}) free = _free(('k', 'uniform_var', 0.0, 10.0, 0.3)) res = assemble_gaussian_gradient(obj, [(sim, exp, routing)], free) @@ -425,8 +452,8 @@ def test_estimated_sigma_on_log_scale_composes_d_with_e(): sim = _sim_with_sensitivities(pred, d_param=dk) exp = _exp_dyn(obs) # no _SD column: sigma is the free parameter, not the data routing = ExperimentRouting(routes={ - 'k': ParamRoute('k', PARAM, 'k', 1.0), - 'noise_sd': ParamRoute('noise_sd', NONE, None, 1.0), + 'k': ParamRoute.single('k', PARAM, 'k', 1.0), + 'noise_sd': ParamRoute.single('noise_sd', NONE, None, 1.0), }) free = _free(('k', 'uniform_var', 0.0, 10.0, 0.3), ('noise_sd', 'uniform_var', 0.01, 100.0, sigma)) @@ -461,7 +488,7 @@ def test_logscale_nonpositive_point_propagates_nonfinite(pred, obs): obj = _logscale_objective(LOG10) sim = _sim_with_sensitivities(np.array(pred, float), d_param=[0.0, -74.0, -110.0, -123.0]) exp = _exp(obs, sigma) - routing = ExperimentRouting(routes={'k': ParamRoute('k', PARAM, 'k', 1.0)}) + routing = ExperimentRouting(routes={'k': ParamRoute.single('k', PARAM, 'k', 1.0)}) free = _free(('k', 'uniform_var', 0.0, 10.0, 0.3)) with np.errstate(all='ignore'): @@ -512,7 +539,7 @@ def test_cumulative_prediction_sensitivity_differences_rows(): obj._cumulative_cols = frozenset({'Stot'}) sim = _sim_with_sensitivities(pred, d_param=dk) exp = _exp(obs, sigma) - routing = ExperimentRouting(routes={'k': ParamRoute('k', PARAM, 'k', 1.0)}) + routing = ExperimentRouting(routes={'k': ParamRoute.single('k', PARAM, 'k', 1.0)}) free = _free(('k', 'uniform_var', 0.0, 10.0, 0.3)) res = assemble_gaussian_gradient(obj, [(sim, exp, routing)], free) @@ -544,8 +571,8 @@ def test_per_measurement_scale_offset_chain_rule(): obj = _per_measurement_objective() sim = _sim_with_sensitivities(pred, d_param=dk) exp = _exp_pm(obs, sigma, scale_token='scale', offset_token=3.0) - routing = ExperimentRouting(routes={'k': ParamRoute('k', PARAM, 'k', 1.0), - 'scale': ParamRoute('scale', NONE, None, 1.0)}) + routing = ExperimentRouting(routes={'k': ParamRoute.single('k', PARAM, 'k', 1.0), + 'scale': ParamRoute.single('scale', NONE, None, 1.0)}) free = _free(('k', 'uniform_var', 0.0, 10.0, 0.3), ('scale', 'uniform_var', 0.0, 10.0, scale)) @@ -570,7 +597,7 @@ def test_per_measurement_numeric_tokens_only_touch_model_axis(): obj = _per_measurement_objective() sim = _sim_with_sensitivities(pred, d_param=dk) exp = _exp_pm([200., 150., 120., 90.], 5.0, scale_token=2.0, offset_token=3.0) - routing = ExperimentRouting(routes={'k': ParamRoute('k', PARAM, 'k', 1.0)}) + routing = ExperimentRouting(routes={'k': ParamRoute.single('k', PARAM, 'k', 1.0)}) free = _free(('k', 'uniform_var', 0.0, 10.0, 0.3)) res = assemble_gaussian_gradient(obj, [(sim, exp, routing)], free) @@ -623,8 +650,8 @@ def test_measurement_model_chain_rule(): _materialize(obj, sim, {'w': w}) # adds obs = w*Stot + 3 into sim obs_pred = w * pred + 3.0 exp = _exp_obs(obs_pred - 10.0, sigma) # a constant 10-unit miss -> non-zero rho - routing = ExperimentRouting(routes={'k': ParamRoute('k', PARAM, 'k', 1.0), - 'w': ParamRoute('w', NONE, None, 1.0)}) + routing = ExperimentRouting(routes={'k': ParamRoute.single('k', PARAM, 'k', 1.0), + 'w': ParamRoute.single('w', NONE, None, 1.0)}) free = _free(('k', 'uniform_var', 0.0, 10.0, 0.3), ('w', 'uniform_var', 0.0, 10.0, w)) @@ -649,7 +676,7 @@ def test_measurement_model_numeric_only_touches_model_axis(): sim = _sim_with_sensitivities(pred, d_param=dk) _materialize(obj, sim, {}) exp = _exp_obs(2.0 * pred + 3.0 - 8.0, 5.0) - routing = ExperimentRouting(routes={'k': ParamRoute('k', PARAM, 'k', 1.0)}) + routing = ExperimentRouting(routes={'k': ParamRoute.single('k', PARAM, 'k', 1.0)}) free = _free(('k', 'uniform_var', 0.0, 10.0, 0.3)) res = assemble_gaussian_gradient(obj, [(sim, exp, routing)], free) @@ -666,8 +693,8 @@ def test_capability_gate_now_accepts_measurement_layer(): sim = _sim_with_sensitivities([100., 74., 55., 41.], d_param=[0., -74., -110., -123.]) _materialize(obj, sim, {'w': 1.5}) exp = _exp_obs([150., 110., 80., 60.], 5.0) - routing = ExperimentRouting(routes={'k': ParamRoute('k', PARAM, 'k', 1.0), - 'w': ParamRoute('w', NONE, None, 1.0)}) + routing = ExperimentRouting(routes={'k': ParamRoute.single('k', PARAM, 'k', 1.0), + 'w': ParamRoute.single('w', NONE, None, 1.0)}) free = _free(('k', 'uniform_var', 0.0, 10.0, 0.3), ('w', 'uniform_var', 0.0, 10.0, 1.5)) res = assemble_gaussian_gradient(obj, [(sim, exp, routing)], free) @@ -702,9 +729,9 @@ def _prediction_sigma_setup(k=0.3, sd_abs=5.0, sd_rel=0.1, sim = _sim_with_sensitivities(pred, d_param=dk) exp = _exp_dyn(obs) routing = ExperimentRouting(routes={ - 'k': ParamRoute('k', PARAM, 'k', 1.0), - 'sd_abs': ParamRoute('sd_abs', NONE, None, 1.0), - 'sd_rel': ParamRoute('sd_rel', NONE, None, 1.0)}) + 'k': ParamRoute.single('k', PARAM, 'k', 1.0), + 'sd_abs': ParamRoute.single('sd_abs', NONE, None, 1.0), + 'sd_rel': ParamRoute.single('sd_rel', NONE, None, 1.0)}) free = _free(('k', 'uniform_var', 0.0, 10.0, k), ('sd_abs', 'uniform_var', 0.01, 100.0, sd_abs), ('sd_rel', 'uniform_var', 0.001, 2.0, sd_rel)) @@ -802,9 +829,9 @@ def test_prediction_sigma_referencing_other_species_chains_that_sensitivity(): d_param=np.stack([dstot, daux], axis=1).reshape(4, 2, 1), d_ic=None) exp = _exp_dyn(obs) routing = ExperimentRouting(routes={ - 'k': ParamRoute('k', PARAM, 'k', 1.0), - 'sd_abs': ParamRoute('sd_abs', NONE, None, 1.0), - 'sd_rel': ParamRoute('sd_rel', NONE, None, 1.0)}) + 'k': ParamRoute.single('k', PARAM, 'k', 1.0), + 'sd_abs': ParamRoute.single('sd_abs', NONE, None, 1.0), + 'sd_rel': ParamRoute.single('sd_rel', NONE, None, 1.0)}) free = _free(('k', 'uniform_var', 0.0, 10.0, 0.3), ('sd_abs', 'uniform_var', 0.01, 100.0, sd_abs), ('sd_rel', 'uniform_var', 0.001, 2.0, sd_rel)) @@ -871,7 +898,7 @@ def test_normalization_peak_closed_form(): sim = _sim_with_sensitivities(raw.copy(), d_param=dk) sim.normalize('peak') # rescales Stot in place, records (N, ref_row) exp = _exp(np.zeros(4), sigma) - routing = ExperimentRouting(routes={'k': ParamRoute('k', PARAM, 'k', 1.0)}) + routing = ExperimentRouting(routes={'k': ParamRoute.single('k', PARAM, 'k', 1.0)}) free = _free(('k', 'uniform_var', 0.0, 10.0, 0.3)) res = assemble_gaussian_gradient(obj := ChiSquareObjective(), [(sim, exp, routing)], free) @@ -897,7 +924,7 @@ def test_normalization_chain_rule_matches_finite_difference(method): sim = _sim_with_sensitivities(raw.copy(), d_param=dk) sim.normalize(method) exp = _exp(np.zeros(4), sigma) - routing = ExperimentRouting(routes={'k': ParamRoute('k', PARAM, 'k', 1.0)}) + routing = ExperimentRouting(routes={'k': ParamRoute.single('k', PARAM, 'k', 1.0)}) free = _free(('k', 'uniform_var', 0.0, 10.0, 0.3)) res = assemble_gaussian_gradient(ChiSquareObjective(), [(sim, exp, routing)], free) @@ -921,7 +948,7 @@ def test_floor_normalization_gradient_is_deferred(): sim.normalize(('floor', 0.03)) # records method='floor' assert sim.normalization['Stot'].method == 'floor' exp = _exp(np.zeros(4), 1.0) - routing = ExperimentRouting(routes={'k': ParamRoute('k', PARAM, 'k', 1.0)}) + routing = ExperimentRouting(routes={'k': ParamRoute.single('k', PARAM, 'k', 1.0)}) free = _free(('k', 'uniform_var', 0.0, 10.0, 0.3)) with pytest.raises(GradientNotSupported, match='[Ff]loor'): assemble_gaussian_gradient(ChiSquareObjective(), [(sim, exp, routing)], free) @@ -934,7 +961,7 @@ def test_analytic_scale_gradient_is_deferred(): raw = np.array([2.0, 9.0, 5.0, 3.0]) sim = _sim_with_sensitivities(raw.copy(), d_param=np.array([0.5, -2.0, 1.3, -0.7])) exp = _exp(np.zeros(4), 1.0) - routing = ExperimentRouting(routes={'k': ParamRoute('k', PARAM, 'k', 1.0)}) + routing = ExperimentRouting(routes={'k': ParamRoute.single('k', PARAM, 'k', 1.0)}) free = _free(('k', 'uniform_var', 0.0, 10.0, 0.3)) obj = ChiSquareObjective() obj._analytic_scale = {'expt': frozenset({'Stot'})} @@ -948,7 +975,7 @@ def test_capability_gate_now_accepts_trajectory_transforms(): ``_prediction`` transforms once refused by the gate -- now assemble a residual/Jacobian like any other MEDIAN Gaussian (the trajectory-transform gate clause is gone).""" sim = _sim_with_sensitivities([100, 74, 55, 41], d_param=[0, -74, -110, -123]) - routing = ExperimentRouting(routes={'k': ParamRoute('k', PARAM, 'k', 1.0)}) + routing = ExperimentRouting(routes={'k': ParamRoute.single('k', PARAM, 'k', 1.0)}) free = _free(('k', 'uniform_var', 0.0, 10.0, 0.3)) cum = ChiSquareObjective() @@ -972,7 +999,7 @@ def test_capability_gate_refuses_unsupported_objectives(factory): obj = factory() sim = _sim_with_sensitivities([100, 74, 55, 41], d_param=[0, -74, -110, -123]) exp = _exp([100, 70, 60, 40], 5.0) - routing = ExperimentRouting(routes={'k': ParamRoute('k', PARAM, 'k', 1.0)}) + routing = ExperimentRouting(routes={'k': ParamRoute.single('k', PARAM, 'k', 1.0)}) free = _free(('k', 'uniform_var', 0.0, 10.0, 0.3)) with pytest.raises(GradientNotSupported): assemble_gaussian_gradient(obj, [(sim, exp, routing)], free) @@ -985,7 +1012,7 @@ def test_capability_gate_now_accepts_log_scale_gaussian(): that once refused it is gone).""" sim = _sim_with_sensitivities([100, 74, 55, 41], d_param=[0, -74, -110, -123]) exp = _exp([100, 70, 60, 40], 0.1) - routing = ExperimentRouting(routes={'k': ParamRoute('k', PARAM, 'k', 1.0)}) + routing = ExperimentRouting(routes={'k': ParamRoute.single('k', PARAM, 'k', 1.0)}) free = _free(('k', 'uniform_var', 0.0, 10.0, 0.3)) for obj in (LogNormalObjective(), _logscale_objective(LOG10)): res = assemble_gaussian_gradient(obj, [(sim, exp, routing)], free) @@ -1021,10 +1048,10 @@ def _formula_sigma_setup(k=0.3, sd_abs=5.0, sd_rel=0.1, scaling=2.0, sim = _sim_with_sensitivities(pred, d_param=dk) exp = _exp_dyn(obs) routing = ExperimentRouting(routes={ - 'k': ParamRoute('k', PARAM, 'k', 1.0), - 'sd_abs': ParamRoute('sd_abs', NONE, None, 1.0), - 'sd_rel': ParamRoute('sd_rel', NONE, None, 1.0), - 'scaling': ParamRoute('scaling', NONE, None, 1.0)}) + 'k': ParamRoute.single('k', PARAM, 'k', 1.0), + 'sd_abs': ParamRoute.single('sd_abs', NONE, None, 1.0), + 'sd_rel': ParamRoute.single('sd_rel', NONE, None, 1.0), + 'scaling': ParamRoute.single('scaling', NONE, None, 1.0)}) free = _free(('k', 'uniform_var', 0.0, 10.0, k), ('sd_abs', 'uniform_var', 0.01, 100.0, sd_abs), ('sd_rel', 'uniform_var', 0.001, 2.0, sd_rel), @@ -1109,8 +1136,8 @@ def test_formula_sigma_single_symbol_matches_free_parameter_sigma(): sim = _sim_with_sensitivities(pred, d_param=dk) exp = _exp_dyn(obs) routing = ExperimentRouting(routes={ - 'k': ParamRoute('k', PARAM, 'k', 1.0), - 'sig': ParamRoute('sig', NONE, None, 1.0)}) + 'k': ParamRoute.single('k', PARAM, 'k', 1.0), + 'sig': ParamRoute.single('sig', NONE, None, 1.0)}) free = _free(('k', 'uniform_var', 0.0, 10.0, 0.3), ('sig', 'uniform_var', 0.01, 100.0, 5.0)) obj_formula = LikelihoodObjective(noise=Gaussian(), sigma_sources={'sigma': FormulaSigma('sig')}) @@ -1169,9 +1196,9 @@ def _per_measurement_sigma_setup(sd_base=4.0, s_hi=3.0, exp = _exp_dyn(obs) exp.measurement_params = {'Stot': {'noiseParameter1_Stot': list(tokens)}} routing = ExperimentRouting(routes={ - 'k': ParamRoute('k', PARAM, 'k', 1.0), - 'sd_base': ParamRoute('sd_base', NONE, None, 1.0), - 's_hi': ParamRoute('s_hi', NONE, None, 1.0)}) + 'k': ParamRoute.single('k', PARAM, 'k', 1.0), + 'sd_base': ParamRoute.single('sd_base', NONE, None, 1.0), + 's_hi': ParamRoute.single('s_hi', NONE, None, 1.0)}) free = _free(('k', 'uniform_var', 0.0, 10.0, 0.3), ('sd_base', 'uniform_var', 0.01, 100.0, sd_base), ('s_hi', 'uniform_var', 0.01, 100.0, s_hi)) @@ -1400,7 +1427,7 @@ def test_laplace_scalar_data_fit_gradient(): obj = _laplace_objective() # fixed scale b from the _SD column sim = _sim_with_sensitivities(pred, d_param=dk) exp = _exp(obs, b) - routing = ExperimentRouting(routes={'k': ParamRoute('k', PARAM, 'k', 1.0)}) + routing = ExperimentRouting(routes={'k': ParamRoute.single('k', PARAM, 'k', 1.0)}) free = _free(('k', 'uniform_var', 0.0, 10.0, 0.3)) res = assemble_gaussian_gradient(obj, [(sim, exp, routing)], free) @@ -1420,7 +1447,7 @@ def test_laplace_kink_takes_the_zero_subgradient(): obj = _laplace_objective() sim = _sim_with_sensitivities(pred, d_param=[10.0, -74.0, -110.0, -123.0]) exp = _exp(pred, 2.0) # obs == pred everywhere -> every point a kink - routing = ExperimentRouting(routes={'k': ParamRoute('k', PARAM, 'k', 1.0)}) + routing = ExperimentRouting(routes={'k': ParamRoute.single('k', PARAM, 'k', 1.0)}) free = _free(('k', 'uniform_var', 0.0, 10.0, 0.3)) res = assemble_gaussian_gradient(obj, [(sim, exp, routing)], free) @@ -1445,7 +1472,7 @@ def test_student_t_exact_sqrt_loss_residual(): obj = _student_t_objective() # fixed sigma (_SD), default df=4 sim = _sim_with_sensitivities(pred, d_param=dk) exp = _exp(obs, sigma) - routing = ExperimentRouting(routes={'k': ParamRoute('k', PARAM, 'k', 1.0)}) + routing = ExperimentRouting(routes={'k': ParamRoute.single('k', PARAM, 'k', 1.0)}) free = _free(('k', 'uniform_var', 0.0, 10.0, 0.3)) res = assemble_gaussian_gradient(obj, [(sim, exp, routing)], free) @@ -1479,7 +1506,7 @@ def test_student_t_residual_jacobian_smooth_through_zero(): obj = _student_t_objective() sim = _sim_with_sensitivities(pred, d_param=dk) exp = _exp(obs, sigma) - routing = ExperimentRouting(routes={'k': ParamRoute('k', PARAM, 'k', 1.0)}) + routing = ExperimentRouting(routes={'k': ParamRoute.single('k', PARAM, 'k', 1.0)}) free = _free(('k', 'uniform_var', 0.0, 10.0, 0.3)) res = assemble_gaussian_gradient(obj, [(sim, exp, routing)], free) @@ -1511,9 +1538,9 @@ def test_student_t_estimated_sigma_and_df_columns(): sim = _sim_with_sensitivities(pred, d_param=dk) exp = _exp_dyn(obs) # no _SD: sigma & df are free parameters routing = ExperimentRouting(routes={ - 'k': ParamRoute('k', PARAM, 'k', 1.0), - 'noise_sd': ParamRoute('noise_sd', NONE, None, 1.0), - 'nu_free': ParamRoute('nu_free', NONE, None, 1.0), + 'k': ParamRoute.single('k', PARAM, 'k', 1.0), + 'noise_sd': ParamRoute.single('noise_sd', NONE, None, 1.0), + 'nu_free': ParamRoute.single('nu_free', NONE, None, 1.0), }) free = _free(('k', 'uniform_var', 0.0, 10.0, 0.3), ('noise_sd', 'uniform_var', 0.01, 100.0, sigma), @@ -1548,8 +1575,8 @@ def test_laplace_estimated_scale_column(): sim = _sim_with_sensitivities(pred, d_param=dk) exp = _exp_dyn(obs) routing = ExperimentRouting(routes={ - 'k': ParamRoute('k', PARAM, 'k', 1.0), - 'b_free': ParamRoute('b_free', NONE, None, 1.0), + 'k': ParamRoute.single('k', PARAM, 'k', 1.0), + 'b_free': ParamRoute.single('b_free', NONE, None, 1.0), }) free = _free(('k', 'uniform_var', 0.0, 10.0, 0.3), ('b_free', 'uniform_var', 0.01, 100.0, b)) @@ -1583,7 +1610,7 @@ def test_mixed_gaussian_and_laplace_objective(): obj = LikelihoodObjective( noise=Gaussian(), sigma_sources={'sigma': DataColumnSigma()}, overrides={'B': (Laplace(), {'scale': ConstantSigma(b)})}) - routing = ExperimentRouting(routes={'k': ParamRoute('k', PARAM, 'k', 1.0)}) + routing = ExperimentRouting(routes={'k': ParamRoute.single('k', PARAM, 'k', 1.0)}) free = _free(('k', 'uniform_var', 0.0, 10.0, 0.3)) res = assemble_gaussian_gradient(obj, [(sim, exp, routing)], free) @@ -1612,7 +1639,7 @@ def test_neg_bin_mean_scalar_data_fit_gradient(): obj = _neg_bin_objective(location=MEAN, dispersion=r) sim = _sim_with_sensitivities(pred, d_param=dk) exp = _exp_dyn(obs) # a PMF self-normalizes: no _SD column - routing = ExperimentRouting(routes={'k': ParamRoute('k', PARAM, 'k', 1.0)}) + routing = ExperimentRouting(routes={'k': ParamRoute.single('k', PARAM, 'k', 1.0)}) free = _free(('k', 'uniform_var', 0.0, 10.0, 0.3)) res = assemble_gaussian_gradient(obj, [(sim, exp, routing)], free) @@ -1639,7 +1666,7 @@ def test_neg_bin_median_gradient_matches_prediction_finite_difference(): obj = _neg_bin_objective(location=MEDIAN, dispersion=r) sim = _sim_with_sensitivities(pred, d_param=dk) exp = _exp_dyn(obs) - routing = ExperimentRouting(routes={'k': ParamRoute('k', PARAM, 'k', 1.0)}) + routing = ExperimentRouting(routes={'k': ParamRoute.single('k', PARAM, 'k', 1.0)}) free = _free(('k', 'uniform_var', 0.0, 100.0, 0.3)) res = assemble_gaussian_gradient(obj, [(sim, exp, routing)], free) @@ -1669,8 +1696,8 @@ def test_neg_bin_estimated_dispersion_column(): sim = _sim_with_sensitivities(pred, d_param=dk) exp = _exp_dyn(obs) routing = ExperimentRouting(routes={ - 'k': ParamRoute('k', PARAM, 'k', 1.0), - 'r_free': ParamRoute('r_free', NONE, None, 1.0)}) + 'k': ParamRoute.single('k', PARAM, 'k', 1.0), + 'r_free': ParamRoute.single('r_free', NONE, None, 1.0)}) free = _free(('k', 'uniform_var', 0.0, 10.0, 0.3), ('r_free', 'uniform_var', 0.01, 100.0, r)) @@ -1690,7 +1717,7 @@ def test_capability_gate_now_accepts_laplace_and_student_t(): exact sqrt-loss residual (#459) and so is ``least_squares_exact`` True.""" sim = _sim_with_sensitivities([100, 74, 55, 41], d_param=[0, -74, -110, -123]) exp = _exp([100, 70, 60, 40], 5.0) - routing = ExperimentRouting(routes={'k': ParamRoute('k', PARAM, 'k', 1.0)}) + routing = ExperimentRouting(routes={'k': ParamRoute.single('k', PARAM, 'k', 1.0)}) free = _free(('k', 'uniform_var', 0.0, 10.0, 0.3)) for obj, expect_exact in ((_laplace_objective(), False), (_student_t_objective(), True)): res = assemble_gaussian_gradient(obj, [(sim, exp, routing)], free) @@ -1704,7 +1731,7 @@ def test_capability_gate_now_accepts_negative_binomial(): clause that once refused every negative-binomial (pointing at #458) is gone.""" sim = _sim_with_sensitivities([10, 7, 5, 4], d_param=[0, -7, -11, -12]) exp = _exp_dyn([10, 6, 6, 4]) - routing = ExperimentRouting(routes={'k': ParamRoute('k', PARAM, 'k', 1.0)}) + routing = ExperimentRouting(routes={'k': ParamRoute.single('k', PARAM, 'k', 1.0)}) free = _free(('k', 'uniform_var', 0.0, 10.0, 0.3)) for obj in (_neg_bin_objective(location=MEAN), _neg_bin_objective(location=MEDIAN)): res = assemble_gaussian_gradient(obj, [(sim, exp, routing)], free) @@ -1730,8 +1757,8 @@ def test_capability_gate_now_accepts_median_negative_binomial_with_estimated_dis sim = _sim_with_sensitivities(pred, d_param=dk) exp = _exp_dyn(obs) routing = ExperimentRouting(routes={ - 'k': ParamRoute('k', PARAM, 'k', 1.0), - 'r_free': ParamRoute('r_free', NONE, None, 1.0)}) + 'k': ParamRoute.single('k', PARAM, 'k', 1.0), + 'r_free': ParamRoute.single('r_free', NONE, None, 1.0)}) free = _free(('k', 'uniform_var', 0.0, 100.0, 0.3), ('r_free', 'uniform_var', 0.01, 100.0, r)) @@ -1767,8 +1794,8 @@ def test_capability_gate_now_accepts_mean_on_log_scale_with_estimated_noise(): sim = _sim_with_sensitivities(pred, d_param=dk) exp = _exp_dyn(obs) routing = ExperimentRouting(routes={ - 'k': ParamRoute('k', PARAM, 'k', 1.0), - 'noise_sd': ParamRoute('noise_sd', NONE, None, 1.0), + 'k': ParamRoute.single('k', PARAM, 'k', 1.0), + 'noise_sd': ParamRoute.single('noise_sd', NONE, None, 1.0), }) free = _free(('k', 'uniform_var', 0.0, 10.0, 0.3), ('noise_sd', 'uniform_var', 0.01, 100.0, sigma)) @@ -1794,7 +1821,7 @@ def test_mean_location_on_linear_scale_matches_median(): dk = np.array([0.0, -74.0, -110.0, -123.0]) sim = _sim_with_sensitivities(pred, d_param=dk) exp = _exp(obs, sigma) - routing = ExperimentRouting(routes={'k': ParamRoute('k', PARAM, 'k', 1.0)}) + routing = ExperimentRouting(routes={'k': ParamRoute.single('k', PARAM, 'k', 1.0)}) free = _free(('k', 'uniform_var', 0.0, 10.0, 0.3)) mean_obj = LikelihoodObjective(noise=Gaussian(location=MEAN), @@ -1815,7 +1842,7 @@ def test_routed_key_absent_from_tensor_refuses(): obj = ChiSquareObjective() sim = _sim_with_sensitivities([100, 74, 55, 41], d_param=[0, -74, -110, -123]) # only 'k' exp = _exp([100, 70, 60, 40], 5.0) - routing = ExperimentRouting(routes={'k2': ParamRoute('k2', PARAM, 'k2', 1.0)}) # 'k2' absent + routing = ExperimentRouting(routes={'k2': ParamRoute.single('k2', PARAM, 'k2', 1.0)}) # 'k2' absent free = _free(('k2', 'uniform_var', 0.0, 10.0, 0.3)) with pytest.raises(GradientNotSupported): assemble_gaussian_gradient(obj, [(sim, exp, routing)], free) @@ -1827,7 +1854,7 @@ def test_missing_sensitivity_payload_refuses(): obj = ChiSquareObjective() sim = Data.from_columns(np.column_stack([TIMES, [100, 74, 55, 41]]), ['time', 'Stot']) exp = _exp([100, 70, 60, 40], 5.0) - routing = ExperimentRouting(routes={'k': ParamRoute('k', PARAM, 'k', 1.0)}) + routing = ExperimentRouting(routes={'k': ParamRoute.single('k', PARAM, 'k', 1.0)}) free = _free(('k', 'uniform_var', 0.0, 10.0, 0.3)) with pytest.raises(GradientNotSupported): assemble_gaussian_gradient(obj, [(sim, exp, routing)], free) @@ -1864,8 +1891,8 @@ def test_log_scale_multiplies_the_right_jacobian_column(): sim = _sim_with_sensitivities(pred, d_param=dk, d_ic=ds0) exp = _exp(obs, sigma) routing = ExperimentRouting(routes={ - 'k': ParamRoute('k', PARAM, 'k', 1.0), - 'S0': ParamRoute('S0', IC, 'S()', 1.0), + 'k': ParamRoute.single('k', PARAM, 'k', 1.0), + 'S0': ParamRoute.single('S0', IC, 'S()', 1.0), }) free = _free(('k', 'loguniform_var', 0.01, 100.0, 0.3), # log10 ('S0', 'uniform_var', 0.0, 500.0, 120.0)) # linear @@ -2053,8 +2080,8 @@ def test_log_scaled_free_sigma_column_takes_the_sampling_factor(): sim = _sim_with_sensitivities(pred, d_param=[0.0, -74.0, -110.0, -123.0]) exp = _exp_dyn(obs) routing = ExperimentRouting(routes={ - 'k': ParamRoute('k', PARAM, 'k', 1.0), - 'noise_sd': ParamRoute('noise_sd', NONE, None, 1.0), + 'k': ParamRoute.single('k', PARAM, 'k', 1.0), + 'noise_sd': ParamRoute.single('noise_sd', NONE, None, 1.0), }) free = _free(('k', 'uniform_var', 0.0, 10.0, 0.3), ('noise_sd', 'loguniform_var', 0.01, 100.0, sigma)) @@ -2693,7 +2720,7 @@ def _constraint_sim(stot, dk, ds0, model='m', suffix='tc'): d_ic=np.asarray(ds0, float).reshape(len(stot), 1, 1)) sdd = {model: {suffix: sim}} routings = {(model, suffix): ExperimentRouting(routes={ - 'k': ParamRoute('k', PARAM, 'k', 1.0), 'S0': ParamRoute('S0', IC, 'S()', 1.0)})} + 'k': ParamRoute.single('k', PARAM, 'k', 1.0), 'S0': ParamRoute.single('S0', IC, 'S()', 1.0)})} free = _free(('k', 'uniform_var', 0.0, 100.0, 0.4), ('S0', 'uniform_var', 0.0, 1000.0, 120.0)) return sdd, routings, free @@ -2717,6 +2744,56 @@ def test_constraint_static_penalty_gradient(): np.testing.assert_allclose(g, [-2.0 * _C_DK[2], -2.0 * _C_DS0[2]]) +def test_constraint_gradient_sums_a_multi_contribution_route(): + """The **constraint** accessor sums a route's contributions too (ADR-0076, #511) -- the + constraint peer of ``test_multi_contribution_route_sums_columns``. One free parameter 'p' + reaches both the parameter axis (k) and the IC axis (S()), so the constraint readout's + sensitivity is their sum.""" + times = np.arange(len(_C_STOT), dtype=float) + sim = Data.from_columns(np.column_stack([times, np.asarray(_C_STOT, float)]), + ['time', 'Stot']) + sim.output_sensitivities = OutputSensitivities( + selectors=['observable:Stot'], param_names=['k'], ic_species=['S()'], + d_param=np.asarray(_C_DK, float).reshape(len(_C_STOT), 1, 1), + d_ic=np.asarray(_C_DS0, float).reshape(len(_C_STOT), 1, 1)) + sdd = {'m': {'tc': sim}} + routings = {('m', 'tc'): ExperimentRouting(routes={'p': ParamRoute('p', ( + RouteContribution(PARAM, 'k', 1.0), RouteContribution(IC, 'S()', 1.0)))})} + free = _free(('p', 'uniform_var', 0.0, 100.0, 0.4)) + + c = AtConstraint('Stot', '>', 90.0, 'm', 'tc', weight=2.0, atvar=None, atval=2.0) + cset = ConstraintSet('m', 'tc') + cset.constraints = [c] + + g = assemble_constraint_gradient([cset], sdd, routings, free) + np.testing.assert_allclose(g, [-2.0 * (_C_DK[2] + _C_DS0[2])]) + + +def test_fisher_hessian_with_a_multi_contribution_route_equals_jtj(): + """The Fisher/EFIM block (gntr's curvature) is built from the same summed Jacobian column as + the gradient for a multi-contribution route (ADR-0076, #511), so it still equals ``J^T J`` + and stays symmetric PSD -- the curvature peer of the gradient sum test.""" + pred = np.array([120.0, 80.0, 53.0, 36.0]) + obs = np.array([118.0, 82.0, 50.0, 38.0]) + sigma = 4.0 + dk = np.array([0.0, -80.0, -106.0, -108.0]) + ds0 = np.array([1.0, 0.66, 0.44, 0.30]) + obj = ChiSquareObjective() + sim = _sim_with_sensitivities(pred, d_param=dk, d_ic=ds0) + exp = _exp(obs, sigma) + routing = ExperimentRouting(routes={'m': ParamRoute('m', ( + RouteContribution(PARAM, 'k', 1.0), RouteContribution(IC, 'S()', 1.0)))}) + free = _free(('m', 'uniform_var', 0.0, 10.0, 0.3)) + + res = assemble_gaussian_gradient(obj, [(sim, exp, routing)], free) + H = assemble_fisher_hessian(obj, [(sim, exp, routing)], free) + + np.testing.assert_allclose(res.jacobian, ((1.0 / sigma) * (dk + ds0)).reshape(-1, 1)) + np.testing.assert_allclose(H, res.jacobian.T @ res.jacobian) + np.testing.assert_allclose(H, H.T) + assert np.all(np.linalg.eigvalsh(H) >= -1e-9) + + def test_constraint_satisfied_has_zero_gradient(): """A satisfied constraint contributes no penalty and no gradient (the penalty is flat 0 in the satisfied region; the boundary kink takes the subgradient 0).""" @@ -3392,7 +3469,7 @@ def test_combined_gradient_fisher_assembles_each_scored_sensitivity_once(monkeyp sim = _sim_with_sensitivities(pred, d_param=dk, d_ic=ds0) exp = _exp(obs, 4.0) routing = ExperimentRouting(routes={ - 'k': ParamRoute('k', PARAM, 'k', 4.0), 'S0': ParamRoute('S0', IC, 'S()', 1.0)}) + 'k': ParamRoute.single('k', PARAM, 'k', 4.0), 'S0': ParamRoute.single('S0', IC, 'S()', 1.0)}) free = _free(('k', 'uniform_var', 0.0, 10.0, 0.3), ('S0', 'uniform_var', 0.0, 500.0, 120.0)) experiments = [(sim, exp, routing)] @@ -3440,7 +3517,7 @@ def test_fisher_hessian_gaussian_equals_jtj(): sim = _sim_with_sensitivities(pred, d_param=dk, d_ic=ds0) exp = _exp(obs, sigma) routing = ExperimentRouting(routes={ - 'k': ParamRoute('k', PARAM, 'k', 4.0), 'S0': ParamRoute('S0', IC, 'S()', 1.0)}) + 'k': ParamRoute.single('k', PARAM, 'k', 4.0), 'S0': ParamRoute.single('S0', IC, 'S()', 1.0)}) free = _free(('k', 'uniform_var', 0.0, 10.0, 0.3), ('S0', 'uniform_var', 0.0, 500.0, 120.0)) res = assemble_gaussian_gradient(obj, [(sim, exp, routing)], free) @@ -3458,10 +3535,10 @@ def test_fisher_hessian_sums_across_experiments(): free = _free(('k', 'uniform_var', 0.0, 10.0, 0.3)) sim1 = _sim_with_sensitivities([100, 74, 55, 41], d_param=[0, -74, -110, -123]) exp1 = _exp([100, 70, 60, 40], 5.0) - r1 = ExperimentRouting(routes={'k': ParamRoute('k', PARAM, 'k', 1.0)}) + r1 = ExperimentRouting(routes={'k': ParamRoute.single('k', PARAM, 'k', 1.0)}) sim2 = _sim_with_sensitivities([100, 55, 30, 17], d_param=[0, -110, -120, -100]) exp2 = _exp([100, 58, 28, 20], 3.0) - r2 = ExperimentRouting(routes={'k': ParamRoute('k', PARAM, 'k', 4.0)}) + r2 = ExperimentRouting(routes={'k': ParamRoute.single('k', PARAM, 'k', 4.0)}) both = assemble_fisher_hessian(obj, [(sim1, exp1, r1), (sim2, exp2, r2)], free) h1 = assemble_fisher_hessian(obj, [(sim1, exp1, r1)], free) @@ -3482,8 +3559,8 @@ def test_fisher_hessian_estimated_sigma_adds_diagonal_noise_block(): sim = _sim_with_sensitivities(pred, d_param=dk) exp = _exp_dyn(obs) routing = ExperimentRouting(routes={ - 'k': ParamRoute('k', PARAM, 'k', 1.0), - 'noise_sd': ParamRoute('noise_sd', NONE, None, 1.0)}) + 'k': ParamRoute.single('k', PARAM, 'k', 1.0), + 'noise_sd': ParamRoute.single('noise_sd', NONE, None, 1.0)}) free = _free(('k', 'uniform_var', 0.0, 10.0, 0.3), ('noise_sd', 'uniform_var', 0.01, 100.0, sigma)) @@ -3506,7 +3583,7 @@ def test_fisher_hessian_laplace_location_block(): obj = _laplace_objective() # fixed _SD column scale sim = _sim_with_sensitivities(pred, d_param=dk) exp = _exp(obs, b) # Stot_SD column == b - routing = ExperimentRouting(routes={'k': ParamRoute('k', PARAM, 'k', 1.0)}) + routing = ExperimentRouting(routes={'k': ParamRoute.single('k', PARAM, 'k', 1.0)}) free = _free(('k', 'uniform_var', 0.0, 10.0, 0.3)) H = assemble_fisher_hessian(obj, [(sim, exp, routing)], free) @@ -3524,7 +3601,7 @@ def test_fisher_hessian_neg_bin_mean_location_block(): obj = _neg_bin_objective(location=MEAN, dispersion=r) sim = _sim_with_sensitivities(pred, d_param=dk) exp = _exp_dyn(obs) # self-normalizing PMF, no _SD column - routing = ExperimentRouting(routes={'k': ParamRoute('k', PARAM, 'k', 1.0)}) + routing = ExperimentRouting(routes={'k': ParamRoute.single('k', PARAM, 'k', 1.0)}) free = _free(('k', 'uniform_var', 0.0, 100.0, 0.3)) H = assemble_fisher_hessian(obj, [(sim, exp, routing)], free) @@ -3538,7 +3615,7 @@ def test_fisher_hessian_refuses_median_negative_binomial(): obj = _neg_bin_objective(location=MEDIAN, dispersion=6.0) sim = _sim_with_sensitivities([50.0, 40.0, 30.0, 22.0], d_param=[0.0, -8.0, -12.0, -11.0]) exp = _exp_dyn([48.0, 42.0, 28.0, 24.0]) - routing = ExperimentRouting(routes={'k': ParamRoute('k', PARAM, 'k', 1.0)}) + routing = ExperimentRouting(routes={'k': ParamRoute.single('k', PARAM, 'k', 1.0)}) free = _free(('k', 'uniform_var', 0.0, 100.0, 0.3)) with pytest.raises(GradientNotSupported): assemble_fisher_hessian(obj, [(sim, exp, routing)], free) @@ -3552,8 +3629,8 @@ def test_fisher_hessian_refuses_mean_on_log_estimated_scale(): sim = _sim_with_sensitivities([100.0, 74.0, 55.0, 41.0], d_param=[0.0, -74.0, -110.0, -123.0]) exp = _exp_dyn([100.0, 70.0, 60.0, 40.0]) routing = ExperimentRouting(routes={ - 'k': ParamRoute('k', PARAM, 'k', 1.0), - 'noise_sd': ParamRoute('noise_sd', NONE, None, 1.0)}) + 'k': ParamRoute.single('k', PARAM, 'k', 1.0), + 'noise_sd': ParamRoute.single('noise_sd', NONE, None, 1.0)}) free = _free(('k', 'uniform_var', 0.0, 10.0, 0.3), ('noise_sd', 'uniform_var', 0.01, 100.0, 0.2)) with pytest.raises(GradientNotSupported): diff --git a/tests/test_gradient_routing.py b/tests/test_gradient_routing.py index d13f3086..e79780ab 100644 --- a/tests/test_gradient_routing.py +++ b/tests/test_gradient_routing.py @@ -19,13 +19,18 @@ from pybnf.gradient import routing as R from pybnf.gradient.routing import ( - PARAM, IC, NONE, ParamRoute, - classify_free_param, condition_factor, route_experiment, + PARAM, IC, NONE, RouteContribution, ParamRoute, + classify_free_param, classify_condition_target, condition_factor, route_experiment, ) from pybnf.pset import Mutation, MutationSet from pybnf.printing import PybnfError +def _pref(*mutations): + """A MutationSet from ``(name, op, value)`` triples, each a parameter-reference (ADR-0076).""" + return MutationSet([Mutation(n, op, v, is_param_ref=True) for n, op, v in mutations], 'c') + + FIXTURES = Path(__file__).resolve().parent / 'bngl_files' # The committed analytic-decay fixture's bind-by-id namespaces: parameters S0 (initial value) @@ -107,9 +112,11 @@ def test_route_experiment_wildtype(): r = route_experiment(FREE, DECAY_PARAMS, DECAY_SPECIES, None) assert r.sensitivity_params == ['k'] assert r.sensitivity_ic == ['S()'] - assert r.routes['k'] == ParamRoute('k', PARAM, 'k', 1.0) - assert r.routes['S0'] == ParamRoute('S0', IC, 'S()', 1.0) - assert r.routes['sigma'] == ParamRoute('sigma', NONE, None, 1.0) + assert r.routes['k'] == ParamRoute('k', (RouteContribution(PARAM, 'k', 1.0),)) + assert r.routes['S0'] == ParamRoute('S0', (RouteContribution(IC, 'S()', 1.0),)) + assert r.routes['sigma'] == ParamRoute('sigma', (RouteContribution(NONE, None, 1.0),)) + # The single-contribution convenience accessors read that sole contribution. + assert (r.routes['k'].target, r.routes['k'].key, r.routes['k'].factor) == (PARAM, 'k', 1.0) def test_route_experiment_pinned_param_drops_request_column(): @@ -148,15 +155,166 @@ def test_route_experiment_preserves_declaration_order_and_dedups(): assert r.sensitivity_params == ['k2', 'k1'] -def test_route_experiment_parameter_reference_condition_raises(): - """A per-condition estimated initial condition (ADR-0076) -- a condition that sets a model - entity to the value of a free parameter -- refuses on the gradient path rather than - silently emitting a zero sensitivity column: the referenced parameter reaches the - trajectory through a *different* id, which bind-by-id routing does not model.""" +# --------------------------- per-condition estimated initial conditions (ADR-0076, #511) ---- + +# ic_seed_map: the model parameter S0 bares species S()'s initial value (a bare +# initialAssignment / .net initializer), so a condition that sets S0 to a free parameter routes +# that free parameter onto the S() initial-condition sensitivity axis. +IC_SEED_MAP = {'S0': 'S()'} + + +def test_classify_condition_target_ic_seed(): + """A condition target that bares a species IC routes to the IC axis, factor 1.""" + assert classify_condition_target('S0', {'S0', 'k'}, {'S()'}, IC_SEED_MAP) == (IC, 'S()', 1.0) + + +def test_classify_condition_target_param(): + """A condition target that is an ordinary global routes to the parameter axis, factor 1.""" + assert classify_condition_target('k', {'S0', 'k'}, {'S()'}, IC_SEED_MAP) == (PARAM, 'k', 1.0) + + +def test_classify_condition_target_non_bare_seed_refuses(): + """A parameter that seeds a species IC non-baruely (map value None) refuses rather than + emitting a parameter-dependent factor.""" from pybnf.gradient import GradientNotSupported - cond = MutationSet([Mutation('S0', '=', 'S0_A', is_param_ref=True)], 'c') - with pytest.raises(GradientNotSupported, match='S0_A'): - route_experiment(['k', 'S0_A'], DECAY_PARAMS, DECAY_SPECIES, cond) + with pytest.raises(GradientNotSupported, match='non-bare'): + classify_condition_target('S0', {'S0', 'k'}, {'S()'}, {'S0': None}) + + +def test_classify_condition_target_unbindable_refuses(): + """A target that is neither a parameter nor a species IC binds no sensitivity column.""" + from pybnf.gradient import GradientNotSupported + with pytest.raises(GradientNotSupported, match='cannot route'): + classify_condition_target('nope', {'S0', 'k'}, {'S()'}, IC_SEED_MAP) + + +def test_route_experiment_param_ref_routes_to_ic(): + """A per-condition estimated initial condition ``S0 = S0_A`` routes the *referenced* free + parameter S0_A onto species S()'s IC axis (chain-rule factor 1), instead of aborting.""" + r = route_experiment(['k', 'S0_A'], DECAY_PARAMS, DECAY_SPECIES, + _pref(('S0', '=', 'S0_A')), ic_seed_map=IC_SEED_MAP) + assert r.routes['S0_A'].contributions == (RouteContribution(IC, 'S()', 1.0),) + assert r.sensitivity_ic == ['S()'] + assert r.sensitivity_params == ['k'] # k still binds by id + + +def test_route_experiment_param_ref_routes_to_param(): + """A condition ``k = kfree`` (target is an ordinary global) routes kfree onto the parameter + axis for k.""" + r = route_experiment(['kfree'], DECAY_PARAMS, DECAY_SPECIES, + _pref(('k', '=', 'kfree')), ic_seed_map=IC_SEED_MAP) + assert r.routes['kfree'].contributions == (RouteContribution(PARAM, 'k', 1.0),) + assert r.sensitivity_params == ['k'] + + +def test_route_experiment_param_ref_multi_target_sums(): + """One free parameter assigned to several targets in one condition (a shared multiplier) + accumulates one contribution per target -- its derivative is their sum.""" + r = route_experiment(['m'], DECAY_PARAMS, DECAY_SPECIES, + _pref(('k', '=', 'm'), ('S0', '=', 'm')), ic_seed_map=IC_SEED_MAP) + assert r.routes['m'].contributions == ( + RouteContribution(PARAM, 'k', 1.0), RouteContribution(IC, 'S()', 1.0)) + assert r.sensitivity_params == ['k'] + assert r.sensitivity_ic == ['S()'] + # A multi-contribution route has no single target/key/factor. + with pytest.raises(ValueError): + _ = r.routes['m'].target + + +def test_route_experiment_param_ref_non_bare_seed_refuses(): + """A per-condition estimated IC through a non-bare seed (map value None) refuses.""" + from pybnf.gradient import GradientNotSupported + with pytest.raises(GradientNotSupported, match='non-bare'): + route_experiment(['S0_A'], DECAY_PARAMS, DECAY_SPECIES, + _pref(('S0', '=', 'S0_A')), ic_seed_map={'S0': None}) + + +def test_route_experiment_param_ref_non_equals_refuses(): + """A parameter reference must be an ``=`` assignment (ADR-0076); a relative op refuses.""" + from pybnf.gradient import GradientNotSupported + with pytest.raises(GradientNotSupported, match="non-'='"): + route_experiment(['m'], DECAY_PARAMS, DECAY_SPECIES, + _pref(('k', '*', 'm')), ic_seed_map=IC_SEED_MAP) + + +def test_route_experiment_param_ref_composes_with_a_non_unit_condition_factor(): + """A free parameter that is BOTH scaled by its own condition perturbation AND param-refed + onto another target keeps each term's own chain-rule factor: the base bind carries the + condition factor (``k*3`` -> 3), the param-ref term carries ``d(target)/d(k) = 1``. Their + sum is the derivative (#511).""" + cond = MutationSet([Mutation('k', '*', 3.0), # scales the base bind + Mutation('S0', '=', 'k', is_param_ref=True)], # k also seeds S()'s IC + 'c') + r = route_experiment(['k'], DECAY_PARAMS, DECAY_SPECIES, cond, ic_seed_map=IC_SEED_MAP) + assert r.routes['k'].contributions == ( + RouteContribution(PARAM, 'k', 3.0), RouteContribution(IC, 'S()', 1.0)) + assert r.sensitivity_params == ['k'] and r.sensitivity_ic == ['S()'] + + +def test_route_experiment_param_ref_survives_a_pinned_base_bind(): + """A free parameter pinned out of its own base bind (factor 0) still contributes through the + condition's param-ref: the zero term drops from the request, the param-ref term does not.""" + cond = MutationSet([Mutation('k', '=', 0.5), # base bind pinned -> 0 + Mutation('S0', '=', 'k', is_param_ref=True)], + 'c') + r = route_experiment(['k'], DECAY_PARAMS, DECAY_SPECIES, cond, ic_seed_map=IC_SEED_MAP) + assert r.routes['k'].contributions == ( + RouteContribution(PARAM, 'k', 0.0), RouteContribution(IC, 'S()', 1.0)) + assert r.sensitivity_params == [] # pinned term dropped from the request + assert r.sensitivity_ic == ['S()'] # param-ref term survives + + +def test_route_experiment_param_ref_directly_bound_free_param_still_binds(): + """A free parameter that binds by id AND is param-referenced accumulates both contributions + (its base bind plus the condition target).""" + r = route_experiment(['k'], DECAY_PARAMS, DECAY_SPECIES, + _pref(('S0', '=', 'k')), ic_seed_map=IC_SEED_MAP) + # k binds param 'k' by id (base) and is routed to IC 'S()' by the condition. + assert r.routes['k'].contributions == ( + RouteContribution(PARAM, 'k', 1.0), RouteContribution(IC, 'S()', 1.0)) + + +class _CapturingModel: + """Records the request :func:`apply_routings` hands to the gradient path.""" + + def __init__(self): + self.applied = None + + def enable_output_sensitivities(self, *, params=None, ic=None): + self.applied = (list(params or []), list(ic or [])) + + +def test_apply_routings_unions_a_column_reached_only_through_a_condition(): + """The applied request is the UNION over routings, not the wildtype's. + + A free parameter routed only through a condition (a per-condition estimated initial + condition, ADR-0076) reaches a column the wildtype never binds, so the wildtype request is + NOT a superset -- the union must carry it or the assembly aborts on a missing column + (#511).""" + free = ['k', 's0_free'] + wildtype = route_experiment(free, DECAY_PARAMS, DECAY_SPECIES, None, ic_seed_map=IC_SEED_MAP) + conditioned = route_experiment(free, DECAY_PARAMS, DECAY_SPECIES, + _pref(('S0', '=', 's0_free')), ic_seed_map=IC_SEED_MAP) + # The wildtype binds no IC column at all: s0_free matches no model id on its own. + assert wildtype.sensitivity_ic == [] + assert conditioned.sensitivity_ic == ['S()'] + + model = _CapturingModel() + params, ic = R.apply_routings(model, [wildtype, conditioned]) + + assert ic == ['S()'] # carried in from the condition alone + assert params == ['k'] + assert model.applied == (['k'], ['S()']) + + +def test_apply_routings_dedups_across_routings(): + """A column several conditions reach is requested once.""" + free = ['k', 'S0'] + r1 = route_experiment(free, DECAY_PARAMS, DECAY_SPECIES, None) + r2 = route_experiment(free, DECAY_PARAMS, DECAY_SPECIES, None) + model = _CapturingModel() + params, ic = R.apply_routings(model, [r1, r2]) + assert params == ['k'] and ic == ['S()'] # ------------------------------------------------- model adapter (bngsim) ---- @@ -177,9 +335,12 @@ def decay_model(): @pytest.mark.bngsim def test_sensitivity_entity_namespace(decay_model): - param_ids, species = decay_model.sensitivity_entity_namespace() + param_ids, species, ic_seed_map = decay_model.sensitivity_entity_namespace() assert set(param_ids) == {'S0', 'k'} assert species == [('S()', 'S0')] + # S0 bares species S()'s initial value, so a condition setting S0 to a free parameter can + # route that free parameter onto the S() IC axis (ADR-0076, #511). + assert ic_seed_map == {'S0': 'S()'} @pytest.mark.bngsim @@ -190,6 +351,31 @@ def test_route_for_model_matches_pure_core(decay_model): assert r.routes['S0'].target == IC and r.routes['S0'].key == 'S()' +@pytest.mark.bngsim +def test_route_for_model_composes_param_ref_on_the_net_backend(decay_model): + """A per-condition estimated initial condition composes on the **net** backend too + (ADR-0076, #511): the .net species block's bare initializer ``S() <- S0`` makes S0 an IC + seed, so a condition setting S0 to free parameter S0_A routes S0_A onto species S()'s IC + axis -- the net peer of the SBML initialAssignment path.""" + cond = MutationSet([Mutation('S0', '=', 'S0_A', is_param_ref=True)], 'c') + r = R.route_for_model(decay_model, ['k', 'S0_A'], cond) + assert r.routes['S0_A'].contributions == (RouteContribution(IC, 'S()', 1.0),) + assert r.sensitivity_ic == ['S()'] + assert r.sensitivity_params == ['k'] + + +@pytest.mark.bngsim +def test_route_for_model_net_backend_multi_target_param_ref_sums(decay_model): + """One free parameter a condition assigns to both a rate parameter and an IC-seeding + parameter accumulates both contributions on the net backend (their sum, #511).""" + cond = MutationSet([Mutation('k', '=', 'm', is_param_ref=True), + Mutation('S0', '=', 'm', is_param_ref=True)], 'c') + r = R.route_for_model(decay_model, ['m'], cond) + assert r.routes['m'].contributions == ( + RouteContribution(PARAM, 'k', 1.0), RouteContribution(IC, 'S()', 1.0)) + assert r.sensitivity_params == ['k'] and r.sensitivity_ic == ['S()'] + + @pytest.mark.bngsim def test_route_for_model_resolves_condition_by_name(decay_model): decay_model.add_mutant(MutationSet([Mutation('k', '*', 4.0)], 'hi')) diff --git a/tests/test_gradient_sbml.py b/tests/test_gradient_sbml.py index 3913ddf3..ca23f9d0 100644 --- a/tests/test_gradient_sbml.py +++ b/tests/test_gradient_sbml.py @@ -41,7 +41,7 @@ import pybnf.bngsim_sbml_model as bngsim_sbml_model from pybnf._bngsim_caps import BNGSIM_HAS_OUTPUT_SENS from pybnf.data import Data -from pybnf.gradient import assemble_gaussian_gradient, route_experiment +from pybnf.gradient import assemble_gaussian_gradient, route_experiment, route_for_model from pybnf.objective import ChiSquareObjective from pybnf.printing import PybnfError from pybnf.pset import ( @@ -332,6 +332,217 @@ def loss_at(u_vec): np.testing.assert_allclose(res.gradient, grad_fd, rtol=1e-3, atol=1e-3) +# --- #511: per-condition estimated initial conditions through a condition ------ # +# A mini-Bruno: species S's initial value is seeded by a bare initialAssignment from global +# param S0, and two decay channels carry rate multipliers kmult1/kmult2. A condition sets S0 +# and BOTH multipliers to the value of free parameters (a per-condition estimated initial +# condition + a shared multiplier, ADR-0076), so the free parameters reach the model ONLY +# through the condition -- S(t) = S0*exp(-(k1*kmult1 + k2*kmult2)*t). +_SEEDED_SBML = """ + + + + + + + + + + + + + + + S0 + + + + + k1kmult1S + + + + k2kmult2S + + + +""" + + +@_needs_output_sens +def test_sbml_ic_seed_map_exposes_bare_initial_assignment(tmp_path): + """The SBML namespace exposes ``{S0 -> S}`` for the bare ``initialAssignment`` ``S = S0``, + so the router can compose a per-condition estimated initial condition (ADR-0076, #511).""" + xml = Path(tmp_path) / 'seeded.xml' + xml.write_text(_SEEDED_SBML) + ps = PSet([FreeParameter('k1', 'uniform_var', 0.0, 1e6, value=0.3)]) + model = bngsim_sbml_model.BngsimSbmlModelNoTimeout( + str(xml), str(xml), pset=ps, actions=(TimeCourse({'time': '10', 'step': '1'}),)) + param_ids, species, ic_seed_map = model.sensitivity_entity_namespace() + assert set(param_ids) == {'S0', 'k1', 'k2', 'kmult1', 'kmult2'} + assert ic_seed_map == {'S0': 'S'} + + +@_needs_output_sens +def test_sbml_fd_oracle_free_params_routed_through_a_condition(tmp_path): + """Central-difference FD of loss(u) vs the assembled gradient(u) when the free parameters + reach the model ONLY through a condition (ADR-0076, #511): ``s0_free`` sets the IC-seeding + param S0 (a per-condition estimated initial condition -> IC axis), and ``m_free`` sets BOTH + rate multipliers at once (a shared multiplier -> the SUM of two parameter-axis columns). + ``k1`` binds by id. The pre-fix gradient path aborted on this; here its gradient must match + finite differences on both composed axes.""" + xml_path = Path(tmp_path) / 'seeded.xml' + xml_path.write_text(_SEEDED_SBML) + xml = str(xml_path) + + def run(s0, m, k1, with_sensitivities, route=None): + # The condition sets S0=s0 (seeds S's IC), kmult1=kmult2=m (both channels), k1=k1. + ps = PSet([FreeParameter('S0', 'uniform_var', 0.0, 1e6, value=s0), + FreeParameter('kmult1', 'uniform_var', 0.0, 1e6, value=m), + FreeParameter('kmult2', 'uniform_var', 0.0, 1e6, value=m), + FreeParameter('k1', 'uniform_var', 0.0, 1e6, value=k1)]) + model = bngsim_sbml_model.BngsimSbmlModelNoTimeout( + xml, xml, pset=ps, actions=(TimeCourse({'time': '10', 'step': '1'}),)) + if with_sensitivities: + model.enable_output_sensitivities( + params=route.sensitivity_params, ic=route.sensitivity_ic) + return model.execute(str(tmp_path), 'fd511', 0)['time_course'] + + obj = ChiSquareObjective() + free = [FreeParameter('s0_free', 'uniform_var', 0.0, 1000.0, value=120.0), + FreeParameter('m_free', 'uniform_var', 0.01, 100.0, value=0.8), + FreeParameter('k1', 'uniform_var', 0.01, 100.0, value=0.4)] + names = [p.name for p in free] + + # The routing comes from the live model's namespaces (the SBML ic_seed_map end-to-end): the + # condition param-refs S0/kmult1/kmult2 to free parameters. + cond = MutationSet([ + Mutation('S0', '=', 's0_free', is_param_ref=True), + Mutation('kmult1', '=', 'm_free', is_param_ref=True), + Mutation('kmult2', '=', 'm_free', is_param_ref=True), + ], 'c') + route_model = bngsim_sbml_model.BngsimSbmlModelNoTimeout( + xml, xml, pset=PSet([FreeParameter('k1', 'uniform_var', 0.0, 1e6, value=0.3)]), + actions=(TimeCourse({'time': '10', 'step': '1'}),)) + route = route_for_model(route_model, names, cond) + # s0_free -> IC(S); m_free -> PARAM(kmult1) + PARAM(kmult2) (summed); k1 -> PARAM(k1). + assert route.sensitivity_ic == ['S'] + assert set(route.sensitivity_params) == {'kmult1', 'kmult2', 'k1'} + assert len(route.routes['m_free'].contributions) == 2 + + # Synthetic data at the true params -> non-zero residuals at the evaluation point. + s0_true, m_true, k1_true, sigma = 100.0, 1.0, 0.3, 5.0 + exp = _exp_from_species(run(s0_true, m_true, k1_true, False), sigma) + + def loss_at(u_vec): + theta = {n: p.from_sampling_space(u) for n, p, u in zip(names, free, u_vec)} + return obj.evaluate(run(theta['s0_free'], theta['m_free'], theta['k1'], False), exp) + + u0 = np.array([p.to_sampling_space(p.value) for p in free]) + h = 1e-6 + grad_fd = np.zeros(len(free)) + for j in range(len(free)): + up, um = u0.copy(), u0.copy() + up[j] += h + um[j] -= h + grad_fd[j] = (loss_at(up) - loss_at(um)) / (2.0 * h) + + sim = run(free[0].value, free[1].value, free[2].value, True, route=route) + res = assemble_gaussian_gradient(obj, [(sim, exp, route)], free) + + assert res.param_names == names + np.testing.assert_allclose(res.gradient, grad_fd, rtol=2e-3, atol=2e-3) + + +_NONBARE_SEED_SBML = """ + + + + + + + + + + + + 2S0 + + + + + kS + + + +""" + + +@_needs_output_sens +def test_sbml_non_bare_initial_assignment_seed_is_refused(tmp_path): + """A **non-bare** initialAssignment ``S = 2*S0`` has a parameter-dependent d(IC)/d(S0), so + S0 is exposed as non-routable (``{'S0': None}``) and routing a per-condition estimated + initial condition through it refuses rather than emitting a wrong factor (ADR-0076, #511).""" + from pybnf.gradient import GradientNotSupported + xml_path = Path(tmp_path) / 'nonbare.xml' + xml_path.write_text(_NONBARE_SEED_SBML) + xml = str(xml_path) + ps = PSet([FreeParameter('k', 'uniform_var', 0.0, 1e6, value=0.3)]) + model = bngsim_sbml_model.BngsimSbmlModelNoTimeout( + xml, xml, pset=ps, actions=(TimeCourse({'time': '10', 'step': '1'}),)) + _, _, ic_seed_map = model.sensitivity_entity_namespace() + assert ic_seed_map == {'S0': None} + cond = MutationSet([Mutation('S0', '=', 's0_free', is_param_ref=True)], 'c') + with pytest.raises(GradientNotSupported, match='non-bare'): + route_for_model(model, ['s0_free'], cond) + + +_AMOUNT_SEED_SBML = """ + + + + + + + + + + + + S0 + + + + + kS + + + +""" + + +@_needs_output_sens +def test_sbml_amount_species_with_non_unit_volume_seed_is_refused(tmp_path): + """A **bare** initialAssignment is still non-routable when the species is amount-based in a + non-unit compartment: the PyBNF species value is an amount but bngsim's IC axis is a + concentration, so d(IC)/d(S0) is 1/size, not 1. The seed is exposed as non-routable and the + per-condition estimated initial condition refuses (ADR-0076, #511).""" + from pybnf.gradient import GradientNotSupported + xml_path = Path(tmp_path) / 'amount.xml' + xml_path.write_text(_AMOUNT_SEED_SBML) + xml = str(xml_path) + ps = PSet([FreeParameter('k', 'uniform_var', 0.0, 1e6, value=0.3)]) + model = bngsim_sbml_model.BngsimSbmlModelNoTimeout( + xml, xml, pset=ps, actions=(TimeCourse({'time': '10', 'step': '1'}),)) + # size-2 compartment on an amount species -> unit factor 0.5, so S0 is NOT a unit-slope seed. + assert model._species_unit_factor['S'] == 0.5 + _, _, ic_seed_map = model.sensitivity_entity_namespace() + assert ic_seed_map == {'S0': None} + cond = MutationSet([Mutation('S0', '=', 's0_free', is_param_ref=True)], 'c') + with pytest.raises(GradientNotSupported, match='not a plain 1'): + route_for_model(model, ['s0_free'], cond) + + # --- setup regression guard (default CI, no full fit) -------------------------- # def _decay_recovery_config(tmp_path, fit_type, *, observables=None, obs_column='S', obs_formula=None, **overrides):