sys-bio / roadrunner

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

Number of timepoints != steps+1 in fixed step size simulation #261

Closed matthiaskoenig closed 8 years ago

matthiaskoenig commented 8 years ago

Simulating for 10 steps must create 11 rows (time points) in the matrix (time=0.0 and 10 steps). Number of timepoints is now == steps. This was different before and must be fixed.

import tellurium as te
r = te.loada ('''
model feedback()
   // Reactions:
   J0: $X0 -> S1; (VM1 * (X0 - S1/Keq1))/(1 + X0 + S1 +   S4^h);
   J1: S1 -> S2; (10 * S1 - 2 * S2) / (1 + S1 + S2);
   J2: S2 -> S3; (10 * S2 - 2 * S3) / (1 + S2 + S3);
   J3: S3 -> S4; (10 * S3 - 2 * S4) / (1 + S3 + S4);
   J4: S4 -> $X1; (V4 * S4) / (KS4 + S4);

  // Species initializations:
  S1 = 0; S2 = 0; S3 = 0;
  S4 = 0; X0 = 10; X1 = 0;

  // Variable initialization:
  VM1 = 10; Keq1 = 10; h = 10; V4 = 2.5; KS4 = 0.5;
end''')

res = r.simulate(0, 40, 10)
print(res.shape)
(10, 5)

Please add tests for these things, i.e. all the following things have to be added to tests to avoid regressions. There should be much more test cases for the various calls to simulate.

r = te.loada ('''
model feedback()
   // Reactions:
   J0: $X0 -> S1; (VM1 * (X0 - S1/Keq1))/(1 + X0 + S1 +   S4^h);
   J1: S1 -> S2; (10 * S1 - 2 * S2) / (1 + S1 + S2);
   J2: S2 -> S3; (10 * S2 - 2 * S3) / (1 + S2 + S3);
   J3: S3 -> S4; (10 * S3 - 2 * S4) / (1 + S3 + S4);
   J4: S4 -> $X1; (V4 * S4) / (KS4 + S4);

  // Species initializations:
  S1 = 0; S2 = 0; S3 = 0;
  S4 = 0; X0 = 10; X1 = 0;

  // Variable initialization:
  VM1 = 10; Keq1 = 10; h = 10; V4 = 2.5; KS4 = 0.5;
end''')

