commons-dev mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From Luc Maisonobe <>
Subject Re: [math] seeking for an implementation advice
Date Sun, 07 Mar 2010 15:41:36 GMT
Phil Steitz a écrit :
> Luc Maisonobe wrote:
>> Hello,
>> As you may have seen, I have implemented a new feature in [math] last
>> few days. It is an extended integrator for ODE that in addition to state
>> y() also provides its jacobians with respect to initial state dy(t)/dy0
>> and with respect to some ODE parameters dy(t)/dp.
>> The underlying algorithm is simply to extend the ODE problem with a few
>> variational equations, so if you start with an y(t) state that has
>> dimension 6 for example, and you have 3 parameters you will end up with
>> a compound state that has dimension 6 + (6 * 6) + (6 * 3) = 60. One the
>> state has been extended and the additional equations have been set up,
>> the problem is a classical Initial Value Problem whci can be handled by
>> all the existing integrators we have. So this feature consists mainly in
>> a bunch of adaptor classes to perform the state versus compound state
>> mapping. There are adaptors for the ODE problem, the step handlers, the
>> step interpolators and the events handlers.
>> For now, the implementation is based on a top level class:
>> FirstOrderIntegratorWithJacobians and five interfaces:
>> StepHandlerWithJacobians, StepInterpolatorWithJacobians,
>> EventHandlerWithJacobians, ParameterizedODE and ODEWithJacobians. The
>> user typically has to implement StepHandlerWithJacobians if he want to
>> monitor the integration process at each step, EventHandlerWithJacobians
>> if he needs to introduce some discrete event in the otherwise continuous
>> integration process and either ParameterizedODE or ODEWithJacobians to
>> describe the problem he wants to integrate. The
>> FirstOrderIntegratorWithJacobians performs all wrapping by itself and
>> mimics many of the original FirstOrderIntegrator interface.
>> Another choice would be to remove the FirstOrderIntegratorWithJacobians
>> class and change the adaptors classes to be user visible
>> (StepHandlerAdaptor, EventHandlerAdaptor, ODEAdaptor). In this case, the
>> user would need to wrap his implementations of the jacobians enabled
>> interface into the adaptors and pass them to a classical integrator. On
>> the one hand, this choice would remove duplication of the integrator
>> methods in a separate class, and would allow better management of
>> handlers by users. Users may also not use adaptors for step handlers or
>> event handlers if they don't explicitly need the jacobians, using only
>> the irst elements of the compound state without any data copyin. On the
>> other hand it will force the user to add the adaptors Here are examples
>> of user code in both cases:
>> Current implementation:
>>   FirstOrderIntegrator rawIntegrator = new XyzIntegrator(...);
>>   ODEWithJacobians ode = new MyProblem(...);
>>   FirstOrderIntegratorWithJacobians integrator =
>>     new FirstOrderIntegratorWithJacobians(rawIntegrator, myProblem);
>>   integrator.addEventHandler(new MyEventHandler(...),
>>                              max, epsilon, count);
>>   integrator.addStepHandler(new MyStepHandler(...));
>>   integrator.integrate(t0, y0, dy0dp, t, y, dydy0, dydp);
>> Implementation with public adaptors:
>>   FirstOrderIntegrator integrator = new XyzIntegrator(...);
>>   ODEWithJacobians ode = new MyProblem(...);
>>   integrator.addEventHandler(new EventAdaptor(new MyEventHandler(...)),
>>                              max, epsilon, count);
>>   integrator.addStepHandler(new StepAdaptor(new MyStepHandler(...)));
>>   StateAdaptor sa = new StateAdaptor(y0, dy0dp);
>>   double[] y = sa.getCompoundStateReference();
>>   integrator.integrate(ode, t0, y t, y);
>>   dydy0 = sa.getDyDy0();
>>   dydp  = sa.getDyDp();
>> My current preference is to have visible adaptors. It seems cleaner from
>> a design perspective, and not too cumbersome for users.
>> What do you think ?
> I agree with the current implementation choice among the two
> described.

OK. So I will stick with the current implementation.

> My only question - which I do not have an answer to as
> this is way out of my domain of practical experience - is is there
> an even simpler setup possible?

I think a simpler alternative would have imply an incompatible change in
user interfaces. The problem is that important methods (from a user
point of view), like FirstOrderIntegrator.integrate and the three
methods in EventHandler use t (the integration variable) and y (the
state vector) only. This extended problem implies providing also dy0/dp
to initialize the integration and it provides the additional results
dy(t)/dy0 and dy(t)/dp, which are the main thing we are interested in.
The two choices above are different ways to deal with this in an upward
compatible way. Current implementation simply create a new classes with
their own methods (integrate, eventOccurred ...) that have the
additional parameters. The public adaptors method provides intermediate
adaptors that pack/unpack (y, dy/dy0, dy/dp) to/from a single array.

A simpler implementation would probably mean getting rid of the y array
and setting up some StateVector class holding (t, y) and extending it to
add the jacobians in some StateVectorWithJacobians. It would also
probably force users to do some casts from StateVector to
StateVectorWithJacobians which would not be elegant either. This would
be a really huge change IMHO.

The examples above show a complete use case when you want the result at
integration end but also the intermediate results (in step handlers) and
have some discrete events to cope with. There are simpler use cases when
only the final state is useful, so only the integrate method is called.

> In any case, examples in the user
> guide and / or class javadoc would be helpful.

Of course! I only did not write them yet because I wanted the API to be
stable. If we go with the current one, I can start writing something.


> Phil
>> Luc
>> ---------------------------------------------------------------------
>> To unsubscribe, e-mail:
>> For additional commands, e-mail:
> ---------------------------------------------------------------------
> To unsubscribe, e-mail:
> For additional commands, e-mail:

To unsubscribe, e-mail:
For additional commands, e-mail:

View raw message