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. 😄

[ankitaggarwal011]

@pjbollinger. You've done great work. Yes, it does need more work considering we're still facing an issue using lsoda solver. Did you get a chance to investigate on this?

By the way, I think we should keep this issue open even after merging your current request, else there is a chance that this might be forgotten. Thanks.

[pjbollinger]

I didn't get a chance to review using the lsoda solver and need to read a little more about how the different integrators work.

One thing that I don't understand is the purpose of t. The function is not dependent on t, but you can't integrate with any value of t. I've tried to use ode.integrate(10) but it does not pass the diagonal edge detection Lenna test. I then tried to use ode.integrate(0.2) and many of the other tests failed but diagonal edge detection test passed.

[ankitaggarwal011]

@pjbollinger, t basically refers to time points for integration. Some processing methods take longer or shorter number of time points to converge and produce results. Let me know if you've more doubts. Thanks.

[ankitaggarwal011]

Even though the current ODE solver has been replaced by the object-oriented solver, this issue still needs investigation of a few parts. So, I think we should keep this issue open until we can be sure about certain parameters. What do you think @pjbollinger?

[pjbollinger]

Do you think those issues can be identified and listed? Depending on the issue, we can probably keep work done in this issue, but if the list gets too big or needs a lot of investigation/discussion, then that list should be broken into separate issues to manage easier. What are you thinking @ankitaggarwal011?

[ankitaggarwal011]

@pjbollinger, at this time, I can only think of the problem of dealing with stiff equations. The lsoda solver automatically takes care of that, but since we're experiencing problem setting it up, this is one of the issues. Other than this, I cannot of think of any major issue in this context right now.

Solving the problem with setting up lsoda solver would probably close this issue.