Fix ReturnData::{preeq_wrms,posteq_wrms} with FSA and check_sensi_steadystate_conv_=True #2071

dweindl · 2023-05-08T17:02:46Z

Documentation says, that ReturnData::{preeq_wrms,posteq_wrms} refers to the weighted root-mean-square of the rhs when steadystate was reached, independent of Solver::check_sensi_steadystate_conv_, but this wasn't what happened. Now it is.

Closes #2066

…adystate_conv_=True Documentation says, that ReturnData::{preeq_wrms,posteq_wrms} refers to the `weighted root-mean-square of the rhs when steadystate was reached`, independent of `Solver::check_sensi_steadystate_conv_`, but this wasn't what happened. Now it is.

codecov · 2023-05-08T17:21:12Z

Codecov Report

Merging #2071 (e55b7e6) into develop (ba70f54) will not change coverage.
The diff coverage is n/a.

Additional details and impacted files

@@           Coverage Diff            @@
##           develop    #2071   +/-   ##
========================================
  Coverage    54.90%   54.90%           
========================================
  Files           30       30           
  Lines         4617     4617           
========================================
  Hits          2535     2535           
  Misses        2082     2082

Flag	Coverage Δ
petab	`54.90% <ø> (ø)`
sbmlsuite	`?`

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

FFroehlich

👍

dweindl · 2023-05-09T07:51:28Z

Any idea why this change would affect state sensitivities? 🤔

Details

