commons-user mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From Luc Maisonobe <...@spaceroots.org>
Subject Re: [math] ODE with Jacobian in commons math 3.5
Date Fri, 15 May 2015 13:28:37 GMT
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/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 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, e-mail: user-unsubscribe@commons.apache.org
For additional commands, e-mail: user-help@commons.apache.org


Mime
View raw message