# tests with the model
res = r.simulate(0, 40, 10)
self.assertEqual(11, res.shape[0])
self.assertEqual(res['time'][0], 0.0)
self.assertEqual(res['time'][-1], 40.0)
self.assertFalse(r.getIntegrator().getSetting('variable_step_size')

# use default steps of 50
res = r.simulate(0, 20)
self.assertEqual(51, res.shape[0])
self.assertEqual(res['time'][0], 0.0)
self.assertEqual(res['time'][-1], 20.0)
self.assertFalse(r.getIntegrator().getSetting('variable_step_size')
self.assertEqual(r.getIntegrator().getSetting('steps', 50)

r.getIntegrator().setSetting('variable_step_size', True)
res = r.simulate(0, 30)
self.assertTrue(res.shape[0]>=2)
self.assertEqual(res['time'][0], 0.0)
self.assertEqual(res['time'][-1], 30.0)
self.assertTrue(r.getIntegrator().getSetting('variable_step_size')
...
0u812 commented 8 years ago

This behavior was changed per @hsauro in #119: Herbert specified that steps should be the number of points, not the number of intervals.

hsauro commented 8 years ago

The reason that the third argument is the number of points is that its been that way in the original roadrunner since it was created (I think there was a short time when intervals was used but this was a mistake). The reason for choosing steps rather than intervals is that many novice and casual modelers especially biologists will find the concept of intervals more difficult to comprehend and will not understand the logic of specifying intervals when points appears to be more straightforward. Most researchers are used to thinking about plotting points on a graph or rows in a array of numbers, not intervals. For example in an excel spreadsheet we don't specify the number of rows in terms of intervals. Hence it is clearer for a user that the value represents how many points they will actually get on a graph or rows in an array.

Its possible we could get an alternative call that uses intervals instead of points or use an named parameter in simulate that specifies intervals, eg

m = r.simulate (0, 100, steps = 9) # yields 10 points

Th later is probably the best solution. Would that work for you Matthias?

Another possibility is to also add another named parameter that specifies the actual interval eg:

m = r.simulate (0, 100, intervalStep = 0.1)

I just remembered we have steps as a simulate option.

Herbert

On Mon, Jan 18, 2016 at 9:07 AM, Kyle Medley notifications@github.com wrote:

This behavior was changed per @hsauro https://github.com/hsauro in #119 https://github.com/sys-bio/roadrunner/issues/119: Herbert specified that steps should be the number of points, not the number of intervals.

— Reply to this email directly or view it on GitHub https://github.com/sys-bio/roadrunner/issues/261#issuecomment-172593041.

matthiaskoenig commented 8 years ago

steps +1 == number_of_points Please do not use this as the same. This will create a multitude of problems and missunderstandings.

If you have an argument to simulate and a integrator setting which is called "steps" it has to mean steps, for instance in r.getIntegrator.setSetting(steps, 50)

If you want to provide functionality for setting number of points, you should support an additional argument/setting number_of_points, which than set steps = number_of_points-1 (with check that >2).

Especially because other solvers understand steps also as (steps +1 == number_of_points), for instance Copasi. They call it intervals, but it is steps.

Is still broken with the integrator setting steps currently in the right direction ;). If you call

s = r.simulate(0,10)
s.shape

you will get 51 points! from steps=50

On Mon, Jan 18, 2016 at 6:07 PM, Kyle Medley notifications@github.com wrote:

This behavior was changed per @hsauro https://github.com/hsauro in #119 https://github.com/sys-bio/roadrunner/issues/119: Herbert specified that steps should be the number of points, not the number of intervals.

— Reply to this email directly or view it on GitHub https://github.com/sys-bio/roadrunner/issues/261#issuecomment-172593041.

Matthias König Email: konigmatt@googlemail.com Tel: + 49 30 20938450 Tel: + 49 176 81168480 https://www.livermetabolism.com/contact.html Humboldt-University Berlin, Institute for Theoretical Biology, Junior Group König Charité Berlin, Institute of Biochemistry, Computational Systems Biochemistry http://www.charite.de/sysbio/people/koenig/

matthiaskoenig commented 8 years ago

Yes. The steps argument is good + clear documentation that 4th positoal argument is number of points. Works for me. On Jan 18, 2016 7:35 PM, "Herbert Sauro" notifications@github.com wrote:

The reason that the third argument is the number of points is that its been that way in the original roadrunner since it was created (I think there was a short time when intervals was used but this was a mistake). The reason for choosing steps rather than intervals is that many novice and casual modelers especially biologists will find the concept of intervals more difficult to comprehend and will not understand the logic of specifying intervals when points appears to be more straightforward. Most researchers are used to thinking about plotting points on a graph or rows in a array of numbers, not intervals. For example in an excel spreadsheet we don't specify the number of rows in terms of intervals. Hence it is clearer for a user that the value represents how many points they will actually get on a graph or rows in an array.

Its possible we could get an alternative call that uses intervals instead of points or use an named parameter in simulate that specifies intervals, eg

m = r.simulate (0, 100, steps = 9) # yields 10 points

Th later is probably the best solution. Would that work for you Matthias?

Another possibility is to also add another named parameter that specifies the actual interval eg:

m = r.simulate (0, 100, intervalStep = 0.1)

I just remembered we have steps as a simulate option.

Herbert

On Mon, Jan 18, 2016 at 9:07 AM, Kyle Medley notifications@github.com wrote:

This behavior was changed per @hsauro https://github.com/hsauro in

119

https://github.com/sys-bio/roadrunner/issues/119: Herbert specified that steps should be the number of points, not the number of intervals.

— Reply to this email directly or view it on GitHub <https://github.com/sys-bio/roadrunner/issues/261#issuecomment-172593041 .

— Reply to this email directly or view it on GitHub https://github.com/sys-bio/roadrunner/issues/261#issuecomment-172618427.

hsauro commented 8 years ago

Not sure if I understand. I need to clarify what the various terms means, we have:

steps intervals number of points

any others?

I assume steps and intervals mean the same thing?

What would be useful is to list the kinds of calls to simulate that you'd like to see. For biologists I still want them them to be able to do following:

m = r.simulate (0, 10, 100) where 100 means the number rows in the returned array or points on the graph. For variable steps size the interval won't necessarily be the same at all points but the number of points returned should be 100.

Other than that I have no fixed ideas and one can have what ever parameters etc you need in the simulate call. We need to set up a comprehensive test suite for simulate given the number of arguments it has. Andy overloaded simulate a lot which doesn't help with the testing. I've created an issue for Kiri to test simulate() and make sure the documentation is consistent. In the mean time you could provide some specifications on how you want simulate to behave.

Herbert

On Mon, Jan 18, 2016 at 10:35 AM, Matthias König notifications@github.com wrote:

steps +1 == number_of_points Please do not use this as the same. This will create a multitude of problems and missunderstandings.

If you have an argument to simulate and a integrator setting which is called "steps" it has to mean steps, for instance in r.getIntegrator.setSetting(steps, 50)

If you want to provide functionality for setting number of points, you should support an additional argument/setting number_of_points, which than set steps = number_of_points-1 (with check that >2).

Especially because other solvers understand steps also as (steps +1 == number_of_points), for instance Copasi. They call it intervals, but it is steps.

Is still broken with the integrator setting steps currently in the right direction ;). If you call

s = r.simulate(0,10)
s.shape

you will get 51 points! from steps=50

On Mon, Jan 18, 2016 at 6:07 PM, Kyle Medley notifications@github.com wrote:

This behavior was changed per @hsauro https://github.com/hsauro in

119

https://github.com/sys-bio/roadrunner/issues/119: Herbert specified that steps should be the number of points, not the number of intervals.

— Reply to this email directly or view it on GitHub <https://github.com/sys-bio/roadrunner/issues/261#issuecomment-172593041 .

Matthias König Email: konigmatt@googlemail.com Tel: + 49 30 20938450 Tel: + 49 176 81168480 https://www.livermetabolism.com/contact.html Humboldt-University Berlin, Institute for Theoretical Biology, Junior Group König Charité Berlin, Institute of Biochemistry, Computational Systems Biochemistry http://www.charite.de/sysbio/people/koenig/

— Reply to this email directly or view it on GitHub https://github.com/sys-bio/roadrunner/issues/261#issuecomment-172618533.

matthiaskoenig commented 8 years ago

Just found this in the SED-ML specification L1V2 under numberOfPoints :/ Nothing to do, just for your information.

"Software interpreting the uniformTimeCourse is expected to produce a first outputPoint at time outputStartTime and then numberOfPoints output points with the results of the simulation. Thus a total of numberOfPoints + 1 output points will be produced." No idea how this got in there and why it was not called steps.

So numberOfPoints in SED-ML are steps everywhere else.

0u812 commented 8 years ago

Wow, odd nomenclature. Thanks for the heads-up.