The weight-penalty gradient of a split-at constraint (`A at t1 < B at t2`) reads B from t2 through t1 instead of at t2, so lbfgs converges to a wrong optimum, or crashes when t2 is after t1

Start in pybnf/constraint.py at lines 858-859, then read _static_penalty_gradient and _difference_argmax around lines 830-886 alongside the split interval logic at line 1200. Compare how the probit and logit branches pass imax2. Reproduce with the unit-level check in the issue; done means split-at weight gradients match finite differences for t2 before, after, and equal to t1, with no empty-slice error.

The weight-penalty gradient of a split-at constraint (A at t1 < B at t2) reads B from t2 through t1 instead of at t2, so lbfgs converges to a wrong optimum, or crashes when t2 is after t1

For a static (weight) penalty, get_penalty_gradient calls _static_penalty_gradient(sim_data_dict, raw_sens, n_param, imin, imax, once, imin2) (pybnf/constraint.py:858-859). It passes imin2 but not imax2, and _static_penalty_gradient (:879) has no imax2 parameter at all. Its call to _difference_argmax (:886) therefore falls back to imax2 = imax (:830-831). SplitAtConstraint._penalty_intervals (:1200) returns imin=fi1, imax=fi1+1, imin2=fi2, imax2=fi2+1, so the second operand is sliced as col2[fi2:fi1+1] instead of the single row fi2. The penalty itself goes through get_difference, which uses the correct imax2 (:635-656). The penalty and its gradient therefore disagree:

• t2 before t1 (fi2 < fi1): B is sliced over several rows. The argmax and the satisfied/violated test run over rows the penalty never reads. The gradient is then read at row fi1 + a of A, which is outside the interval, and at row fi2 + a of B. A satisfied constraint (penalty 0) gets a nonzero gradient.
• t2 after t1 (fi2 > fi1): the slice is empty and np.argmax raises ValueError.

The probit (tolerance) and logit (scale) branches pass imax2 (:861-867) and agree with finite differences.

Failure scenario

A is produced at a constant rate kA, and B decays from 20 at rate 0.5, so A(t) = kA·t and B(t) = 20·e^(-t/2). The data are A = t, so the true fit is kA = 1 with objective 0. The constraint tc.Obs_A at time=6 < tc.Obs_B at time=2 weight 1 holds there (A(6) = 6 < B(2) = 7.358), so it is inactive and should not move the fit.

In the gradient, B is sliced over rows t=2..6. Its minimum is at t=6, where 6 − 0.996 > 0, so the constraint is scored as violated. The gradient is then read at A(t=10), which adds a spurious +10 to d/dkA. L-BFGS-B gets the true objective values but this wrong gradient. It reports "Stop criterion satisfied" at the point where the spurious term cancels the data gradient: kA = 1 − 10/385 = 0.974026, objective 0.12987. With weight 100, kA is driven to its lower bound.

When the times are swapped (tc.Obs_A at time=2 < tc.Obs_B at time=6), the fit dies instead. That constraint is violated at kA = 1 (A(2) = 2 > B(6) = 0.996). Its correct fit balances the data against the penalty: kA = 1 − 2/385 = 0.994805, objective 0.999064. The same constraint written with B(6) as a constant fits to exactly that point.

Reproduction

ab.bngl:

begin model
begin parameters
  kA 1
  kB 0.5
  B0 20
end parameters
begin molecule types
  A()
  B()
end molecule types
begin seed species
  A() 0
  B() B0
end seed species
begin observables
  Molecules Obs_A A()
  Molecules Obs_B B()
end observables
begin reaction rules
  0 -> A() kA
  B() -> 0 kB
end reaction rules
end model

a.exp: header # time Obs_A, then rows t t for t = 0, 1, ..., 10.

split.prop: one line, tc.Obs_A at time=6 < tc.Obs_B at time=2 weight 1

fit.conf:

output_dir = out
edition = 2
model: ab.bngl
bngl_backend = bngsim
job_type = lbfgs
objective = sos
experiment: tc, data: a.exp
experiment: qual, data: split.prop, t_end: 10
uniform_var = kA 0.2 3.0
population_size = 2
max_iterations = 50
parallel_count = 1
random_seed = 7

Run pybnf -c fit.conf. Only the .prop line changes between the five runs below. The two control rows write the same constraint with B's value as a constant, so they show the correct fit.

