Repository navigation
fix: time-integral and joined sensitivities, sparse output variables, and first_state - #5794
Conversation
3a85437 to
c4f4039
Compare
|
| Project | PyBaMM |
| Branch | a3b/processed-variable-followups |
| Testbed | bare-metal |
🐰 View full continuous benchmarking report in Bencher
⚠️ WARNING: Truncated view!The full continuous benchmarking report exceeds the maximum length allowed on this platform.
Codecov Report✅ All modified and coverable lines are covered by tests. Additional details and impacted files@@ Coverage Diff @@
## a3b/processed-variable-fixes #5794 +/- ##
================================================================
- Coverage 98.10% 98.10% -0.01%
================================================================
Files 346 346
Lines 35458 35482 +24
================================================================
+ Hits 34786 34809 +23
- Misses 672 673 +1 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
c4f4039 to
cee7afb
Compare
martinjrobins
left a comment
There was a problem hiding this comment.
Thanks @BradyPlanden. I must admit I didn't quite understand the sparsity fix, could you explain it a bit more?
rtimms
left a comment
There was a problem hiding this comment.
The fixes in this PR look right. Leaving the initial condition out of time-integral sensitivities matches finite differences for I and I**2 on both the full-state and output_variables paths. The shared-boundary sensitivity join works where #5785 raised a CasADi error. The sparse dense fill is correct, and _y0 is carried through +, from_sub_solutions, copy and reduce_solution.
Two sensitivity bugs sit underneath this PR (both reproduced on #5785's code); details are inline. Missing tests:
first_state.sensitivitieson anoutput_variablessolve. This would catch the first inline issue.- An
output_variablestime integral whose post-sum expression depends on an input. - A sparse output variable combined with sensitivities.
81052db to
eb2ff88
Compare
eb2ff88 to
c629964
Compare
c629964 to
1ec388e
Compare
1ec388e to
ae38764
Compare
martinjrobins
left a comment
There was a problem hiding this comment.
one more query @BradyPlanden:
rtimms
left a comment
There was a problem hiding this comment.
No blocking issues. I checked:
- the time-integral chain rule, which is evaluated at the integral and leaves out the initial condition
jacp_stackedagainst the sensitivity inputs only, in the C++ parameter orderrows_repeatedfor time-major data- the
first_stateslicing fix - the sparse scatter:
find()is column-major and every output is a column - the new C++
y_init/yS_init/yS_termlifetimes across repeated and threaded solves
One minor point is inline. A follow-up, not for this PR: on the DAE path that stops on an event at t0 (IDAKLUSolverOpenMP.inl ~L375), IDACalcIC doesn't write back to yy/yyS unless IDAGetConsistentIC/IDAGetSensConsistentIC is called, so y_init, yS_init and yS_term hold the pre-initialisation guesses there.
| new_sol.closest_event_idx = other.closest_event_idx | ||
| new_sol._y0 = self._y0 | ||
| new_sol._y0_sensitivities = self._y0_sensitivities | ||
| new_sol._y_event_sensitivities = other._y_event_sensitivities |
There was a problem hiding this comment.
[P3] The one-point shortcut above (L1303-1314) doesn't do this. It copies other._y_event but keeps self._y_event_sensitivities. If an outputs-only step stops at once on an event at the shared time, last_state pairs the new y_event with the previous step's state sensitivities, which then seed the next restart. Suggest also setting new_sol._y_event_sensitivities = other._y_event_sensitivities in that branch.
There was a problem hiding this comment.
This is essentially the same as Martin's thread above, so I've kept self's here. other can be that carried-over last_state, so copying its event sensitivities would break the restart (d1eacc2 tests this). In this case, an event at t0, yS_term holds the seeded values because of the IDACalcIC gap you noted, and those are self's anyway.
9e52a62 to
d1eacc2
Compare
martinjrobins
left a comment
There was a problem hiding this comment.
LGTM!, thanks @BradyPlanden
postfix_sensitivities summed the sensitivities with postfix_sum, which adds the time integral's initial condition. The initial condition does not depend on the inputs, so every sensitivity of ExplicitTimeIntegral(c, 5) came back 5 too large, on both the full-state and output_variables paths. The sum over time is now kept apart from the initial condition, which only the value adds. A new test checks the sensitivity against its analytic value on both paths, and the output_variables join test now uses a nonzero initial condition, which the joined value counts once and its sensitivity not at all.
IDAKLUSolver returns only the structural nonzeros of an output variable, such as the zero-flux boundary edges of "Electrolyte current density [A.m-2]", while as_computed() data holds every entry. ProcessedVariableComputed inferred one layout for a joined variable from its first CasADi function and filled every segment from the first segment's data, so joining two such solutions, or the steps of an experiment, raised an IndexError, and joining one to a full-state solution raised a ValueError. The solver now scatters the nonzeros into every entry, and _unroll_nnz is gone. The output_variables matrix gains a 1D and a 2D sparse case, and its join test also joins two output_variables solutions. Its comparisons take an atol, since the full-state path rounds values within about 1e-16 of zero to zero.
…gments Adding solutions solved with sensitivities dropped the repeated time point at a shared boundary from the states but not from their sensitivities, so reading the sensitivities of any variable raised a CasADi RuntimeError. _merge_sensitivities now drops that time point's rows, in __add__ and from_sub_solutions alike. A time integral's value sums over every segment at once, but its sensitivities were summed per segment and stacked, one row per segment. They now go through a single postfix over all segments.
postfix_sensitivities evaluated the derivative of a post-sum expression, such as ExplicitTimeIntegral(c, 0) ** 2, at the expression's own value rather than at the integral, so the square gave 2 I^2 dI/da instead of 2 I dI/da. It now takes the summed variable at each time point, as postfix does, and forms the integral itself; both callers pass those raw values. test_explicit_time_integral and test_discrete_data_sum encoded the wrong derivative and are corrected. At the discrete test's a = 2 the model reproduces its data, and the square root of a zero sum has no derivative, so that one sensitivity assertion is skipped. The joined time-integral sensitivity test gains a post-sum case.
first_state sliced the last n_states rows of the first segment's sensitivities, which belong to its last time point, so it carried the sensitivities of the wrong state. It now takes the first rows, as the states themselves do.
first_state of a solve that returns only output variables read the model's y0full when first accessed. That vector is the solver's input: it carries the initial sensitivities after the states, holds the algebraic states from before consistent initialization (a DFN's potentials were off by up to 0.08 V), and is overwritten by the model's next solve. In an experiment whose cycle repeats a step, the cycle's first state and every "Change in ..." summary variable came from the last repeat, halving "Change in throughput capacity [A.h]". pybammsolvers now returns y_init alongside y_term when it saves only output variables: the states at t0 after consistent initialization, from which the outputs at t0 are computed. IDAKLUSolver keeps it on the Solution, and __add__, from_sub_solutions, copy and reduce_solution carry it. It ships in the unreleased pybammsolvers 0.11.0, which pybamm already requires. Solution gains a _y0 slot, which a solution pickled before it existed loads as None.
The fill of a sparse output variable indexed the dense entries by the sparsity's row indices, which are the flat indices only because every output variable is a column vector. It now uses find(), the flat indices for any shape, and reads the sparsity from sparsity_out() instead of evaluating the variable at zero.
An output variable's parameter Jacobian, which pybammsolvers uses to compute its sensitivities during an output_variables solve, was taken with respect to every input, but the solver reads one column per sensitivity input. With calculate_sensitivities=["b"] and inputs a and b, the solver took the derivative with respect to a, returned it as "a", left "b" empty, and wrote past the end of its buffer for the column of b. The Jacobian is now taken with respect to the sensitivity inputs, in their order, and the output variables' sensitivities are named by them. The post-sum expression of a time integral was differentiated with respect to every input too, so its (1 x inputs) derivative broadcast against the (1 x sensitivities) chain-rule term: "b" came back as the derivative with respect to a. postfix_sensitivities now takes the sensitivity names, keeps the other inputs as values, and differentiates with respect to those names only.
An output_variables solve filled the states' sensitivities from sol.yS, which then holds the outputs' sensitivities as (n_t, n_outputs, n_p): indexing it by parameter picked a time point. With more sensitivity inputs than time points the solve raised an IndexError; otherwise first_state and last_state held one entry per output, of the wrong quantity, and solver.step started the next step from them, raising "y0 has wrong number of cols" or, when the counts matched, giving wrong sensitivities (du/da = 0.23 for an exact 0.86). pybammsolvers now returns yS_init and yS_term alongside y_init and y_term: the states' sensitivities at the first and last time points, one row per sensitivity parameter. IDAKLUSolver keeps them on the Solution for first_state and last_state, and __add__, from_sub_solutions, copy and reduce_solution carry them as they carry y_init and y_event. The states' sensitivities over time are empty, as the states are. Solution gains _y0_sensitivities and _y_event_sensitivities slots, which a solution pickled before they existed loads as None. BaseSolver._check_restart_sensitivities rejected every sensitivity restart from an output_variables solve, whose step boundary had no state sensitivities to seed from. It now rejects only a solve without them, and the step test that expected the rejection checks the restarted sensitivities instead.
Each changelog entry is now one line, leaving the detail to the PR, and eleven two-line comments become one line.
changelog entry
from_sub_solutions takes _y_event_sensitivities from the last kept segment, not sols[-1], matching __add__'s single-sample short-circuit, which keeps self's. The comment above it covered only closest_event_idx, so it now says why the event sensitivities follow it too. Simulation joins a trailing single-sample duplicate when every step of a cycle is infeasible: it carries last_state over, which holds no event sensitivities of its own. Taking them from sols[-1] would drop them, so last_state.sensitivities would come from the outputs' rows and a restart would raise. A new test joins a solve with its last_state through both __add__ and from_sub_solutions and checks them.
d1eacc2 to
fba9d10
Compare
Description
Fixes eight older bugs in computed variables' sensitivities and a solution's first and last states, one commit each:
"Electrolyte current density [A.m-2]", raised. IDAKLU returns only their structural nonzeros, whichIDAKLUSolvernow scatters into place.first_statetook its sensitivities from the first segment's last time point.first_stateof anoutput_variablessolve read the model's latesty0, halving"Change in throughput capacity [A.h]"in experiments that repeat a step.output_variablessolves stored the outputs' sensitivities as the states', sostepraised or continued from wrong sensitivities.Two tests that asserted the wrong derivative for an expression of a time integral are corrected.
pybammsolversnow returnsy_init,yS_initandyS_termfrom outputs-only solves. They ship in the unreleasedpybammsolvers0.11.0, whichpybammalready requires, so the two packages need releasing together.Type of change
Add an entry under
# [Unreleased]in CHANGELOG.md, in one of## Breaking changes,## Deprecated,## Features, or## Bug fixes(include PR number; breaking and deprecation entries need a one-line migration note). Internal-only PRs — refactor, docs, CI, tests — can skip this. See RELEASE.md for the full policy.Important checks:
Please confirm the following before marking the PR as ready for review:
nox -s pre-commitnox -s testsnox -s doctests