Hi Luc,
Thank you for your detailed answer.
JIRA issue MATH1225 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 écrit :
> > Dear all,
>
> Hi Bernard,
>
> >
> > The user guide on ordinary differential equations
> > http://commons.apache.org/proper/commonsmath/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 like
> > it can only handle one parameter at a time and is computed independently
> > 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 model
> > 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 for
> > 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ørsett 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 follows:
>
> 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 = jacobian[nameToIndex.get(name)];
> System.arraycopy(precomputed, 0, dFdP, 0, precomputed.length);
> }
>
> }
>
>
> Hope this helps,
> Luc
>
> >
> >
> > Thank you,
> >
> > Kind regards,
> >
> >
> > Bernard Godard
> >
>
>
> 
> To unsubscribe, email: userunsubscribe@commons.apache.org
> For additional commands, email: userhelp@commons.apache.org
>
>
