Mailing-List: contact user-help@commons.apache.org; run by ezmlm
Precedence: bulk
Reply-To: "Commons Users List" <user@commons.apache.org>
MIME-Version: 1.0
In-Reply-To: <5555F485.5040004@spaceroots.org>
References: 
 <CAFxvmmeMyE9Z-4MvJ9Pm8GYFpoKnxoQpSbQ1bgu91U_omf6hyw@mail.gmail.com>
	<5555F485.5040004@spaceroots.org>
Date: Fri, 15 May 2015 14:45:25 +0000
Message-ID: 
 <CAFxvmmdZqVnPrReUjiQzCWL+mPrWUHQ4fcPP2Mz9O_SA9TqS+Q@mail.gmail.com>
Subject: Re: [math] ODE with Jacobian in commons math 3.5
From: Bernard GODARD <bernard.godard@gmail.com>
To: Commons Users List <user@commons.apache.org>
Content-Type: multipart/alternative; boundary=089e0141a8dcfbf00705161fe463

--089e0141a8dcfbf00705161fe463
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: quoted-printable

Hi Luc,

Thank you for your detailed answer.

JIRA issue MATH-1225 was created.


On Fri, May 15, 2015 at 1:28 PM, Luc Maisonobe <luc@spaceroots.org> wrote:

> Le 14/05/2015 16:40, Bernard GODARD a =C3=A9crit :
> > Dear all,
>
> Hi Bernard,
>
> >
> > The user guide on ordinary differential equations
> > http://commons.apache.org/proper/commons-math/userguide/ode.html
> > is very useful to understand how to use ODEWithJacobians and
> > ParameterizedODE.
> > but not up to date.
>
> You are right, we need to update this. Thanks for pointing the problem.
> Could you open a JIRA issue at
> <https://issues.apache.org/jira/browse/MATH> so we don't forget about it?
>
> >
> >>From the latest documentation, it does not seem clear to me how to use
> > ParameterizedODE and ParameterJacobianProvider in math 3.5.
> >
> > Moreover  ParameterJacobianProvider-> computeParameterJacobian looks li=
ke
> > it can only handle one parameter at a time and is computed independentl=
y
> > from the main dynamic equation. For the use case I am interested in
> > (integration of the orbit of a satellite with complex dynamic model of
> > several hundred parameters) I prefer to compute together the force mode=
l
> > and its partial derivatives with respect to all parameters together
> > because these computations share a lot of intermediate results.
> Otherwise I
> > have to either recompute or use a cache mechanism with a test on input
> > parameter change (which also may be invalidated by the ODE integrator
> > depending on the order of the evaluation calls). This design then seems
> > inefficient to me but there might be a good reason for it. I would
> > appreciate if someone could explain this choice.
> >
> > I would rather have a method like:
> >
> > void computeDerivativesWithJacobian(double t, double[] y,
> > parameterList,double[] yDot,double[][] Jacobian )
> > with parameterList an input specifying the ordered list of parameters f=
or
> > which to compute the Jacobian.
>
> The rationale for splitting the computation in several parts (mainly the
> yDot from the Jacobian) was twofold. First, it allowed design of a
> simpler class when Jacobians are not desired, which could then be
> extended by subclasses that add Jacobians support. The second reason was
> that it allowed the Apache Commons Math library to provide by itself  a
> generic implementation using finite differences for computing the
> Jacobians when the user would still provide only the base differential
> equation (this is a case of what Hairer, N=C3=B8rsett and Wanner called
> internal differentiation).
>
> The rationale for splitting the Jacobians itself with one column at a
> time for one parameter is also twofold. First it frees the user (who
> provides this Jacobian) from worrying about putting the right parameter
> in the right column. The matrix is built automatically with proper
> columns assignment even if for some runs all parameters were selected
> and in other runs only a subset of the known parameters are used for
> differentiation, the other parameters being considered fixed. The second
> reason is that the providers for the various columns (i.e. for the
> various parameters) can be different objects. You can even have one
> different class for each parameter if you want. This allows for example
> to handle ODE for which some parameters are computed analytically and
> some other parameters are computed by finite differences. In such a
> case, user would call JacobianMatrices.addParameterJacobianProvider for
> the parameters for which he knows how to differentiate, and he would
> rather call JacobianMatrices.setParameterStep for the parameters for
> which he doesn't know how to differentiate. The use case for these two
> aspect was also orbit determination: user can decide at run time which
> parameters are estimated and which are fixed, hence increasing or
> decreasing the number of columns in the Jacobian, and some force model
> may be complex enough that no analytical expression is available, hence
> needing finite differences.
>
> The current design can be used in your case and benefit from already
> computed values. This is not clear in the documentation, though and
> should be better explained. In fact, the various calls to the
> computeParameterJacobian that are used to retrieve all the columns one
> at a time are all done in a small loop and without any change to the
> state, i.e. they correspond to the same step. The scheduling is as follow=
s:
>
> for each evaluation (i.e. for one value of t and y[]):
>  - one call to the computeMainStateJacobian method
>    (which belongs to either the MainStateJacobianProvider interface
>     or the MainStateJacobianProvider interface, depending on
>     which JacobianMatrices constructor you use)
>  - for each selected parameter, iterating in the exact same
>    order the parameters were specified in JacobianMatrices constructor:
>      * one call to computeParameterJacobian for the current parameter,
>        using the proper provider (which may be a finite differences
>        provider for some parameters and an analytical provider for
>        other parameters)
>  - end for
> end for
>
> With the calls sequence above, you can for example perform all the
> computation early (say when computeMainStateJacobian is called) and
> store the Jacobians columns in an internal array, and then only copy
> the column back when the computeParameterJacobian method is called in
> the internal loop. Of course, this implies the implementations of
> MainStateJacobianProvider and ParameterJacobianProvider know each other.
> This is even easier if they correspond to the same object, but it is not
> required.
>
> So in your case, if you already know how to compute the full Jacobian by
> yourself, I would suggest to put everything in a top level class that
> would implement:
>    - either MainStateJacobianProvider or MainStateJacobianProvider
>    - AND ParameterJacobianProvider
>
> This could be done as follows (note that the computeParameterJacobian
> is complete here, it can really be restricted to two lines of code):
>
> MyClass
>   implements MainStateJacobianProvider, ParameterJacobianProvider {
>
>    private double [][] jacobian;
>    private Map<String, int> nameToIndex;
>
>    public MyClass() {
>      // set up the nameToIndex map
>    }
>
>    public void computeDerivatives(double t, double[] y, double[] yDot) {
>
>      // compute yDot
>
>      // compute and store Jacobian
>      // (beware to store it in transposed form ...)
>
>    }
>
>    void computeParameterJacobian(double t, double[] y, double[] yDot,
>                                   String paramName, double[] dFdP) {
>      // this method is known to be called after computeDerivatives
>      // we simply pick up the column that was already computed before
>      // and that was stored transposed
>      final double[] precomputed =3D jacobian[nameToIndex.get(name)];
>      System.arraycopy(precomputed, 0, dFdP, 0, precomputed.length);
>    }
>
> }
>
>
> Hope this helps,
> Luc
>
> >
> >
> > Thank you,
> >
> > Kind regards,
> >
> >
> >     Bernard Godard
> >
>
>
> ---------------------------------------------------------------------
> To unsubscribe, e-mail: user-unsubscribe@commons.apache.org
> For additional commands, e-mail: user-help@commons.apache.org
>
>

--089e0141a8dcfbf00705161fe463--