## Specifying systems of equations

Help
2010-10-02
2013-11-04
• Eric Pedersen - 2010-10-02

Is there any way to include  external scipy python functions when defining systems of equations in PyCont? I want to make one of the state variables in my continuation a function of a spline, computed using the UnivariateSpline from sci.interpolate.
The code I'm using to define the system looks like this:

pred_spline = scipy.interpolate.UnivariateSpline(x,y,s=0)

DSargs.varspecs = { 'N': '(delta * (Nin - N) - (pred_spline(N) * D))',
'D': '(pred_spline(N) * D) - (delta * D))' }

However, all I get is this error:

Undeclared or illegal token `scipy` in spec string `pred_spline`

I'm not sure how to proceed from here.

Regards,
Eric

• Rob Clewley - 2010-10-02

Please email me your entire script so that I can reproduce the problem and explore whether I can resolve it for you.

• Rob Clewley - 2010-10-05

Yes, for continuing equilbria points only, the solution is use the less well-known functionality of inserting your own custom code/calls into a right-hand side function. First, make sure you download the latest updates to PyDSTool from here: http://www2.gsu.edu/~matrhc/CodeTopics.html. (For now the SVN repository is broken.)

Then, you basically add two arguments to your DSargs to declare that you are referring to a special function:

```DSargs.varspecs = { 'N': '(delta * (Nin - N) - (pred_spline_val * D))',
'D': '(pred_spline_val * D) - (( ( (bB * D)/(Kb + D)) * B / epsilon)) - (delta * D)',
'S': '((((bB * D)/(Kb + D)) * S ) - (delta + mort + lambd) * S)',
'B': '((((bB * D)/(Kb + D)) * S ) - (delta + mort) * B)'}
```

This version of your RHS uses the special name 'pred_spline_val' which is not yet declared. Now add:

```DSargs.vfcodeinsert_start = 'pred_spline_val = ds.pred_spline(N)'
DSargs.ignorespecial = ['pred_spline_val']
```

which will insert the definition of 'pred_spline_val' at the beginning of the python spec function for the RHS, and also tells the parser to ignore instances of that definition in the body of the varspecs. Notice that instead of 'self' you use 'ds' to refer to the final Generator object. Lastly, you need to attach a function to the Generator so that ds.pred_spline refers to something:

```from scipy.interpolate import UnivariateSpline
x = linspace(0,80,20)
y = 3.3*x/(4.3+x)
pred_spline = scipy.interpolate.UnivariateSpline(x,y,s=0)
ode  = PyDSTool.Generator.Vode_ODEsystem(DSargs)
ode.pred_spline = pred_spline
```

Unfortunately, this solution won't work for continuing periodic orbits as that requires an RHS to be built in pure C code. This is normally done automatically by PyCont but will fail with these additional definitions. You can still use a code insert but you'll have to use the 'nobuild'=False argument to DSargs and then add C code to make a spline (e.g., from a C library for splines) and it won't be identical to the scipy code. There's no way to call back to the python level here, that's something I wanted to avoid for the sake of speed, but in principle this feature could be added in the future, with the understanding of the performance hit involved.

See tests/HH_model_Cintegrator for an example of using 'nobuild' and then completing the object with the makeLib method call (where in between creating the incomplete Generator and making this call you would go and edit the .c vector field file by hand to insert your spline call and add your spline library code to the include list in the call).

• Anthony Anderson - 2013-07-03

I'm having trouble adapting this example to a similar situation where I've defined a function. Is there existing documentation on using custom code calls in Vode?

Suppose my model has state variables (x, y) and parameters (p,q) and I've defined a custom function

```def my_function(x,  y, p, q):
...
```

which has some complicated structure. Following the example above, I would specify the problem as follows:

```pars = {'p': 1, 'q': 2}
ics = {'x': 0, 'y': 1}
xdot = 'my_function_value'
ydot = 'p*my_function_value**2 - q'
dsargs = dst.args(name='example')
dsargs.pars = pars
dsargs.tdata = [0, 1]
dsargs.ics = ics
dsargs.varspecs = {'x': xdot, 'y': ydot}
dsargs.vfcodeinsert_start = 'my_function_value = my_function(x, y, p, q)'
dsargs.ignorespecial = ['my_function_value']
```

How should I attach the function to the Generator?

Also…  How would I modify vfcodeinsert_start if I have multiple custom functions? Is there a schema for using custom functions with Events?

Thanks for you help

• Rob Clewley - 2013-07-03

Hi,

The answer to the first part of your post is in my earlier reply to a similar question. You should have

```dsargs.vfcodeinsert_start = 'my_function_value = ds.my_function(x, y, p, q)'
```

and then once you've built the Generator (as "ode", say) you attach the function with

```ode.my_function = my_function  # or any callable object
```

You can have as many such calls to as many custom functions as you like. The inserted string can be triple-quoted to allow multiple lines to be inserted, which can do anything (the ODE object will be available as object "ds"). You have to indent subsequent lines with four spaces. If I'm not remembering that right then try indenting with TAB, i.e. '\t'.

