Learn how easy it is to sync an existing GitHub or Google Code repo to a SourceForge project!

## Is PyDSTool right for me?

Help
Brian M
2014-06-19
2014-06-20

• Brian M
2014-06-19

I have been attempting to model a moderately simple system where there are different types of particles. Particles of a certain type can change type depending on certain rules which can be represented by ODEs. I do not model their position, only the number of particles of each type at some particular locations.

The main thing is that the total number of particles must be conserved.

I was initially using scipy.integrate.odeint, but ran into the issue I describe here: http://scicomp.stackexchange.com/questions/12938/scipy-odeint-sum-of-conservative-ode-equations-does-not-remain-zero-as-it-is-be

Would PyDSTool be able to help by allowing me to perform symplectic integration with it?

• Rob Clewley
2014-06-19

Hi Brian,

We don't have a symplectic integrator in PyDSTool, unfortunately. Are you sure there is not a way to replace the definition of one of your ODEs with an algebraic condition on all the others that ensures the whole system satisfies the constraint that you want to conserve?

e.g.
dx/dt = f(x,y,z)
dy/dt = g(x,y,z)
z = N - x - y

And, does the total number of particles actually have to be an integer, or is it more like density of particles as a real number? It would be helpful to see some of the ODEs written out to be able to help more.

• Brian M
2014-06-19

Hi Rob,

Thanks for taking out the time to help.

Let us say there are 3 "nodes" (points where we are tracking the number of particles). Let us say that we have particles that can convert between two species: Pa and Pb. Pb can further interact with a "global pool" (pool) accessible to all nodes.

Let the rate at which species a becomes b be: k_ab.
Let the rate at which species b becomes a be: k_ba.
Let the rate at which species b goes into the pool be: k_bp.
Let the rate at which particles from the pool become species b at a node be: k_pb.
Let the rate at which particles "diffuse" between nodes be: D

Then, for this simple system, the following differential equations describe interactions (the numbers after the underscores represent the node number, . represents multiplication because the usual asterisk was being interpreted as the italics decorator):

d Pa_1/dt = -k_ab.Pa_1 + k_ba.Pb_1 + D.(Pa_2 + Pa_3 - 2.Pa_1)
d Pb_1/dt = k_ab.Pa_1 - k_ba.Pb_1 + D.(Pb_2 + Pb_3 - 2.Pa_1) - k_bp.Pb_1 + k_pb.pool

Similar equations are written for the other nodes. In general, if we have three nodes i, j, k and we are looking at node i, the diffusion term is for species x (where x can be a or b):

D.(Px_j + Px_k - 2.Px_i)

The following equation describes the time evolution of the pool:

d pool/dt = k_bp.Pb_1 + k_bp.Pb_2 + kbp_3.Pb_3 - 3.kpb.pool

Let the total number of particles be P = Pa_1 + Pa_2 + Pa_3 + Pb_1 + Pb_2 + Pb_3 + pool. Importantly:

d P/dt = 0,

which one can confirm the above equations satisfy.

The total number of particles doesn't have to be an integer, and is not always an integer in my implementation either. However, the odes should still sum to zero, or, the the total number of particles should remain constant, which does not happen if naively use odeint.

After a couple of days of poring over my code in order to find a bug, I am fairly certain that I don't have a mistake in the definition of the odes themselves. I found that changing odeint's error tolerance values (atol, rtol) changed the answer somewhat, but not necessarily significantly. However, I found that writing something like this:

```initialValues = np.array([.....])

''' 5000 steps because 5000*0.0001 = 0.5 '''
for x in xrange(5000):
result = odeint(deriv_function, initial_values, [0, 0.0001])
initial_values = result[-1]
```

conserved things quite well (even though there is a truly massive speed reduction) compared to:

```odeint(derivative_function, initialValues, [0, 0.5])
```

or

```odeint(derivative_function, initial_values, np.linspace(0, 0.5, 5000))
```

So I am now quite certain that the issue is indeed numerical.

The other option I am going to consider, based on someone else's suggestion is taking a symbolic route.

Last edit: Brian M 2014-06-19

• Rob Clewley
2014-06-19

Hi Brian, as the respondents on stack exchange explained, it is almost certainly unavoidable numerical discretization error. However, your type of conservation is a very easy kind to deal with, and from a brief look it seems to fall exactly into the case that I suggested you could solve with an algebraic constraint.

The simplest way to do this is to pick one variable, remove its ODE, and add a PyDSTool auxiliary variable defined as P-(sum of P_i), where you'll have to write out that sum longhand. The resulting output from a simulation will include that coordinate, although it will have a different status than the others. So, if you want to do certain analyses, such as bifurcation analysis, on that variable, you'd have to switch which variable is the algebraic one.

Making an aux var is simple in PyDST. Just define that variable (call it x here) along with the others in the usual way (see the tutorials), then add a declaration to the FuncSpec arguments:

DSargs.auxvars = ['x']

Follow a simple tutorial example on the web to get started before you try to understand this step.

A slightly more sophisticated solution is to use the Radau DAE solver, provided that a C compiler is set up. See DAE_example.py in the tests/ directory (or examples/ in the github version, if you use that).

-Rob

• Brian M
2014-06-19

