Upgrade to object-oriented ODE solver

[ankitaggarwal011]

The current ODE solver odeint is inflexible and suffers from problem of stiffness as explained here in some cases.

It is better to upgrade to ODE solver object-oriented ode, as it is more flexible and allows to set various options.

[pjbollinger]

I would like to look into this and see if I can replace the current ODE solver for the more flexible one.

[ankitaggarwal011]

@pjbollinger, please go ahead. Your contributions are most welcome. Thank you. Let me know if you need any assistance.

[pjbollinger]

@ankitaggarwal011 Thanks. I'm trying to integrate it but running into memory error. You can see it in my latest build on Travis CI https://travis-ci.org/pjbollinger/PyCNN/jobs/167522797. It doesn't always happen but it occurs when setting initial value of integrator. I did a quick search and found this post about using a banded jacobian or simpler ODE. What are your thoughts?

[ankitaggarwal011]

Good work @pjbollinger. While I need to try this on system myself, from the cursory reading I would like to point this out from here:

Note: The first two arguments of f(t, y, ...) are in the opposite order of the arguments in the system definition function used by scipy.integrate.odeint.

Are you taking care of this? I'll try this myself and update you at my earliest. Thank you for your contribution.

[pjbollinger]

Good catch. I will have to correct that and see the results.

[pjbollinger]

So I'm debugging the difference between scipy.integrate.ode and scipy.integrate.odeint. From the attached screenshot you can see that there is a difference in the object returned after integration. I will continue to investigate but this is harder than I initially thought. 😅

[ankitaggarwal011]

@pjbollinger, I think your doing an excellent job. Keep it up! I will try it out myself today and update you. Thank you.

[ankitaggarwal011]

@pjbollinger, I tried out your code with minor modifications and I'm getting good results. You can try out this code:

tf, dt = 10, 1
ode = sint.ode(self.f) \
    .set_integrator("vode") \
    .set_initial_value(z0, 0) \
    .set_f_params(Ib, Bu, tempA)
while ode.successful and ode.t < tf:
  ode_result = ode.integrate(ode.t + dt)
z = self.cnn(ode_result)
l = z[:].reshape((self.n, self.m))

The problems come when I try to use lsoda integrator (which automatically switchs between Adams and BDF method for non-stiff and stiff equations respectively). This is will require a banded jacobian function for stiff equations, as mentioned in this link that you shared.

Using this with lsoda integrator would require more investigation, but I think we can go ahead with vode integrator at this point. What do you think?

Any suggestions are welcome. Thanks.

[pjbollinger]

@ankitaggarwal011 I opened a pull request with the new addition. Thank you for the help. I think there is still more work to be done but this is a good step. 😄