[Montana]

Closes #1286.

What was wrong

A rate law shaped as a smooth narrow bump of a clock is stepped over entirely, with no error and no warning. Examples are a Gaussian exp(-((T - c)/s)^2) and a sech² 1/cosh((T - c)/s)^2. Away from the bump the rate underflows to exactly 0 (exp(-10^4) is 0 in doubles). So the right-hand side is exactly what it would be with no bump, CVODE's error estimate is 0, and its step grows until one step spans the bump. Corners and conditions get stops or roots (#714, #1068, #1127), but a smooth bump has neither.

The issue's model: X → Y at a Gaussian of area 0.5 in a counter T fed at rate 2, centred at T = 50 (t = 25), so Y(30) = 10(1 − e^(−1/4)) = 2.211992. On main c0877fa, Y(30) = 0 at w = 0.01 and w = 1, at rtol 1e-6 and 1e-12, at 301 output points, and in a sensitivity run. The same holds for sech² and for time() in place of the counter.

The fix

New python/bngsim/_smooth_pulses.py. smooth_pulse_stops(core, t_start, t_end) finds the bumps and returns a stop at each one's centre. A step that would have spanned the bump ends on its peak instead, where the rate is largest, so the error test sees it.

• Bumps are read from each functional rate law with functions and derived parameters inlined, parsed by _exprtk_to_sympy, with the run's parameter values put in:
  • exp(q), where q is a quadratic in one clock with a negative leading coefficient (a Gaussian at any scale or offset); its centre is −b/2a;
  • cosh(u), sech(u) or tanh(u), where u is affine in one clock. That covers sech² spelled 1/cosh(u)^2 or 1 - tanh(u)^2 (tanh saturates to exactly 1, so the second is also 0 away from the bump). A plain tanh step gets a stop it does not need, which costs one restart.
• Clocks are literal time(), or any observable all of whose species change at one constant rate. That is found by two probes of the right-hand side at shifted states and times, as unit_rate_clock_candidates does, but at any rate: the issue's counter is fed at k_clock = 2, which the unit-rate machinery does not take for a clock. A counter's centre is converted to a time from its value at t_start and its rate.
• It is wired in Simulator._apply_crossing_stops only. The stops join the condition, periodic and table-function stops in merge_crossing_stops. They are not added to time_discontinuity_conditions(), so the sensitivity switch-time machinery, the state-switch roots and the other consumers of that list never see a bump as a switch.
• Only the value of a rate law, not its conditions. A bump inside a condition (if(exp(...) > 0.5, ...), a threshold written with tanh) is a crossing, which the condition scan places or declines (A narrow if() window is still stepped over when written as a number, as one non-monotone comparison, or reached through steady_state() (after #897) #904, A floor() or max() pulse on a counter that is not unit-rate gets no stop and no root, so it is stepped over: Y = 0 for 3.9347, with no error (the if() spelling is rooted) #1068). A stop from here would override what that scan decided. Two existing tests caught that in a first version: test_clock_condition_scan.py::test_a_text_the_evaluator_cannot_be_asked_about_is_declined and the moving-threshold window test in test_net_state_switch_stepping.py. Both pass unchanged now.
• Cheap where it does nothing. A model none of whose functions spells exp, cosh, sech or tanh is answered from the raw texts, without building the inlined context (0.01 ms measured). The answer is cached on everything it reads: function texts, parameter values, state and window. A repeat costs 0.03 ms.

One stop, at the peak. I first placed stops at the centre and four widths either side. Against a max_step reference I measured a Gaussian of widths 10, 3, 0.5 and 0.005, two bumps in one run, and a run that starts on a bump's rise, at rtol 1e-6 and 1e-12. The centre stop alone was within the tolerance in every case, as the three were, and is two restarts fewer per bump.

Tests (test_smooth_pulse_stops.py, 20)

test	without the stops (as on main)
test_a_bump_in_a_counter_is_integrated [Gaussian, sech² (cosh), sech² (1-tanh²)] × [w 0.01, 1] × [rtol 1e-6, 1e-12]: Y(30) = 2.211992	12 fail (Y = 0)
test_a_bump_in_literal_time_is_integrated	fail
test_dense_output_does_not_step_over_it_either (301 points)	fail
test_a_sensitivity_run_integrates_the_bump: Y and dY/dkx = 2.5·e^(−1/4) to 1e-5	fail
test_the_stop_is_at_the_centre (t = 25 on the rate-2 counter); test_a_stop_outside_the_run_is_not_placed	pass
controls, no stop placed: a Lorentzian (never underflows; caught today), an exponential decay of the clock, a bump of the state	pass

A test this changes

test_net_state_switch_stepping.py::test_a_released_hold_logs_what_main_logs failed with the bump stops in. That test checks that a state-switch hold gives CVODE's error log back once it ends, by counting the batch-end messages logged after the release. The count fell to 0 with a bump stop, and also on main with an ordinary condition stop.

I first took this for a bug and filed it as #1389. Traced, it is not one:

• the hold unmutes at the same point with a crossing stop as without one;
• a run with a crossing stop takes the stretch it has not looked inside one step at a time (An ODE run with a stop on every wrap of a remainder of a far clock stops adding to X after about 72 stops, with no error (100*mod(1e12 + time(), 0.125): 0.72 of the integral over 100 periods) #998), and a one-step call never spends a 20-step batch, so there was no batch end to log.

The test now runs with the bump stops patched out (pulse_stops=False), the run it was written for, and its comment says why. #1389 can be closed.

Mutations

mutant	fails
P2: no counter clocks	every counter test
P3: tanh not recognized	the four 1-tanh² cases
P4: cosh not recognized	the four cosh sech² cases
P5: conditions scanned for bumps too (_value_calls walks every node)	the two condition tests above

Measured

• The corpus: every .net model in the repo and the 1,323 BioModels SBML files, 1,468 of which load, over t = 0–100. No model gets a bump stop, so no corpus answer moves. The change is reached only by a model that writes such a bump.
• Full suite: 16,436 passed, 137 skipped, 1 xfailed (a test since removed, see above). pre-commit run --all-files (new files staged) and sphinx -W with docs/requirements-ci.txt are clean.

Not in this PR

• A bump of a clock read through anything but these spellings. A logistic's derivative exp(-u)/(1+exp(-u))^2, for example, has an exp of an affine argument, which is not taken for a bump. max_step under the bump's width remains the workaround there.
• The sech² sensitivity failure the issue notes at the end (the analytic derivative of 1/cosh(u)^2 is ∞/∞ far from the centre). That one fails loudly. The issue comment's 1 - tanh(u)^2 spelling avoids it, and this PR integrates that spelling.