As for events, support for custom functions that has not yet been included in the main release, but as it happens I recently completed support for it and it is now available if you update PyDSTool to the latest Bazaar repository version. Information for obtaining this version can be found at: http://sourceforge.net/scm/?type=bzr&group_id=140858

Usage of this new code looks like this:

```Mev3to4 = Events.makeZeroCrossEvent('minf(2*V/widthNa - offsetNa)-m', 1,
{'name': 'ev_3to4',
'eventtol': 1e-4,
'precise': True,
'term': True}, varnames=['V', 'm'], parnames=['widthNa', 'offsetNa'],
extra_funcspec_args={'ignorespecial': ['minf'],
'codeinsert_start': "minf = ds._minf"}
)
Mev3to4._minf = my_custom_callable_object  # e.g., can be re-referenced in your Generator later
```

• Anthony Anderson - 2013-07-08

Thanks for your reply. I've successfully inserted a custom function into the RHS of my vector equations. However, I'm having some trouble doing the same with Events (I've updated to the latest version, rev 378, on the Bazaar repository).

As before, I want to detect events from a function of the form:

```def my_function(x, y, p, q):
...
return val
```

with state variables (x,y) and parameters (p,q). My event specification is as follows:

```event = dst.Events.makeZeroCrossEvent('my_val - val0', -1,
{'name': 'terminate',
'eventtol': 1e-4,
'precise': True,
'active': True,
'term': True},
varnames = ['x', 'y'],
parnames = ['p', 'q', 'val0'],
extra_funcspec_args = {'ignorespecial': ['my_val'],
'codeinsert_start':  'my_val = ds.my_func(x, y, p, q)'})
```

I get an error message: KeyError: 'vars dictionary not present'. Any thoughts?

Thanks in advance for any help.

• Rob Clewley - 2013-07-08

This happens when the original event definition string uses neither variable or input names explicitly. I didn't foresee this use case. Fortunately, I'm fairly sure you can fix this by including the following key-dictionary mapping inside the argument dictionary in the third position of the makeZeroCrossEvent, alongside 'name' etc.

```'vars': {'x': None, 'y': None}
```

etc. for any variables x, y, etc. that are used in your function. For this use case, this dictionary technically just has to be non-empty and keyed by at least one variable name used in your system.

I haven't tested this myself yet so let me know if there's a problem.

• Anthony Anderson - 2013-07-08

OK, I've changed the vars dictionary in the Events, giving:

```event = dst.Events.makeZeroCrossEvent('my_val - val0', -1,
{'name': 'terminate',
'vars': {'x': None, 'y': None},
'eventtol': 1e-4,
'precise': True,
'active': True,
'term': True},
varnames = ['x', 'y'],
parnames = ['p', 'q', 'val0'],
extra_funcspec_args = {'ignorespecial': ['my_val'],
'codeinsert_start': 'my_val = ds.my_function(x, y, p, q)'})
event.my_function = my_function
```

Running the script gives me the error:

```Traceback (most recent call last):
File "integrate.py", line 206, in <module>
traj = odes.compute('test')
File "/Users/Anthony/PyDSTool/Generator/Vode_ODEsystem.py", line 398, in compute
eventslist)
File "/Users/Anthony/PyDSTool/Events.py", line 254, in pollHighLevelEvents
parDict=parDict), eventlist)
File "/Users/Anthony/PyDSTool/Events.py", line 254, in <lambda>
parDict=parDict), eventlist)
File "/Users/Anthony/PyDSTool/Events.py", line 730, in __call__
self.fval = self._fn(varDict_temp, parDict)
File "<string>", line 2, in _f_terminate_ud
NameError: global name 'x' is not defined
```

I'm guessing this is related to the above issue? Perhaps I'm not correctly attaching the callable function to the events object? That aspect is slightly confusing to me.

• Rob Clewley - 2013-07-08

OK, sorry about that. It would save time if you could just email me your entire script so that I can track in detail how the 'vars' dictionary is being used internally in this case. I'll post the reply and fix here for others to see.

• Rob Clewley - 2013-07-08

Yes, I needed to do more to support the more general cases of these custom functions. I'd only tested it on a simple kind in my own code so far. I've updated Events.py that will run your code and I think it also covers the general case properly now. I've uploaded a fix to changeset 379: http://pydstool.bzr.sourceforge.net/bzr/pydstool/revision/379

• Anthony Anderson - 2013-07-12

Yep, that did it. The Event detection is working with custom functions. Thanks for your help!

I see what you mean about the speed of the Vode solver. Especially with 'precise' event detection, which slows the integration down by a factor of 10. There's obvious advantages in using interpolated events as with the C integrators. Unfortunately, I've not been able to get those working on my system (OS X 10.8). I wish I could contribute some ideas here, but at the moment, it's a little beyond me….

Keep up the good work. PyDSTool is a great addition to Python!

• Rob Clewley - 2013-07-12

Thanks. Your install problem is unfortunate as you could probably get the custom functions to work in the C code too, maybe with a little more help from me (I have never tried using them in events at the C level). If you open a different forum thread about your OS X 10.8 problems I might be able to help. It would at least be a helpful data point to know about. I'm pretty sure I know people who have it working on 10.8 so I can look into that.