```python __________________________________________________________________________________________ test_manual_preequilibration ___________________________________________________________________________________________

preeq_fixture = (<amici.amici.ModelPtr; proxy of <Swig Object of type 'std::unique_ptr< amici::Model > *' at 0x7f58c37fb750> >, <Swig ... (2/2 measurements & 2/2 sigmas set)
10x0 event-resolved datapoints
(0/0 measurements & 0/0 sigmas set)

, ...)

def test_manual_preequilibration(preeq_fixture):
    """Manual preequilibration"""

    model, solver, edata, edata_preeq, \
        edata_presim, edata_sim, pscales, plists = preeq_fixture

    settings = itertools.product(pscales, plists)

    for pscale, plist in settings:

        model.setInitialStates([])
        model.setInitialStateSensitivities([])
        model.setParameterList(plist)
        model.setParameterScale(pscale)

        # combined
        rdata_auto = amici.runAmiciSimulation(model, solver, edata)
        assert rdata_auto.status == amici.AMICI_SUCCESS

        # manual preequilibration
        rdata_preeq = amici.runAmiciSimulation(model, solver, edata_preeq)
        assert rdata_preeq.status == amici.AMICI_SUCCESS

        # manual reinitialization + presimulation
        x0 = rdata_preeq['x'][0, :]
        x0[1] = edata_presim.fixedParameters[0]
        x0[2] = edata_presim.fixedParameters[1]
        sx0 = rdata_preeq['sx'][0, :, :]
        sx0[:, 1] = 0
        sx0[:, 2] = 0
        model.setInitialStates(x0)
        model.setInitialStateSensitivities(sx0.flatten())
        rdata_presim = amici.runAmiciSimulation(model, solver, edata_presim)
        assert rdata_presim.status == amici.AMICI_SUCCESS

        # manual reinitialization + simulation
        x0 = rdata_presim['x'][0, :]
        x0[1] = edata_sim.fixedParameters[0]
        x0[2] = edata_sim.fixedParameters[1]
        sx0 = rdata_presim['sx'][0, :, :]
        sx0[:, 1] = 0
        sx0[:, 2] = 0
        model.setInitialStates(x0)
        model.setInitialStateSensitivities(sx0.flatten())
        rdata_sim = amici.runAmiciSimulation(model, solver, edata_sim)
        assert rdata_sim.status == amici.AMICI_SUCCESS

        for variable in ['x', 'sx']:

          assert_allclose(

                rdata_auto[variable],
                rdata_sim[variable],
                atol=1e-6, rtol=1e-6,
                err_msg=str(dict(pscale=pscale, plist=plist, variable=variable))
            )

python/tests/test_preequilibration.py:118:

args = (<function assert_allclose..compare at 0x7f58fd350b80>, array([[[-1.23878552, 0.9791515 , -0.65434205, -1.235...3906, 0.20831407],
[ 0.03262003, -0.01297652, 1.10028436, 0.00754647,
-1.34715986, 1.30699336]]]))
kwds = {'equal_nan': True, 'err_msg': "{'pscale': <ParameterScaling.log10: 2>, 'plist': [3, 1, 2, 4], 'variable': 'sx'}", 'header': 'Not equal to tolerance rtol=1e-06, atol=1e-06', 'verbose': True}

@wraps(func)
def inner(*args, **kwds):
    with self._recreate_cm():

      return func(*args, **kwds)

E AssertionError:
E Not equal to tolerance rtol=1e-06, atol=1e-06
E {'pscale': <ParameterScaling.log10: 2>, 'plist': [3, 1, 2, 4], 'variable': 'sx'}
E Mismatched elements: 1 / 48 (2.08%)
E Max absolute difference: 2.52982694e-06
E Max relative difference: 4.61690539e-06
E x: array([[[-1.238786, 0.979151, -0.654342, -1.235868, 2.097963,
E 0.37669 ],
E [-3.794804, -1.890767, 0.526177, 4.689174, -0.811134,...
E y: array([[[-1.238786, 0.979152, -0.654342, -1.235868, 2.097963,
E 0.37669 ],
E [-3.794805, -1.890765, 0.526177, 4.689176, -0.811135,...

/usr/lib/python3.11/contextlib.py:81: AssertionError

</p>
</details>

FFroehlich · 2023-05-09T08:41:53Z

Any idea why this change would affect state sensitivities? 🤔

__________________________________________________________________________________________ test_manual_preequilibration ___________________________________________________________________________________________

preeq_fixture = (<amici.amici.ModelPtr; proxy of <Swig Object of type 'std::unique_ptr< amici::Model > *' at 0x7f58c37fb750> >, <Swig ...   (2/2 measurements & 2/2 sigmas set)
  10x0 event-resolved datapoints
    (0/0 measurements & 0/0 sigmas set)
>, ...)

    def test_manual_preequilibration(preeq_fixture):
        """Manual preequilibration"""
    
        model, solver, edata, edata_preeq, \
            edata_presim, edata_sim, pscales, plists = preeq_fixture
    
        settings = itertools.product(pscales, plists)
    
        for pscale, plist in settings:
    
            model.setInitialStates([])
            model.setInitialStateSensitivities([])
            model.setParameterList(plist)
            model.setParameterScale(pscale)
    
            # combined
            rdata_auto = amici.runAmiciSimulation(model, solver, edata)
            assert rdata_auto.status == amici.AMICI_SUCCESS
    
            # manual preequilibration
            rdata_preeq = amici.runAmiciSimulation(model, solver, edata_preeq)
            assert rdata_preeq.status == amici.AMICI_SUCCESS
    
            # manual reinitialization + presimulation
            x0 = rdata_preeq['x'][0, :]
            x0[1] = edata_presim.fixedParameters[0]
            x0[2] = edata_presim.fixedParameters[1]
            sx0 = rdata_preeq['sx'][0, :, :]
            sx0[:, 1] = 0
            sx0[:, 2] = 0
            model.setInitialStates(x0)
            model.setInitialStateSensitivities(sx0.flatten())
            rdata_presim = amici.runAmiciSimulation(model, solver, edata_presim)
            assert rdata_presim.status == amici.AMICI_SUCCESS
    
            # manual reinitialization + simulation
            x0 = rdata_presim['x'][0, :]
            x0[1] = edata_sim.fixedParameters[0]
            x0[2] = edata_sim.fixedParameters[1]
            sx0 = rdata_presim['sx'][0, :, :]
            sx0[:, 1] = 0
            sx0[:, 2] = 0
            model.setInitialStates(x0)
            model.setInitialStateSensitivities(sx0.flatten())
            rdata_sim = amici.runAmiciSimulation(model, solver, edata_sim)
            assert rdata_sim.status == amici.AMICI_SUCCESS
    
            for variable in ['x', 'sx']:
>               assert_allclose(
                    rdata_auto[variable],
                    rdata_sim[variable],
                    atol=1e-6, rtol=1e-6,
                    err_msg=str(dict(pscale=pscale, plist=plist, variable=variable))
                )

python/tests/test_preequilibration.py:118: 
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _

args = (<function assert_allclose.<locals>.compare at 0x7f58fd350b80>, array([[[-1.23878552,  0.9791515 , -0.65434205, -1.235...3906,  0.20831407],
        [ 0.03262003, -0.01297652,  1.10028436,  0.00754647,
         -1.34715986,  1.30699336]]]))
kwds = {'equal_nan': True, 'err_msg': "{'pscale': <ParameterScaling.log10: 2>, 'plist': [3, 1, 2, 4], 'variable': 'sx'}", 'header': 'Not equal to tolerance rtol=1e-06, atol=1e-06', 'verbose': True}

    @wraps(func)
    def inner(*args, **kwds):
        with self._recreate_cm():
>           return func(*args, **kwds)
E           AssertionError: 
E           Not equal to tolerance rtol=1e-06, atol=1e-06
E           {'pscale': <ParameterScaling.log10: 2>, 'plist': [3, 1, 2, 4], 'variable': 'sx'}
E           Mismatched elements: 1 / 48 (2.08%)
E           Max absolute difference: 2.52982694e-06
E           Max relative difference: 4.61690539e-06
E            x: array([[[-1.238786,  0.979151, -0.654342, -1.235868,  2.097963,
E                     0.37669 ],
E                   [-3.794804, -1.890767,  0.526177,  4.689174, -0.811134,...
E            y: array([[[-1.238786,  0.979152, -0.654342, -1.235868,  2.097963,
E                     0.37669 ],
E                   [-3.794805, -1.890765,  0.526177,  4.689176, -0.811135,...

/usr/lib/python3.11/contextlib.py:81: AssertionError

Looking at the code, I am wondering how this was even working before. while(true), but no break?

FFroehlich · 2023-05-09T08:58:31Z

l 668/9

    /* If run after Newton's method checks again if it converged */
    wrms_ = getWrms(model, sensitivityFlag);

also appears to be defunct. And with the current implementation, we still need some way of exiting the loop in case of adjoint sensitivies or no FSA convergence check.

And it might be worthwhile investigating how this even worked in the first place? Were we always just using up all integration steps and then converging via the convergence check for newton updates?

dweindl · 2023-05-09T09:45:27Z

l 668/9
    /* If run after Newton's method checks again if it converged */
    wrms_ = getWrms(model, sensitivityFlag);
also appears to be defunct. And with the current implementation, we still need some way of exiting the loop in case of adjoint sensitivies or no FSA convergence check.

And it might be worthwhile investigating how this even worked in the first place? Were we always just using up all integration steps and then converging via the convergence check for newton updates?

Not sure I can follow. We still have

AMICI/src/steadystateproblem.cpp

Line 704 in 2ea8281

wrms_ = getWrms(model, sensitivityFlag);

AMICI/src/steadystateproblem.cpp

Lines 715 to 716 in 2ea8281

    
           if (wrms_ < conv_thresh) 
        
               break; // converged

to not exhaust the number of steps, not?

dweindl · 2023-05-09T10:11:20Z

Okay, the failing tests were just because the getWrmsFSA check would be ignored in subsequent iterations. That should be fixed now.

Not sure what to do about

    /* If run after Newton's method checks again if it converged */
    wrms_ = getWrms(model, sensitivityFlag);

Simply remove? Should we even end up here if Newton's method already converged?

Move the convergence check to the beginning of the loop, so the wrms_ computed there is actually checked?

FFroehlich · 2023-05-09T10:41:46Z

Not sure what to do about
    /* If run after Newton's method checks again if it converged */
    wrms_ = getWrms(model, sensitivityFlag);
Simply remove? Should we even end up here if Newton's method already converged?

No we shouldn't, but we might start at steady state and not have run Newton's method.

Move the convergence check to the beginning of the loop, so the wrms_ computed there is actually checked?

Yes.

FFroehlich · 2023-05-09T10:46:26Z

Okay, the failing tests were just because the getWrmsFSA check would be ignored in subsequent iterations. That should be fixed now.

Not sure I follow, but I do think this fixes the issue that I mentioned above about not having a break statement in case there is no getWrmsFSA check.

FFroehlich

👍

dweindl requested a review from a team as a code owner May 8, 2023 17:02

FFroehlich approved these changes May 8, 2023

View reviewed changes

Merge branch 'develop' into fix_2066

77fb5fd

dweindl self-assigned this May 9, 2023

Prettify test output

7dc70b9

Merge branch 'develop' into fix_2066

c0a9119

dweindl and others added 2 commits May 9, 2023 11:06

rtol

f79b937

Merge branch 'develop' into fix_2066

2ea8281

don't skip getWrmsFSA check in subsequent iterations...

0ee586c

dweindl added 3 commits May 9, 2023 12:49

check before stepping

35339b6

remove redundant wrms

1beca97

style

e55b7e6

dweindl added this to the v0.17.0 milestone May 9, 2023

dweindl requested a review from FFroehlich May 9, 2023 12:13

FFroehlich approved these changes May 9, 2023

View reviewed changes

dweindl merged commit 252da79 into develop May 9, 2023

dweindl deleted the fix_2066 branch May 9, 2023 13:09

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Fix ReturnData::{preeq_wrms,posteq_wrms} with FSA and check_sensi_steadystate_conv_=True #2071

Fix ReturnData::{preeq_wrms,posteq_wrms} with FSA and check_sensi_steadystate_conv_=True #2071

dweindl commented May 8, 2023 •

edited

Loading

codecov bot commented May 8, 2023 •

edited

Loading

FFroehlich left a comment

dweindl commented May 9, 2023 •

edited

Loading

FFroehlich commented May 9, 2023

FFroehlich commented May 9, 2023

dweindl commented May 9, 2023

dweindl commented May 9, 2023 •

edited

Loading

FFroehlich commented May 9, 2023 •

edited

Loading

FFroehlich commented May 9, 2023

FFroehlich left a comment

Fix ReturnData::{preeq_wrms,posteq_wrms} with FSA and check_sensi_steadystate_conv_=True #2071

Fix ReturnData::{preeq_wrms,posteq_wrms} with FSA and check_sensi_steadystate_conv_=True #2071

Conversation

dweindl commented May 8, 2023 • edited Loading

codecov bot commented May 8, 2023 • edited Loading

Codecov Report

FFroehlich left a comment

Choose a reason for hiding this comment

dweindl commented May 9, 2023 • edited Loading

FFroehlich commented May 9, 2023

FFroehlich commented May 9, 2023

dweindl commented May 9, 2023

dweindl commented May 9, 2023 • edited Loading

FFroehlich commented May 9, 2023 • edited Loading

FFroehlich commented May 9, 2023

FFroehlich left a comment

Choose a reason for hiding this comment

dweindl commented May 8, 2023 •

edited

Loading

codecov bot commented May 8, 2023 •

edited

Loading

dweindl commented May 9, 2023 •

edited

Loading

dweindl commented May 9, 2023 •

edited

Loading

FFroehlich commented May 9, 2023 •

edited

Loading