AMICI-dev / AMICI

Advanced Multilanguage Interface to CVODES and IDAS
https://amici.readthedocs.io/
Other
108 stars 30 forks source link

Incorrect handling of event assignments to compartments #2426

Closed dweindl closed 5 months ago

dweindl commented 5 months ago

I think there is a problem with event assignments to compartments (targeted by rate rules, others aren't supported yet), if 1) there are multiple event assignments for a single event, and 2a) the assignment to the compartment is not the last assignment, or 2b) there are assignments to multiple compartments

In this case, here: https://github.com/AMICI-dev/AMICI/blob/da19a550688e8232b23ace6d06722c66deb2703f/python/sdist/amici/sbml_import.py#L1675-L1677 formula has the value of the most recently processed event assignments, independent on whether that was the respective compartment change.

Reproducer for case 2a:

    import amici
    from amici.antimony_import import antimony2amici
    from amici.testing import TemporaryDirectoryWinSafe as TemporaryDirectory
    from numpy.testing import assert_allclose

    ant_model = """
    model test_events_multiple_assignments
        compartment event_target = 1
        event_target' = 0
        species species_in_event_target in event_target = 1
        unrelated = 1

        # those two events produce identical results, but they don't
        #  this one fails
        at (time > 10): unrelated = 2, event_target = 10

        #  this one works
        #at (time > 10): event_target = 10, unrelated = 2

        # (note that antimony adds the event assignments to the sbml model 
        # from right to left)
    end
    """
    module_name = "test_multiple_event_assignment_with_compartment"
    with TemporaryDirectory(prefix=module_name, delete=False) as outdir:
        antimony2amici(
            ant_model,
            model_name=module_name,
            output_dir=outdir,
            verbose=True,
        )
        model_module = amici.import_model_module(
            module_name=module_name, module_path=outdir
        )
        amici_model = model_module.getModel()
        assert amici_model.ne == 1
        assert amici_model.ne_solver == 0
        assert amici_model.nx_rdata == 3
        amici_model.setTimepoints(np.linspace(0, 15, 200))
        amici_solver = amici_model.getSolver()
        rdata = amici.runAmiciSimulation(amici_model, amici_solver)
        assert rdata.status == amici.AMICI_SUCCESS
        idx_event_target = amici_model.getStateIds().index("event_target")
        idx_unrelated = amici_model.getStateIds().index("unrelated")
        idx_species_in_event_target = amici_model.getStateIds().index(
            "species_in_event_target"
        )
        assert (rdata.x[rdata.ts < 10, idx_event_target] == 1).all()
        assert (rdata.x[rdata.ts > 10, idx_event_target] == 10).all()
        assert (rdata.x[rdata.ts < 10, idx_unrelated] == 1).all()
        assert (rdata.x[rdata.ts > 10, idx_unrelated] == 2).all()
        assert (rdata.x[rdata.ts < 10, idx_species_in_event_target] == 1).all()
        # fails, it's 0.5, because of division by the new value of `unrelated`
        assert_allclose(rdata.x[rdata.ts > 10, idx_species_in_event_target], 0.1, rtol=0, atol=1e-15)