Fix piecewise/Heaviside handling #2234

dweindl · 2023-12-12T20:46:14Z

Fixes a bug where Heaviside functions weren't correctly substituted, potentially resulting in DiracDelta's in the model equations. The problem was, the substitution targets were collected from an expanded expression, but the actual substitution was attempted on the non-expanded expression, which did not always succeed.

Fixed by not expanding the expressions beforehand. Closes NaNs in model expressions with DiracDelta due to unfortunate term ordering #2231
Fixes redundant trigger-function processing

* Substitute non-time-dependent Heavisides * Substitute inside the expanded expression, otherwise not all substitution targets are found Closes AMICI-dev#2231

codecov · 2023-12-12T21:07:16Z

Codecov Report

Merging #2234 (4e727f6) into develop (eba9557) will increase coverage by 27.98%.
The diff coverage is 100.00%.

Additional details and impacted files

@@             Coverage Diff              @@
##           develop    #2234       +/-   ##
============================================
+ Coverage    48.65%   76.64%   +27.98%     
============================================
  Files           91       91               
  Lines        15109    15112        +3     
============================================
+ Hits          7352    11583     +4231     
+ Misses        7757     3529     -4228

Flag	Coverage Δ
cpp	`73.11% <ø> (?)`
cpp_python	`37.08% <ø> (ø)`
petab	`53.77% <100.00%> (+0.02%)`	⬆️
python	`77.45% <100.00%> (+15.13%)`	⬆️
sbmlsuite	`∅ <ø> (∅)`

Flags with carried forward coverage won't be shown. Click here to find out more.

Files	Coverage Δ
python/sdist/amici/de_export.py	`91.87% <100.00%> (+3.66%)`	⬆️
python/sdist/amici/import_utils.py	`88.23% <100.00%> (+15.34%)`	⬆️
python/sdist/amici/sbml_import.py	`79.03% <ø> (+13.37%)`	⬆️

... and 56 files with indirect coverage changes

dweindl · 2023-12-12T21:18:00Z

Okay, this might have been a more widespread problem - this also fixes frequent Brannmark_JBC2010 issues (i.e. breaks my example at https://amici.readthedocs.io/en/latest/example_errors.html#Steady-state-computation-failed :-|)

dweindl · 2023-12-12T21:50:56Z

Worth checking whether this change also affects the issues encountered in #1539

dweindl · 2023-12-12T22:10:18Z

Okay, this might have been a more widespread problem - this also fixes frequent Brannmark_JBC2010 issues (i.e. breaks my example at https://amici.readthedocs.io/en/latest/example_errors.html#Steady-state-computation-failed :-|)

No. Brannmark_JBC2010 didn't have any stray Heavisides. Maybe just calling sympy.expand() did the trick there...

dweindl · 2023-12-13T08:02:16Z

=========================== short test summary info ============================
FAILED test_petab_benchmark.py::test_benchmark_gradient[Giordano_Nature2020-True] - AssertionError:                                                          direction  ...                       analysis_results
  directi...   [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...  ...  Empty DataFrame
  Columns: []
  Index: []
  
  [50 rows x 6 columns]
assert False
============= 1 failed, 36 passed, 7 skipped in 941.58s (0:15:41) ==============

Unable to reproduce locally 🙄

FFroehlich · 2023-12-13T14:46:40Z

python/sdist/amici/de_export.py

@@ -2767,9 +2767,11 @@ def _process_heavisides(
        if heavisides:
            # only apply subs if necessary
            for heaviside_sympy, heaviside_amici in heavisides:
-                dxdt = dxdt.subs(heaviside_sympy, heaviside_amici)
+                dt_expanded = dt_expanded.subs(


was the issue that expand led to manipulation of the trigger function which was then missed during this substitution?

Not super happy that we manipulate the rhs by applying expand here as I think we discussed at some point we want to keep simplification optional. I don't see why the expand is necessary here in the first place. Looks like this was introduced in #1352

was the issue that expand led to manipulation of the trigger function which was then missed during this substitution?

Correct.

Not super happy that we manipulate the rhs by applying expand here as I think we discussed at some point we want to keep simplification optional. I don't see why the expand is necessary here in the first place. Looks like this was introduced in #1352

# expanding the rhs will in general help to collect the same # heaviside function

So I guess, not expanding the expression may lead to more root functions to be tracked than strictly necessary in certain cases. I don't know how common this issue really is, and how much of a performance impact each trigger function has.

I am open to dropping expand.

Looks like this was introduced in #1352

Specifically in f27004f
Doesn't really provide more context, though.

was the issue that expand led to manipulation of the trigger function which was then missed during this substitution?

Correct.

Not super happy that we manipulate the rhs by applying expand here as I think we discussed at some point we want to keep simplification optional. I don't see why the expand is necessary here in the first place. Looks like this was introduced in #1352

# expanding the rhs will in general help to collect the same # heaviside function

So I guess, not expanding the expression may lead to more root functions to be tracked than strictly necessary in certain cases. I don't know how common this issue really is, and how much of a performance impact each trigger function has.

I am open to dropping expand.

But we are anyways calling sp.simplify in _get_unique_root. If expand really makes a difference, we could do sp.simplify(root_found.expand() - root.get_val()), but I would be surprised if that's really the case. We can map multiple Heaviside functions to the same root finding function, so I would prefer keeping all symbolic manipulations in _get_unique_root.

If expand really makes a difference

I don't really have any evidence.

I find it hard to imagine though that we'd get a relevant performance hit from tracking e.g. 3 equivalent root functions instead of just one.

I'm fine with either. Just dropping the expand requires fewest changes.

Well, I don't think the expand ever really improved anything in terms of reducing the number of root functions. Either it would simply have had no effect at all (i.e. the same roots would have been extracted with or without expand), or it would have prevented replacing the Heavisides because the to-be-substituted subexpressions weren't found in the original expression. So I'd say it's safe to drop it completely.

This reverts commit fd3f076.

dweindl added 4 commits December 12, 2023 21:45

Fix piecewise/Heavisides handling

e0fbeec

* Substitute non-time-dependent Heavisides * Substitute inside the expanded expression, otherwise not all substitution targets are found Closes AMICI-dev#2231

perf: Handle unique roots only

fd227ec

..

fcc36cc

..

bbfe5f7

dweindl requested a review from FFroehlich December 12, 2023 21:20

dweindl linked an issue Dec 12, 2023 that may be closed by this pull request

NaNs in model expressions with DiracDelta due to unfortunate term ordering #2231

Closed

Remove broken error example

fd3f076

dweindl mentioned this pull request Dec 12, 2023

Use Model.setAlwaysCheckFinite(True) in tests #2235

Closed

deterministic unique

31518e2

Merge branch 'develop' into fix_2231

a0b061b

dweindl changed the title ~~Fix piecewise/Heavisides handling~~ Fix piecewise/Heaviside handling Dec 13, 2023

dweindl added 3 commits December 13, 2023 09:25

seed++?

1f7dc9a

seed++

f76678e

skip

c854930

dweindl marked this pull request as ready for review December 13, 2023 10:30

dweindl requested a review from a team as a code owner December 13, 2023 10:30

Merge branch 'develop' into fix_2231

a8d382c

FFroehlich reviewed Dec 13, 2023

View reviewed changes

dweindl added 3 commits December 13, 2023 15:59

don't expand

2771ba0

try unskip Giordano_Nature2020

bcbea7e

Revert "Remove broken error example"

4e727f6

This reverts commit fd3f076.

dweindl requested a review from FFroehlich December 13, 2023 16:13

FFroehlich approved these changes Dec 14, 2023

View reviewed changes

dweindl merged commit 594b07e into AMICI-dev:develop Dec 14, 2023
23 checks passed

dweindl deleted the fix_2231 branch December 14, 2023 09:45

dweindl mentioned this pull request Dec 15, 2023

Non-deterministic code generation #2241

Closed

dweindl mentioned this pull request Oct 17, 2024

SBML import: Incorrect processing of piecewise/Heaviside function #2545

Closed

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Fix piecewise/Heaviside handling #2234

Fix piecewise/Heaviside handling #2234

dweindl commented Dec 12, 2023 •

edited

Loading

codecov bot commented Dec 12, 2023 •

edited

Loading

dweindl commented Dec 12, 2023

dweindl commented Dec 12, 2023

dweindl commented Dec 12, 2023 •

edited

Loading

dweindl commented Dec 13, 2023

FFroehlich Dec 13, 2023

dweindl Dec 13, 2023

dweindl Dec 13, 2023

FFroehlich Dec 13, 2023

dweindl Dec 13, 2023

dweindl Dec 13, 2023 •

edited

Loading

Fix piecewise/Heaviside handling #2234

Fix piecewise/Heaviside handling #2234

Conversation

dweindl commented Dec 12, 2023 • edited Loading

codecov bot commented Dec 12, 2023 • edited Loading

Codecov Report

dweindl commented Dec 12, 2023

dweindl commented Dec 12, 2023

dweindl commented Dec 12, 2023 • edited Loading

dweindl commented Dec 13, 2023

FFroehlich Dec 13, 2023

Choose a reason for hiding this comment

dweindl Dec 13, 2023

Choose a reason for hiding this comment

dweindl Dec 13, 2023

Choose a reason for hiding this comment

FFroehlich Dec 13, 2023

Choose a reason for hiding this comment

dweindl Dec 13, 2023

Choose a reason for hiding this comment

dweindl Dec 13, 2023 • edited Loading

Choose a reason for hiding this comment

dweindl commented Dec 12, 2023 •

edited

Loading

codecov bot commented Dec 12, 2023 •

edited

Loading

dweindl commented Dec 12, 2023 •

edited

Loading

dweindl Dec 13, 2023 •

edited

Loading