Send SciPy-Dev mailing list submissions to

scipy-dev <at> scipy.orgTo subscribe or unsubscribe via the World Wide Web, visit

https://mail.scipy.org/mailman/listinfo/scipy-devor, via email, send a message with subject or body 'help' to

scipy-dev-request <at> scipy.orgYou can reach the person managing the list at

scipy-dev-owner <at> scipy.orgWhen replying, please edit your Subject line so it is more specific

than "Re: Contents of SciPy-Dev digest..."

Today's Topics:

1. Re: ode and odeint callable parameters (Marcel Oliver)

2. Re: ode and odeint callable parameters (Benny Malengier)

3. Re: ode and odeint callable parameters (Marcel Oliver)

----------------------------------------------------------------------

Message: 1

Date: Mon, 25 Jan 2016 14:29:20 +0100

From: Marcel Oliver <

m.oliver <at> jacobs-university.de>

To: SciPy Developers List <

scipy-dev <at> scipy.org>

Subject: Re: [SciPy-Dev] ode and odeint callable parameters

Message-ID: <

22182.9008.344108.842230 <at> xps13.localdomain>

Content-Type: text/plain; charset="us-ascii"

Anne Archibald writes:

> There is periodically discussion of the mess of interfaces and

> solvers for ODEs in scipy (see the archives about six months ago,

> for example). One of the concerns is that people want to do very

> different things, so settling on a good interface is not at all

> easy. I don't just mean the underlying algorithms, but also the

> calling interface. It's easy to support someone who drops in an RHS

> and asks for the solution at a hundred predefined points. It's also

> clear that a general toolkit shouldn't support my use case: a

> 22-dimensional solver with a compiled RHS generated by a symbolic

> algebra package, with internal precision switchable between 80 and

> 128 bits, with a root-finder attached to the output to define

> stopping places, that needs to all run with compiled speed. So

> where along that scale do you aim scipy's interface? I have my own

> ideas (see aforementioned archive thread) but don't have the energy

> to implement it myself. And anyway, some people need very different

> things from their ODE solvers (for example, solution objects

> evaluable anywhere).

>

> Anne

Thanks for all the comments. I had not been aware of the existence of

scikits odes - the additional capabilities are very nice, but API-wise

it just adds to the mess...

Being unsure what has previously been discussed, I'd just like to

share some random thoughts, as I don't claim to understand the problem

domain completely. Basically, the API should make simple thing simple

and natural, and complicated things possible. So I think your

requirements stated above should not be completely out of reach - in

principle - even though high precision arithmetic would require

completely new backends and whether one can get close to compiled

speed is probably very problem dependent.

So, here an unordered list of thoughts:

1. I don't think there is a big difference in practice between the

ode, odes, and odeint interfaces. Personally, think that

odeint (f, y0, t)

is just fine. But the object oriented approach of ode/odes is

good, too; however, the various parameters should be controlled by

(keyword) arguments as a statement like

r.set_initial_value(y0, t0).set_f_params(2.0).set_jac_params(2.0)

is just not pretty.

So I like the odes approach best, with a basic object oriented

structure, but the parameters set as arguments, not via methods.

Not sure if odes allows to set an initial time different from zero,

that is essential. I also dislike that specifying the solver is

mandatory and first in the argument list; I think this is

specialist territory as many of the existing solvers perform

sufficiently well over a large range of problems so that changing

the solver is more of an exception.

If the interface is object oriented (as in ode or odes), then at

least the default solver should be reentrant, otherwise the

API raises false expectations.

I would also suggest that, independent of the rest of the API, the

final time t can be an array or a scalar; the result is either

returned as a vector of shape y0 when t is scalar and as an array

with one additional axis for time when t is a 1-d-array.

2. It would be very nice to control output time intervals or exit

times. Mathematica has a fairly rich "WhenEvent" API, although it

seems to me it may also be no more than a fancy wrapper around a

standard solver.

The Mathematica API feels almost a bit too general, but there are

certainly two independent ways triggers could be implemented, both

seem useful to me:

(a) Simple Boolean trigger for "output event" and "termination

event" based on the evaluation of some function of (y,t,p)

which is passed to the solver and which is evaluated after each

time step.

(b) Trigger for "output event" and/or "termination" based on

finding the root of a scalar function of the (y,t,p). This

needs either a separate root finder or a modification to the

implicit solver (the latter may be hard as that would require

nontrivial modification of the embedded solver code).

(odes "rootfn" may be doing just that, but there seems to be no

detailed documentation...)

3. Returning interpolating function objects as, e.g., Mathematica

