Constraint gradient/Hessian use the raw sensitivity tensor while the penalty is scored on the normalized column, so the ADR-0053 chain rule is dropped

[wshlavacek]

What happens

_constraint_sensitivity_accessor's raw_sens body (assembly.py:1082-1108) goes straight from _selector_for(sens, observable) to _sensitivity(sens, selector, c, row, name) on the #447 forward tensor. It never reads sim_data.normalization, has no _normalized_sensitivity wrapper, and raises no GradientNotSupported when the column it is differentiating has been normalized.

The objective counterpart does exactly that, deliberately: _raw_sensitivity_accessor (assembly.py:724, 748-754) does

norm = sim_data.normalization or {}
...
records = norm.get(col_name)
if not records:
    return tensor_sens(col_name, row)
...
normalized[col_name] = _normalized_sensitivity(records, col_name, sim_data, tensor_sens)

Meanwhile the penalty itself is read out of the in-place-normalized column: Constraint._difference_argmax -> Constraint.index(sim_data_dict, keys) (constraint.py:485-490) returns sim_data_dict[model][suffix][observable], the array Data.normalize already rescaled. So for normalization = <method> plus an active .con/.prop constraint on a measured observable under job_type lbfgs/trf/gntr, the penalty is P(normalized q) while assemble_constraint_gradient (line 1029) and assemble_constraint_hessian (line 1059) return P'(q) or P''(q) times d(RAW q)/dtheta.

The error is not a uniform scale that a line search could absorb. For peak the normalized column is raw/N with N = max over the same theta-dependent column, so d(raw/N)/dtheta = (1/N)*d(raw)/dtheta - (n_i/N)*d(raw_ref)/dtheta. The code keeps only N times the first term and drops the second entirely, so the returned direction is wrong: measured above, a pure initial-condition scale whose true normalized-column derivative is exactly 0 gets a spurious -1.0976 component while the k component comes back 100x too large.

Data.normalize_to_peak (data.py:517-542) rescales self.data and appends a NormalizationRecord but leaves Data.output_sensitivities in raw units, which is why the objective path has to fold the rule and why the constraint path silently gets it wrong. The Hessian shares the same accessor, so gntr's constraint curvature block inherits the same error squared (outer(grad_q, grad_q)).

Reproduction

Read-only script (no repo files touched), run as:
cd /Users/l119605/Code/PyBNF && uv run --extra tests --extra petab python /private/tmp/claude-503/-Users-l119605-Code-PyBNF/eb55c5a3-1006-4c5d-9d1e-9d14fe8045f0/scratchpad/fd.py

The script builds an exactly-known model Stot(t) = S0exp(-kt) on t = [0,1,2,3] with the matching OutputSensitivities (d_param = -tS0exp(-kt), d_ic = exp(-kt)), routes k->PARAM 'k' and S0->IC 'S()' with factor 1, then in two configurations:

A (control, no normalization): AtConstraint('Stot', '>', 90.0, 'm', 'tc', weight=2.0, atvar=None, atval=2.0)
B (the bug): sim.normalize('peak') first, then AtConstraint('Stot', '>', 0.9, ...) — the same constraint expressed in normalized units.

In each, it compares assemble_constraint_gradient([cset], sdd, routings, free) at k=0.3, S0=100 against a central difference (h = 1e-6, relative for S0) of ConstraintSet.total_penalty over the SAME construct-then-normalize-then-score path.

Observed:
normalize=False thresh=90.0 penalty 70.23767278
analytic [219.52465444, -1.09762327]
central [219.52465443, -1.09762327] ratio [1.0, 1.0]
normalize=True thresh=0.9 penalty 0.702376728
analytic [219.52465444, -1.09762327]
central [2.19524654, 0.0] ratio [100.0, inf]

Expected: analytic == central in both cases. The k column is off by exactly the peak normalizer N = 100 and the S0 column is -1.0976 where the true derivative is exactly 0.

A second script confirms the two accessors disagree on the same normalized Data:
_raw_sensitivity_accessor(...)('Stot', 2) -> [-1.1, 0.0] (objective side, normalization folded)
_constraint_sensitivity_accessor(...)('m','tc','Stot',2) -> [-110.0, 0.55] (constraint side, raw tensor)

End-to-end config that reaches this (not run; traced in config.py): edition 2, job_type = lbfgs (or trf/gntr), an experiment whose data: lists both E.exp and E.prop, normalization Stot = peak (or a whole-fit normalization = peak) where Stot is measured in E.exp and constrained in E.prop.

Verification notes

CONFIRMED by a central-difference test on the real code path (normalize in place -> total_penalty -> assemble_constraint_gradient), which is exactly the order a gradient fit runs (AlgorithmBase.score_result at pybnf/algorithms/base.py:857 calls res.normalize(...) and then scores res.simdata in place; GradientBase.gradient_at at pybnf/algorithms/optimizers/gradient_base.py:712-752 then hands that same res.simdata to assemble_constraint_gradient).

