sys-bio / roadrunner

libRoadRunner: A high-performance SBML simulator
http://libroadrunner.org/
Other
39 stars 24 forks source link

Problem with RK45 integrator? #877

Open CiaranWelsh opened 2 years ago

CiaranWelsh commented 2 years ago
from roadrunner import RoadRunner

sbml_file = "AModel.sbml"
r = RoadRunner(sbml_file)

integrators = ('cvode', 'gillespie', 'rk4', 'rk45', 'euler')
data = {}
for integrator in integrators:
    r.setIntegrator(integrator)
    r.resetAll()
    arr = r.simulate(0, 10, 11)
    data[integrator] = arr

import matplotlib.pyplot as plt
import seaborn as sns

sns.set_context('talk')

fig = plt.figure()
for integrator, da in data.items():
    plt.plot(da['time'], da['[S1]'], label=integrator)

plt.legend(loc='best')
sns.despine(fig=fig, top=True, right=True)
plt.show()

with AModel.sbml

produces :

image

luciansmith commented 2 years ago

I can confirm that while the rk45 results in the rrtests still do what they used to, but it fails 100% of the SBML Test Suite tests, many of which are ridiculously simple.

jpsluka commented 1 year ago

There appears to be an inconsistency in what RK45 actually is versus what it is described as. In the documents at https://libroadrunner.readthedocs.io/en/latest/PythonAPIReference/cls_Integrator.html the sections on Gillespie and RK45 have the same text. But RK45 is a continuous method and not stochastic like the documentation says:

Gillespie RoadRunner’s implementation of the standard Gillespie Direct Method SSA. The granularity of this simulator is individual molecules and kinetic processes are stochastic. Results will, in general, be different in each run, but a sufficiently large ensemble of runs should be statistically correct.

RK45 RoadRunner’s implementation of the standard Gillespie Direct Method SSA. The granularity of this simulator is individual molecules and kinetic processes are stochastic. Results will, in general, be different in each run, but a sufficiently large ensemble of runs should be statistically correct.

So the two are described with the same text.

https://www.hindawi.com/journals/complexity/2017/1232868/ says:

2.4. Runge-Kutta-Fehlberg Algorithm (RK45) "The forward Euler algorithm introduced in the previous section is a first-order numerical method provided for didactic purposes. Several numerical methods have been introduced to increase the accuracy of deterministic simulations and to decrease simulation runtime. A popular algorithm is the Runge-Kutta-Fehlberg algorithm (RK45). This algorithm represents the standard choice of several numerical integrators when the model does not exhibit stiffness. It is also often used in hybrid simulation strategies for the simulation of the fast part of the network. The reader can find a comprehensive guide to numerical methods of ODEs in [7, 13]."

So the roadrunner documentation seems wrong, and/or the RK45 integrator implemented might be something other than RK45.

hsauro commented 1 year ago

Something to look at.

Herbert

On Tue, May 9, 2023 at 11:16 AM James (Jim) Sluka @.***> wrote:

There appears to be an inconsistency in what RK45 actually is versus what it is described as. In the documents at https://libroadrunner.readthedocs.io/en/latest/PythonAPIReference/cls_Integrator.html https://urldefense.com/v3/__https://libroadrunner.readthedocs.io/en/latest/PythonAPIReference/cls_Integrator.html__;!!K-Hz7m0Vt54!i3XqvYJU4MILtWdq7Z3O2gAjlTS5beIMjNtLOJp8HB3U8rH-D2KJ8c1t5lVSbObgREuX9zTCkONDU0-zJevrPmxs1mpR0Q$ the sections on Gillespie and RK45 have the same text. But RK45 is a continuous method and not stochastic like the documentation says:

Gillespie RoadRunner’s implementation of the standard Gillespie Direct Method SSA. The granularity of this simulator is individual molecules and kinetic processes are stochastic. Results will, in general, be different in each run, but a sufficiently large ensemble of runs should be statistically correct.

RK45 RoadRunner’s implementation of the standard Gillespie Direct Method SSA. The granularity of this simulator is individual molecules and kinetic processes are stochastic. Results will, in general, be different in each run, but a sufficiently large ensemble of runs should be statistically correct.

So the two are described with the same text.

https://www.hindawi.com/journals/complexity/2017/1232868/ https://urldefense.com/v3/__https://www.hindawi.com/journals/complexity/2017/1232868/__;!!K-Hz7m0Vt54!i3XqvYJU4MILtWdq7Z3O2gAjlTS5beIMjNtLOJp8HB3U8rH-D2KJ8c1t5lVSbObgREuX9zTCkONDU0-zJevrPmwNzSHP6Q$ says:

2.4. Runge-Kutta-Fehlberg Algorithm (RK45) "The forward Euler algorithm introduced in the previous section is a first-order numerical method provided for didactic purposes. Several numerical methods have been introduced to increase the accuracy of deterministic simulations and to decrease simulation runtime. A popular algorithm is the Runge-Kutta-Fehlberg algorithm (RK45). This algorithm represents the standard choice of several numerical integrators when the model does not exhibit stiffness. It is also often used in hybrid simulation strategies for the simulation of the fast part of the network. The reader can find a comprehensive guide to numerical methods of ODEs in [7, 13]."

So the roadrunner documentation seems wrong, and/or the RK45 integrator implemented might be something other than RK45.

— Reply to this email directly, view it on GitHub https://urldefense.com/v3/__https://github.com/sys-bio/roadrunner/issues/877*issuecomment-1540651864__;Iw!!K-Hz7m0Vt54!i3XqvYJU4MILtWdq7Z3O2gAjlTS5beIMjNtLOJp8HB3U8rH-D2KJ8c1t5lVSbObgREuX9zTCkONDU0-zJevrPmwQgplGGQ$, or unsubscribe https://urldefense.com/v3/__https://github.com/notifications/unsubscribe-auth/AAIBSDQR4SOZK5X7X2WLWGTXFKCXPANCNFSM5FTUTX7A__;!!K-Hz7m0Vt54!i3XqvYJU4MILtWdq7Z3O2gAjlTS5beIMjNtLOJp8HB3U8rH-D2KJ8c1t5lVSbObgREuX9zTCkONDU0-zJevrPmxYGw9dCw$ . You are receiving this because you are subscribed to this thread.Message ID: @.***>

-- Herbert Sauro, Professor University of Washington, Bioengineering 206-685-2119, www.sys-bio.org Mobile: 206-880-8093 @.*** Books: http://books.analogmachine.org/

luciansmith commented 1 year ago

Thanks for the note! Fixed with https://github.com/sys-bio/roadrunner/pull/1108