does, seems a bit of overkill to me. It would be easier to set up

an independent general interpolating function library (as possibly

even exists already) and pass a discrete time series to it. It is

likely that by deriving the interpolation function directly from

the interpolation on which the solver is based one could gain some

efficiency and give better control of the interpolation accuracy,

but to me that seems to require a lot of effort for relatively

little improvement.

4. odes includes DAE solvers. I wonder if it would not be more

natural to pass the algebraic part of the equation as a separate

keyword argument and let the solver wrapper decide, based on the

presence or absence of this argument, whether to call a DAE

backend.

5. As I wrote before, the API should deal transparently with

d-dimensional vector fields.

Just a question about solver backends: it seems that the SUNDIALS

solvers are particularly strong for large stiff systems which come

from parabolic PDEs? (And of course the ability to propagate

sensitivities, is that ability made available through odes, too?)

On smaller system of ODEs, I have found that ode "dop853" is extremely

accurate and seems capable of dealing with rather stiff problems even

though I have no clue how it does it. (E.g., for the Van der Pol

oscillator with large stiffness parameter where other explicit solvers

either blow up or grind to a halt, dop853 and also dopri5 are doing

just fine while Matlab's dopri5 solver fails as expected. It's a

mystery to me.)

In any case, it would be very nice to develop a clear plan for

cleaning up ODE solver backend in Scipy. Once there is a clear

target, it might be easier to see what is easy to do and what is

missing on the backend side.

Regards,

Marcel

------------------------------

Message: 2

Date: Mon, 25 Jan 2016 15:22:09 +0100

From: Benny Malengier <

benny.malengier <at> gmail.com>

To: SciPy Developers List <

scipy-dev <at> scipy.org>

Subject: Re: [SciPy-Dev] ode and odeint callable parameters

Message-ID:

<

CAOyNq8qd-4FrS4QDn861F0bSwsOaEwx9sHs0oAwCENLOhYQTGA <at> mail.gmail.com>

Content-Type: text/plain; charset="utf-8"

2016-01-25 14:29 GMT+01:00 Marcel Oliver <

m.oliver <at> jacobs-university.de>:

> Anne Archibald writes:

> > <snip>

> > to implement it myself. And anyway, some people need very different

> > things from their ODE solvers (for example, solution objects

> > evaluable anywhere).

> >

> > Anne

>

> <snip>

>

> 2. It would be very nice to control output time intervals or exit

> times. Mathematica has a fairly rich "WhenEvent" API, although it

> seems to me it may also be no more than a fancy wrapper around a

> standard solver.

>

> The Mathematica API feels almost a bit too general, but there are

> certainly two independent ways triggers could be implemented, both

> seem useful to me:

>

> (a) Simple Boolean trigger for "output event" and "termination

> event" based on the evaluation of some function of (y,t,p)

> which is passed to the solver and which is evaluated after each

> time step.

>

> (b) Trigger for "output event" and/or "termination" based on

> finding the root of a scalar function of the (y,t,p). This

> needs either a separate root finder or a modification to the

> implicit solver (the latter may be hard as that would require

> nontrivial modification of the embedded solver code).

>

> (odes "rootfn" may be doing just that, but there seems to be no

> detailed documentation...)

>

Yes rootfn does that.

A nicer interface with onroot functionality to that was added by a

contributor to CVODE, I just extended it to the IDA solver. Only available

through the solve method of the solvers (as STEP methods just stop),

example:

https://github.com/bmcage/odes/blob/master/docs/src/examples/ode/freefall.pySame with fixed tstop added in a list of output times:

https://github.com/bmcage/odes/blob/master/docs/src/examples/ode/freefall_tstop.py>

> 3. Returning interpolating function objects as, e.g., Mathematica

> does, seems a bit of overkill to me. It would be easier to set up

> an independent general interpolating function library (as possibly

> even exists already) and pass a discrete time series to it. It is

> likely that by deriving the interpolation function directly from

> the interpolation on which the solver is based one could gain some

> efficiency and give better control of the interpolation accuracy,

> but to me that seems to require a lot of effort for relatively

> little improvement.

>

agree. Also, just reinit the solver and solve again from closest start time

is an option too if for some reason interpolation error must be avoided.

CVODE can actually give you output backward if you have jumped too far

ahead, but that requires you to know where you want interpolation output

immediately after doing a STEP. All very problem depending.

>

> 4. odes includes DAE solvers. I wonder if it would not be more

> natural to pass the algebraic part of the equation as a separate

> keyword argument and let the solver wrapper decide, based on the