I tried hard to refute it and could not:

1. "The factor is applied one layer up" — it is not. Constraint.penalty_gradient / penalty_curvature (pybnf/constraint.py:978, 1035) only multiply raw_sens by the scalar penalty slope P'(q) / P''(q); they never consult sim_data.normalization. assemble_constraint_gradient (assembly.py:1029) and assemble_constraint_hessian (assembly.py:1059) apply only _sampling_scale_factors (the native->sampling d theta/d u), which is orthogonal. There is no other layer.

2. "The tensor is already normalized" — it is not. Data.normalize_to_peak (pybnf/data.py:517-542) rescales self.data in place and appends a NormalizationRecord; it never touches Data.output_sensitivities. So the A — Gradient plumbing: capability gate + preserve the output-sensitivity tensor through net execution #447 tensor stays in raw units while the scored column is rescaled. The objective path exists precisely because of this: _raw_sensitivity_accessor (assembly.py:724, 748-754) reads norm = sim_data.normalization and folds the chain via _normalized_sensitivity. _constraint_sensitivity_accessor (assembly.py:1069-1109) has no such lookup and no refusal.

3. "It is unreachable" — it is reachable on the supported new-era surface. Config._load_experiment_constraints (pybnf/config.py:2192-2216) builds ConstraintSet(base, data_key) where data_key is the EXPERIMENT's data key, and _partition_experiment_data (config.py:2060) explicitly allows one experiment to list both .exp and .con/.prop files. _resolve_normalization_grid (config.py:4200-4280) keys the resolved grid by that same data_key and by the experiment's measured column names, and Result.normalize (pybnf/algorithms/core.py:75-81) applies it to simdata[m][data_key]. So an experiment with both .exp data and a .prop/.con constraint on a measured observable normalizes exactly the column the constraint reads. Constraints can also reach a normalized column from any other suffix through the documented 'suffix.Observable' form (Constraint.get_key, constraint.py:445-483). No config validation gates constraints against normalization, and nothing gates gradient job types against constraints (gradient_base.py:751-755 supports them).

4. "A test pins it deliberately" — no. tests/test_gradient_assembly.py has ~15 constraint-gradient tests (lines 3324-4280, all on _constraint_sim) and none of them normalizes the sim Data; the file's 58 normalize mentions are all in the objective section.

Measured numbers (scratchpad script, read-only, uv run --extra tests --extra petab):

• control, no normalization: analytic [219.52465444, -1.09762327] vs central difference [219.52465443, -1.09762327] — agrees to 1e-9.
• peak-normalized column, same model, same constraint expressed in normalized units: analytic [219.52465444, -1.09762327] vs central difference [2.19524654, 0.0]. The k column is 100x too large (exactly the peak scale N = S0 = 100) and the S0 column is -1.0976 where the truth is exactly 0, because the quotient rule's second term (-n_i * s_ref / N) is dropped — a pure IC scale cancels under peak normalization, so its true normalized-column derivative is zero.
  I also confirmed the two accessors disagree on the same Data: _raw_sensitivity_accessor -> [-1.1, 0.0]; _constraint_sensitivity_accessor -> [-110.0, 0.55].

One correction to the claim as filed: its parenthetical "the returned gradient is bit-identical with and without normalize('peak')" is not reproducible as written — without normalization that constraint is satisfied and the gradient is [0, 0] (the threshold has to be re-expressed in the other units). The substance of the claim is unaffected; the correct statement is that the gradient ignores the normalization records entirely.

Where

pybnf/gradient/assembly.py:1082 — severity medium, confidence certain. Repro executed: True.

[wshlavacek]

Line-number drift after #727 (which changed data.py's normalizations and gradient/assembly.py's z-score rule). This issue is not fixed by it — I re-checked _constraint_sensitivity_accessor's raw_sens on f3bd79f0: it still goes straight to _sensitivity(...) with no sim_data.normalization read and no _normalized_sensitivity wrapper, so the diagnosis above is unchanged.

Current locations:

cited above	now
data.py:517-542 (normalize_to_peak)	data.py:596
assembly.py:1082-1108 (_constraint_sensitivity_accessor)	assembly.py:1079
assembly.py:1029 (assemble_constraint_gradient)	assembly.py:1013
assembly.py:1059 (assemble_constraint_hessian)	assembly.py:1048

One thing worth knowing for whoever picks this up: #727 made _zscore_sensitivity mask its reductions to a column's measured rows, so if the fix here routes the constraint path through _normalized_sensitivity (which is what the objective path does), it inherits that nan-awareness for free rather than needing its own.