Hi Rob. Let me give that a whirl. I hope to report back with success.

• Brian M
2014-06-19

This is just for fun, but I tried running the DAE_example.py. I am fairly certain I setup things correctly (SWIG etc.), and I am getting the following error:

```Error occurred in generating Radau system...
<type 'exceptions.SystemExit'> error: [Error 2] The system cannot find the file specified

Traceback (most recent call last):
File "C:\PyDSTool\tests\DAE_example.py", line 41, in <module>
File "C:\PyDSTool\Generator\Radau_ODEsystem.py", line 550, in __init__
self.compileLib()
File "C:\PyDSTool\Generator\Radau_ODEsystem.py", line 940, in compileLib
raise RuntimeError
RuntimeError
```

• Rob Clewley
2014-06-19

Well, it's likely that you don't quite have it set up right, unless maybe you're on a 64 bit system. You should read the recent discussions here that include a lot of dialogue about setting up SWIG, etc. and troubleshooting.

If you get that set up then you'll get a speed up of 10x or more from using Dopri/Radau for solving your equations. You may want to adjust the step size for Vode in the meantime, in case you can speed it up that way. The default error tolerances are very small so those can probably also be relaxed compared to what you've been using. I'm glad you got the basic setup working.

• Brian M
2014-06-19

I am on a 64-bit system, but I am not using 64-bit Python or anything...

I'll give the other thread a good look through. Profiling has shown that most of my time is indeed spent in the solver.

• Brian M
2014-06-19

Hi Rob! Good news again: I got the functions that user external integrators to work! Yay!

EDIT: actually, the later tests had this error:

```Traceback (most recent call last):
File "C:\PyDSTool\tests\DAE_example.py", line 41, in <module>
File "C:\PyDSTool\Generator\Radau_ODEsystem.py", line 550, in __init__
self.compileLib()
File "C:\PyDSTool\Generator\Radau_ODEsystem.py", line 937, in compileLib
rout.stop()
File "C:\PyDSTool\Redirector.py", line 82, in stop
os.remove(self.tmpfn)
WindowsError: [Error 32] The process cannot access the file because it is being used by another process: 'c:\\users\\pkbs\\appdata\\local\\temp\\tmploz9hr.pyout'
Traceback (most recent call last):
File "C:\PyDSTool\tests\DAE_example.py", line 41, in <module>
File "C:\PyDSTool\Generator\Radau_ODEsystem.py", line 550, in __init__
self.compileLib()
File "C:\PyDSTool\Generator\Radau_ODEsystem.py", line 937, in compileLib
rout.stop()
File "C:\PyDSTool\Redirector.py", line 82, in stop
os.remove(self.tmpfn)
WindowsError: [Error 32] The process cannot access the file because it is being used by another process: 'c:\\users\\pkbs\\appdata\\local\\temp\\tmploz9hr.pyout'
```

Last edit: Brian M 2014-06-19

• Brian M
2014-06-20

Hi Rob, I haven't fixed the Radau error, but I am not getting better performance with Dopri. The opposite in fact: marginally worse than VODE and it gives me an error

~~~
integrator/memory.c: Assertion failed '

initIntegData == 0
~~~

Last edit: Brian M 2014-06-20

• Rob Clewley
2014-06-20

When the C integrators are correctly set up they are almost always a lot quicker than VODE. Your problem may be stiff, which would benefit more from Radau, or you may not have the algorithmic params set up well. Something with you setup is unhealthy if you've got a memory error from the start. I suggest you email me your script privately (click on my name in this thread and then click on "send message" in top right) and I can check it over.

• Rob Clewley
2014-06-20

When the C integrators are correctly set up they are almost always a lot quicker than VODE. Your problem may be stiff, which would benefit more from Radau, or you may not have the algorithmic params set up well. Something with you setup is unhealthy if you've got a memory error from the start. I suggest you email me your script privately (click on my name in this thread and then click on "send message" in top right) and I can check it over.

• Brian M
2014-06-19

Hi Rob. I just wanted to report that I was able to successfully use PyDSTool as we discussed -- and indeed, the computation is correct (everything sums up to the right value in the end)!

There is a slight increase in the time it takes to run the script (1.25x longer than before), but that is a small price to pay for correctness, no? Plus, I am sure after I spend some time profiling, I will find that there are little things I can clean up on my end.

Anyway, thanks very much for your help!

A small hitch I encountered (that I quickly circumvented), is that I could not refer to the auxiliary variables in the RHS of the other equations. It was easy to circumvent this problem by simply writing out what the auxiliary variables stand for where I had to refer to them...but...yeah. It would have been simpler to just be able to refer to the auxiliaries themselves?

Last edit: Brian M 2014-06-19

• Rob Clewley
2014-06-20

That's the reason you'd want to use Radau's DAE solver, so that all the variables are primary variables and can be referenced in each other's RHS. The quick way I showed you computes aux vars at the end of the simulation. That's why they are auxiliary, not primary. Alternatively, you can still tweak this by adding a "reused term" definition like "x_temp" which holds the definition of your final variable. You would still use "x" as the aux var, and use "x_temp" in your RHS defs to cross-reference without repeatedly re-computing the sum. That temp value is computed once at the top of the vector field function, so it's more efficient. Reused terms are explained in the docs.