[SciPy-Dev] Consensus on the future of integrate.ode

[Geoff Oxberry]

On Mon, Sep 9, 2013 at 5:53 AM, Benny Malengier
<benny.malengier at gmail.com>wrote:

>
>
>
> 2013/9/9 Benny Malengier <benny.malengier at gmail.com>
>
>>
>>
>>
>> 2013/9/8 boyfarrell at gmail.com <boyfarrell at gmail.com>
>>
>> Hello Juan Luis,
>>>
>>> Yes, I saw your pull request and this, in part, motivated me to start
>>> this discussion. Looks like you put a lot of work into your odeintmodifications!
>>>
>>> I completely agree with you, that ode should be fixed first, then odeintshould be a nice interface wrapped around it.
>>> integrate.ode becomes a package aimed at people who understand the
>>> complexities of solving ODE problems and integrate.odeint is for people
>>> who just want an answer.
>>>
>>> Also, thank you for doing a nice summary, that was really helpful.
>>>
>>> Remark 1: Backend.
>>>
>>> I don't think we need to reinvent the wheel, Benny's work looks
>>> excellent so I think it make sense to base a new integrate.ode off
>>> scikits.odes. It has a nice array of modern solvers, in both C and
>>> Fortran.
>>>
>>
>> Thanks for the thumbs up.
>> However, odes is only one of 3 implementations, and then the one that
>> exposes least of the C solvers. Odes was based originally on the previous
>> pysundials which was no longer maintained and was a non-cython wrapper.
>> The other implementations are geared towards a specific problem domain
>> though.
>>
>> One is: https://github.com/casadi/casadi/wiki (LGPL tough)
>>
>> The other: http://www.jmodelica.org/assimulo (annoying copyright
>> transfer contribution license though to a company,
>> http://www.jmodelica.org/page/14)
>>
>> The non sundials solvers in odes are also not as state of the art as
>> some that are added in above two interfaces.
>> The problem with above packages is their license and the fact that they
>> are packages with parsing language, ..., see eg:
>> https://github.com/casadi/casadi/blob/master/examples/python/dae_single_shooting_in_30_lines.py.
>>
>> Note that you can pass an array of times to sundials at which you want
>> output, and then all iteration is done inside sundials, so the point of
>> time step outside of sundials is not needed (though in many problems you'll
>> want to do stepping controlled in a python loop).
>>
>>
> What would really be interesting is find some way to use the parallel
> implementation of sundals in python. Some mapping of a numpy array to the
> parallel sundial vector (
> https://github.com/bmcage/odes/blob/master/scikits/odes/sundials/nvector/nvector_parallel.h)
> would then be needed. Has there ever been work to contruct a parallel
> aware array?
>

petsc4py does this already, using data structures from PETSc. You can
instantiate some of these data structures with NumPy arrays.

I'm not sure if such work has already been done in SciPy, nor am I sure if
SciPy is "the right place" for it. I don't know a ton about parallel
computing, but my limited experience with it backs up the common warning
people give about parallelizing code -- you really have to parallelize your
data structures in order for most algorithms to be effective. My worry is
that people will write SciPy code as they have done for the previous
releases -- that is, in serial -- and then expect parallel ODE
functionality to just work, when continuing to use serial data structures
won't let you take advantage of parallelism.

The other thing about wrapping the parallel implementation of SUNDIALS I'd
be concerned about is the preconditioning of iterative linear solvers; in
most cases, these solvers have to be preconditioned in order to converge
quickly. For that, do you expect users to provide their own
preconditioners, like
SUNDIALS does?

Something that is useful about solvers like PETSc/petsc4py or
Trilinos/PyTrilinos is that they provide a number of methods that can also
be used to precondition iterative linear solvers (variants of incomplete
LU/Cholesky factorization, block Jacobi, algebraic multigrid, and others),
which makes it easier for users to select from a number of "black box"
methods in the event that they cannot easily construct a
preconditioner themselves.
These sorts of things would probably be better implemented in other parts
of SciPy, since they could be used for other purposes besides solving ODEs
(namely, for any application involving the solution of a linear system
using iterative methods).

Geoff



> Would be nice to find funding for that somewhere
>
> Benny
>
>
>
>> Benny
>>
>>
>>> Remark 2 Frontend.
>>>
>>> We can define a new integrate.ode interface that wraps scikits.odes. It
>>> seems that you already have some ideas in that area, mentioning the MATLAB
>>> style. I looked at the MATLAB docs just now and spend 20 minutes thinking
>>> how this might look if it were more pythonic. You can take a look here,
>>> fork your own version to make changes or add comments,
>>>
>>> https://gist.github.com/danieljfarrell/6482713
>>>
>>> Best wishes,
>>>
>>> Dan
>>> _______________________________________________
>>> SciPy-Dev mailing list
>>> SciPy-Dev at scipy.org
>>> http://mail.scipy.org/mailman/listinfo/scipy-dev
>>>
>>
>>
>
> _______________________________________________
> SciPy-Dev mailing list
> SciPy-Dev at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-dev
>
>


-- 
Geoffrey Oxberry, Ph.D., E.I.T.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/scipy-dev/attachments/20130909/14f2313a/attachment.html>