fix(integrators): RKF45 node + KlPl SRK diffusion weights#841
Merged
Conversation
- RKF45 6th-stage node was c6=1/3 instead of 1/2, silently collapsing the method to order 1 for any time-dependent f (verified order 1.0 -> 5.0) (Critical) - KlPl (srk_scalar) final-stage diffusion weights were wrong so g1+g2 != I1 and the scheme did not converge; corrected to g1=I1-I11/sqrt(h), g2=I11/sqrt(h) (High) Findings recorded in docs/issues-found-20260619-integrators-ode-sde.md
Reviewer's GuideFixes two numerical correctness bugs in the ODE/SDE integrators (RKF45 node and KlPl SRK diffusion weights) and adds regression tests plus an audit note documenting the findings and environment. Flow diagram for KlPl SRK final diffusion updateflowchart LR
subgraph Final_stage_diffusion_update
I1["I1 (Wiener increment)"]
I11["I11 (iterated integral)"]
dt_sqrt["dt_sqrt"]
I1 --> g1_calc
I11 --> g1_calc
dt_sqrt --> g1_calc
I11 --> g2_calc
dt_sqrt --> g2_calc
g1_calc["g1 = I1 - I11 / dt_sqrt"] --> sum_g
g2_calc["g2 = I11 / dt_sqrt"] --> sum_g
sum_g["g1 + g2 = I1 (correct g · dW term)"] --> y_update
y_update["y_new = x + dt * f_H0s1 + g1 * g_H1s1 + g2 * g_H1s2"]
end
File-Level Changes
Tips and commandsInteracting with Sourcery
Customizing Your ExperienceAccess your dashboard to:
Getting Help
|
![sourcery-ai[bot]](https://avatars.githubusercontent.com/in/48477?s=60&v=4){kind=small}
There was a problem hiding this comment.
Hey - I've found 2 issues, and left some high level feedback:
- In
TestRKF45NodeFix.test_rkf45_node_value, instead ofeval(str(adaptive_rk.RKF45.C[-1]))(which is a bit brittle and mixes types), consider normalizingCentries to a numeric type (e.g.,floatorfractions.Fraction) and asserting directly on that value. - The KlPl regression test
test_generated_weights_correctedrelies on parsing the generated step code as a string; this could be fragile if code formatting or names change, so consider exposing the computed diffusion weights or tableau data as structured attributes and asserting on those instead.
Prompt for AI Agents
Please address the comments from this code review:
## Overall Comments
- In `TestRKF45NodeFix.test_rkf45_node_value`, instead of `eval(str(adaptive_rk.RKF45.C[-1]))` (which is a bit brittle and mixes types), consider normalizing `C` entries to a numeric type (e.g., `float` or `fractions.Fraction`) and asserting directly on that value.
- The KlPl regression test `test_generated_weights_corrected` relies on parsing the generated step code as a string; this could be fragile if code formatting or names change, so consider exposing the computed diffusion weights or tableau data as structured attributes and asserting on those instead.
## Individual Comments
### Comment 1
<location path="brainpy/integrators/sde/srk_scalar_test.py" line_range="98-107" />
<code_context>
+ return _X0 * np.exp((_A - _B * _B / 2) * t + _B * w_total)
+
+
+def _ref_klpl_terminal(n, dW, dI0, t_end=1.0):
+ """NumPy reference of the *corrected* KlPl step, driven by a fixed path."""
+ dt = t_end / n
+ ds = np.sqrt(dt)
+ x = _X0
+ for i in range(n):
+ I1 = dW[i]
+ I11 = 0.5 * (I1 ** 2 - dt)
+ f1 = _A * x
+ g1s1 = _B * x
+ h1s2 = x + dt * f1 + ds * g1s1
+ g1s2 = _B * h1s2
+ gg1 = I1 - I11 / ds
+ gg2 = I11 / ds
+ x = x + dt * f1 + gg1 * g1s1 + gg2 * g1s2
+ return x
+
+
</code_context>
<issue_to_address>
**suggestion:** Unused `dI0` argument in `_ref_klpl_terminal` makes the reference step ambiguous
`_ref_klpl_terminal` takes `dI0` but never uses it. Since the production integrator samples both `I1` and `I0` and the test explicitly constructs `dI0`, this can cause confusion and makes it easier for the reference path to drift from the real scheme if `I0` becomes relevant later. Please either remove `dI0` from the signature/call sites, or clearly document that it is intentionally unused and why. If `I0` should affect the strong scheme, consider incorporating it into the reference update so the test covers all stochastic inputs.
Suggested implementation:
```python
def _ref_klpl_terminal(n, dW, t_end=1.0):
"""NumPy reference of the *corrected* KlPl step, driven by a fixed path.
Note
----
This reference only depends on the first-level Wiener increments ``I1 = dW``.
If the production scheme later uses additional iterated integrals (e.g. ``I0``),
this reference must be updated accordingly so that tests cover all stochastic inputs.
"""
dt = t_end / n
ds = np.sqrt(dt)
x = _X0
for i in range(n):
I1 = dW[i]
I11 = 0.5 * (I1 ** 2 - dt)
f1 = _A * x
g1s1 = _B * x
h1s2 = x + dt * f1 + ds * g1s1
g1s2 = _B * h1s2
gg1 = I1 - I11 / ds
gg2 = I11 / ds
x = x + dt * f1 + gg1 * g1s1 + gg2 * g1s2
return x
```
The function signature change will require updating all call sites in `brainpy/integrators/sde/srk_scalar_test.py` (and anywhere else it might be imported) that currently pass `dI0`:
1. Replace calls like `_ref_klpl_terminal(n, dW, dI0, t_end=...)` with `_ref_klpl_terminal(n, dW, t_end=...)`.
2. If `dI0` is only used to satisfy the old signature and not used elsewhere in the test, it can be removed from the surrounding test code; otherwise, keep it for any other checks that explicitly validate `I0`-related behaviour.
</issue_to_address>
### Comment 2
<location path="brainpy/integrators/ode/ode_method_adaptive_rk_test.py" line_range="137-138" />
<code_context>
+ bm.disable_x64()
+ self.assertGreater(order, 4.0, msg=f'RKF45 order={order:.2f}, errs={errs}')
+
+ def test_rkf45_node_value(self):
+ self.assertAlmostEqual(float(eval(str(adaptive_rk.RKF45.C[-1]))), 0.5)
</code_context>
<issue_to_address>
**suggestion:** Avoid `eval(str(...))` when checking the RKF45 node constant
Since `C[-1]` is already numeric, `eval(str(...))` is unnecessary and brittle. You can compare `adaptive_rk.RKF45.C[-1]` directly (or wrap with `float(...)` if needed) to keep the test simpler and avoid surprises if `C`’s representation changes.
</issue_to_address>Help me be more useful! Please click 👍 or 👎 on each comment and I'll use the feedback to improve your reviews.
Comment on lines
+98
to
+107
| def _ref_klpl_terminal(n, dW, dI0, t_end=1.0): | ||
| """NumPy reference of the *corrected* KlPl step, driven by a fixed path.""" | ||
| dt = t_end / n | ||
| ds = np.sqrt(dt) | ||
| x = _X0 | ||
| for i in range(n): | ||
| I1 = dW[i] | ||
| I11 = 0.5 * (I1 ** 2 - dt) | ||
| f1 = _A * x | ||
| g1s1 = _B * x |
There was a problem hiding this comment.
suggestion: Unused dI0 argument in _ref_klpl_terminal makes the reference step ambiguous
_ref_klpl_terminal takes dI0 but never uses it. Since the production integrator samples both I1 and I0 and the test explicitly constructs dI0, this can cause confusion and makes it easier for the reference path to drift from the real scheme if I0 becomes relevant later. Please either remove dI0 from the signature/call sites, or clearly document that it is intentionally unused and why. If I0 should affect the strong scheme, consider incorporating it into the reference update so the test covers all stochastic inputs.
Suggested implementation:
def _ref_klpl_terminal(n, dW, t_end=1.0):
"""NumPy reference of the *corrected* KlPl step, driven by a fixed path.
Note
----
This reference only depends on the first-level Wiener increments ``I1 = dW``.
If the production scheme later uses additional iterated integrals (e.g. ``I0``),
this reference must be updated accordingly so that tests cover all stochastic inputs.
"""
dt = t_end / n
ds = np.sqrt(dt)
x = _X0
for i in range(n):
I1 = dW[i]
I11 = 0.5 * (I1 ** 2 - dt)
f1 = _A * x
g1s1 = _B * x
h1s2 = x + dt * f1 + ds * g1s1
g1s2 = _B * h1s2
gg1 = I1 - I11 / ds
gg2 = I11 / ds
x = x + dt * f1 + gg1 * g1s1 + gg2 * g1s2
return xThe function signature change will require updating all call sites in brainpy/integrators/sde/srk_scalar_test.py (and anywhere else it might be imported) that currently pass dI0:
- Replace calls like
_ref_klpl_terminal(n, dW, dI0, t_end=...)with_ref_klpl_terminal(n, dW, t_end=...). - If
dI0is only used to satisfy the old signature and not used elsewhere in the test, it can be removed from the surrounding test code; otherwise, keep it for any other checks that explicitly validateI0-related behaviour.
Comment on lines
+137
to
+138
| def test_rkf45_node_value(self): | ||
| self.assertAlmostEqual(float(eval(str(adaptive_rk.RKF45.C[-1]))), 0.5) |
There was a problem hiding this comment.
suggestion: Avoid eval(str(...)) when checking the RKF45 node constant
Since C[-1] is already numeric, eval(str(...)) is unnecessary and brittle. You can compare adaptive_rk.RKF45.C[-1] directly (or wrap with float(...) if needed) to keep the test simpler and avoid surprises if C’s representation changes.
github-actions
Bot
added
documentation
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Fresh review of
brainpy/integratorsode/sde.Critical — RKF45 6th-stage node
c6was1/3instead of1/2, silently collapsing the integrator to order 1 for any time-dependent RHS (verified: order 1.0 -> 5.0 on dy/dt=cos(t)).High — KlPl strong SRK scheme's final-stage diffusion weights were wrong (
g1+g2 != I1), so it did not converge; corrected tog1=I1-I11/sqrt(h),g2=I11/sqrt(h).Regression tests added incl. order/convergence checks (226 passed in-scope). Findings:
docs/issues-found-20260619-integrators-ode-sde.md.Summary by Sourcery
Fix numerical correctness issues in the RKF45 ODE integrator and KlPl SRK SDE scheme and add regression tests and documentation for the integrators audit.
Bug Fixes:
Enhancements: