Fix piecewise/Heaviside handling

[dweindl]

• 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

[codecov]

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]

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]

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

[dweindl]

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]

=========================== 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]

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

[dweindl]

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.

[dweindl]

Looks like this was introduced in #1352

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

[FFroehlich]

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.

[dweindl]

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.

[dweindl]

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.