> presence or absence of this argument, whether to call a DAE

> backend.

>

These things all require a lot of boilerplate code. In the end, my

experience is that the documentation of the original solvers is the best to

learn the details, so keeping variable names the same in options is more

important than something that is more pythonic. Eg, the rootfn above is a

bad name, but fits the original documentation at

https://computation.llnl.gov/casc/sundials/documentation/documentation.htmlCVODE doc is 164 pages, in the end, for problems you want to optimize you

will grab that documentation to learn what options you want to tweak.

About transparant algebraic part. The backends need a lot of info on the

algebraic part. You also need to structure your variables specifically if

you eg want banded Jacobian. And the calling sequence of all functions

becomes different with an extra parameter in and out (ydot).

Quite a lot of boilerplate that slows your wrapper down will be needed I'm

afraid.

> 5. As I wrote before, the API should deal transparently with

> d-dimensional vector fields.

>

I'm not really following. You mean unwrap and wrap again like the

complex_ode solver does now in scipy?

http://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.integrate.complex_ode.htmlWriting your own wrapper for your problem does not seem like a big task to

me.

>

> Just a question about solver backends: it seems that the SUNDIALS

> solvers are particularly strong for large stiff systems which come

> from parabolic PDEs? (And of course the ability to propagate

> sensitivities, is that ability made available through odes, too?)

>

They are used for that too, but it is not the core. Large complex kinetic

systems with different time scales seems like the largest use case on the

mailing list, see also as first example in the example doc.

Sensitivities are not present in odes at the moment, but I think some

wrappers have that covered already. Let's say it is asked once or twice a

year to add.

> On smaller system of ODEs, I have found that ode "dop853" is extremely

> accurate and seems capable of dealing with rather stiff problems even

> though I have no clue how it does it. (E.g., for the Van der Pol

> oscillator with large stiffness parameter where other explicit solvers

> either blow up or grind to a halt, dop853 and also dopri5 are doing

> just fine while Matlab's dopri5 solver fails as expected. It's a

> mystery to me.)

>

Some extra tests seem in order. How does method=adams in ode succeed?

Should you be able to quickly make an ipython notebook ...

> In any case, it would be very nice to develop a clear plan for

> cleaning up ODE solver backend in Scipy. Once there is a clear

> target, it might be easier to see what is easy to do and what is

> missing on the backend side.

>

Yes. We don't have our own backends, we only wrap backends. In the end,

those backends with continued development will win out. Wrapping solvers is

a tricky business, as it is very hard to avoid design of the solver to not

slip through

>

> Regards,

> Marcel

> _______________________________________________

> SciPy-Dev mailing list

>

SciPy-Dev <at> scipy.org>

https://mail.scipy.org/mailman/listinfo/scipy-dev>

-------------- next part --------------

An HTML attachment was scrubbed...

URL: <

https://mail.scipy.org/pipermail/scipy-dev/attachments/20160125/e34dbc3b/attachment-0001.html>

------------------------------

Message: 3

Date: Tue, 26 Jan 2016 11:43:26 +0100

From: Marcel Oliver <

m.oliver <at> jacobs-university.de>

To: SciPy Developers List <

scipy-dev <at> scipy.org>

Subject: Re: [SciPy-Dev] ode and odeint callable parameters

Message-ID: <

22183.19918.593905.530551 <at> xps13.localdomain>

Content-Type: text/plain; charset="us-ascii"

Benny Malengier writes:

[...]

> Yes rootfn does that. A nicer interface with onroot functionality

> to that was added by a contributor to CVODE, I just extended it to

> the IDA solver. Only available through the solve method of the

> solvers (as STEP methods just stop), example:

>

https://github.com/bmcage/odes/blob/master/docs/src/examples/ode/freefall.py > Same with fixed tstop added in a list of output times:

>

https://github.com/bmcage/odes/blob/master/docs/src/examples/ode/freefall_tstop.pyThanks, I'll be looking at that!