.prop line	best kA	objective	correct
tc.Obs_A < 7.357588823428847 at time=6 weight 100 (control: B(2) as a constant)	1.0000000000000007	2.96e-31	kA = 1, objective 0
tc.Obs_A at time=6 < tc.Obs_B at time=2 weight 1	0.9740259740259751	0.12987	kA = 1, objective 0
tc.Obs_A at time=6 < tc.Obs_B at time=2 weight 100	0.2 (lower bound)	123.2	kA = 1, objective 0
tc.Obs_A < 0.9957413673572789 at time=2 weight 1 (control: B(6) as a constant)	0.9948051948051941	0.999064	kA = 1 − 2/385 = 0.994805, objective 0.999064
tc.Obs_A at time=2 < tc.Obs_B at time=6 weight 1	Sorry, an unknown error occurred: ValueError: attempt to get argmax of an empty sequence		kA = 0.994805, objective 0.999064

Both lbfgs starts give the same result in each run. The traceback of the last run goes through gradient_at -> assemble_constraint_gradient -> pybnf/constraint.py:886.

Unit-level check, with no simulator involved. Each "parameter" is one cell of the simulated output, so the gradient shows which rows are read:

import os, tempfile
import numpy as np
from pybnf.constraint import ConstraintSet
from pybnf.data import Data

t = np.arange(11.0)
B = np.full(11, 10.0); B[4] = 1.0                     # A(t) = t; B = 10 except B(4) = 1
sim = {'m': {'s': Data.from_columns(np.column_stack([t, t, B]), ['time', 'A', 'B'])}}

def raw_sens(model, suffix, obs, row):                # one "parameter" per cell: A[0..10], B[0..10]
    v = np.zeros(22); v[(0 if obs == 'A' else 11) + row] = 1.0
    return v

for line in ['A at 6 < B at 2 weight 1', 'A at 2 < B at 6 weight 1', 'A at 6 < B at 6 weight 1']:
    prop = os.path.join(tempfile.mkdtemp(), 's.prop')
    open(prop, 'w').write(line + '\n')
    cs = ConstraintSet('m', 's'); cs.load_constraint_file(prop)
    c = cs.constraints[0]
    try:
        g = c.penalty_gradient(sim, raw_sens, {}, 22)
        g = {f"{'AB'[j // 11]}(t={j % 11})": float(g[j]) for j in np.flatnonzero(g)}
    except ValueError as e:
        g = f'ValueError: {e}'
    print(f'{line}: penalty = {c.penalty(sim)}, gradient = {g}')

Output:

A at 6 < B at 2 weight 1: penalty = 0.0, gradient = {'A(t=8)': 1.0, 'B(t=4)': -1.0}
A at 2 < B at 6 weight 1: penalty = 0.0, gradient = ValueError: attempt to get argmax of an empty sequence
A at 6 < B at 6 weight 1: penalty = 0.0, gradient = {}

In both split cases the penalty is 0 (A(6) = 6 < B(2) = 10, and A(2) = 2 < B(6) = 10). Central finite differences of the penalty are therefore zero in every cell. The correct gradient is {} in all three lines. The keyword-free form (A at 6 < B at 2, which defaults to weight 1) gives the same wrong result. The tolerance and scale forms of the first line match finite differences.

Reachability

• Job type: job_type = lbfgs, or any fit on the scalar-gradient path. gradient_at adds assemble_constraint_gradient for every constraint set, with no filter by constraint type (pybnf/algorithms/optimizers/gradient_base.py:750-754).
• Constraint: any split-at constraint with the default or weight penalty. The split syntax is in the grammar (pybnf/constraint.py:261) and documented (docs/config.rst:224, A at 5 < B at C=6). qualitative_loss = hinge coerces tolerance/scale lines to the static model, so it reaches the same code.
• Observables: they must be suffix-qualified to an experiment with measured data (e.g. tc.Obs_A). A constraint on a data-less qualitative suffix is refused on the gradient path before this code runs.
• Consequences: when t2 is before t1, the fit converges without complaint to a wrong optimum. When t2 is after t1, the ValueError is not a GradientNotSupported, so it escapes gradient_at's handler and ends the fit with the generic unknown-error message.

The fix is to pass imax2 through _static_penalty_gradient to _difference_argmax, as the probit and logit branches already do.

Where

• pybnf/constraint.py:858-859 (static dispatch drops imax2)
• pybnf/constraint.py:879 (_static_penalty_gradient signature has no imax2)
• pybnf/constraint.py:886 (_difference_argmax called without imax2)
• pybnf/constraint.py:830-831 (fallback imax2 = imax)
• pybnf/constraint.py:1200 (SplitAtConstraint._penalty_intervals)

Related: #887, #890.

Found in a whole-codebase audit for silently wrong results (2026-09-23); the reproduction above was re-run independently of the original finding.