(Side remark: tried to compile odes today on a Fedora 23 machine, but

it fails with

gcc: error: /usr/lib/rpm/redhat/redhat-hardened-cc1: No such file or directory

gcc: error: /usr/lib/rpm/redhat/redhat-hardened-cc1: No such file or directory

error: Command "gcc -pthread -Wno-unused-result -DDYNAMIC_ANNOTATIONS_ENABLED=1 -DNDEBUG -O2 -g -pipe -Wall -Werror=format-security -Wp,-D_FORTIFY_SOURCE=2 -fexceptions -fstack-protector-strong --param=ssp-buffer-size=4 -grecord-gcc-switches -specs=/usr/lib/rpm/redhat/redhat-hardened-cc1 -m64 -mtune=generic -D_GNU_SOURCE -fPIC -fwrapv -fPIC -Ibuild/src.linux-x86_64-3.4 -I/usr/lib64/python3.4/site-packages/numpy/core/include -I/usr/include/python3.4m -c build/src.linux-x86_64-3.4/fortranobject.c -o build/temp.linux-x86_64-3.4/build/src.linux-x86_64-3.4/fortranobject.o" failed with exit status 1

Seems there are some hard-coded paths somewhere, did not have time to

track this down other than noting that I get this error, so I may have

done something stupid - don't worry if this is the case.)

[...]

> 4. odes includes DAE solvers. I wonder if it would not be more

> natural to pass the algebraic part of the equation as a separate

> keyword argument and let the solver wrapper decide, based on the

> presence or absence of this argument, whether to call a DAE

> backend.

>

> These things all require a lot of boilerplate code. In the end, my experience

> is that the documentation of the original solvers is the best to learn the

> details, so keeping variable names the same in options is more important than

> something that is more pythonic. Eg, the rootfn above is a bad name, but fits

> the original documentation at

>

https://computation.llnl.gov/casc/sundials/documentation/documentation.html > CVODE doc is 164 pages, in the end, for problems you want to optimize you will

> grab that documentation to learn what options you want to tweak.

> About transparant algebraic part. The backends need a lot of info on the

> algebraic part. You also need to structure your variables specifically if you

> eg want banded Jacobian. And the calling sequence of all functions becomes

> different with an extra parameter in and out (ydot).

> Quite a lot of boilerplate that slows your wrapper down will be needed I'm

> afraid.

I see the point. But then one can do this only for the default solver

suite (I assume you'd go for sundials) and if it's necessary to

replace the backend, then the non-default solver should get the

wrapper code to make this transparent.

> 5. As I wrote before, the API should deal transparently with

> d-dimensional vector fields.

>

> I'm not really following. You mean unwrap and wrap again like the

> complex_ode solver does now in scipy?

>

http://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.integrate.complex_ode.html >

> Writing your own wrapper for your problem does not seem like a big

> task to me.

I have done this many times, but if feels like something the computer

should be doing for me. I am coming from the user perspective, not as

a scipy developer, so I don't want the ugliness in MY code...

The other question is that best practice is not clear to me. What do

I need to do to avoid extra copies of the data, which should be

possible in general. A harder problem is probably what to do e.g. for

systems of reaction diffusion equations where one would want to keep

the sparse structure. Not sure if this is even possible to address in

a sufficiently general way.

> On smaller system of ODEs, I have found that ode "dop853" is extremely

> accurate and seems capable of dealing with rather stiff problems even

> though I have no clue how it does it. (E.g., for the Van der Pol

> oscillator with large stiffness parameter where other explicit solvers

> either blow up or grind to a halt, dop853 and also dopri5 are doing

> just fine while Matlab's dopri5 solver fails as expected. It's a

> mystery to me.)

>

> Some extra tests seem in order. How does method=adams in ode succeed? Should

> you be able to quickly make an ipython notebook ...

I attach a small code snippet. In fact, my memory failed me in that

dop853 indeed complains about insufficient nmax (nsteps in the Python

API). What breaks is the accuracy of the solution, but it's still

remarkably good and can be controlled by nsteps. It's obvious that

one should not solve this kind of problem with dop853, but I am still

amazed that it does not fall flat on its face.

In fact, vode/Adams dies completely (as it should) and vode/BDF gives

lots of warnings and computes a completely wrong solution (BDF schemes

are supposed to work here). lsoda is fine, though.

Best,

Marcel

PS.: My motivation for looking at the solver API comes from a nascent

book project where the idea is to include computer experiments on

nonlinear dynamical systems done in Scientific Python. Thus, I am

exploring best practice examples that are generic for various problem

domains.

-------------- next part --------------

A non-text attachment was scrubbed...

Name: prob4.py

Type: application/octet-stream

Size: 1032 bytes

Desc: not available

URL: <

https://mail.scipy.org/pipermail/scipy-dev/attachments/20160126/ae6bd65f/attachment.obj>

------------------------------

Subject: Digest Footer

_______________________________________________

SciPy-Dev mailing list

SciPy-Dev <at> scipy.orghttps://mail.scipy.org/mailman/listinfo/scipy-dev------------------------------

End of SciPy-Dev Digest, Vol 147, Issue 20

******************************************