View Javadoc
1   /*
2    * Licensed to the Apache Software Foundation (ASF) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * The ASF licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *      http://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  
18  /*
19   * This is not the original file distributed by the Apache Software Foundation
20   * It has been modified by the Hipparchus project
21   */
22  package org.hipparchus.migration.ode;
23  
24  import java.lang.reflect.Array;
25  import java.lang.reflect.Constructor;
26  import java.lang.reflect.InvocationTargetException;
27  import java.util.ArrayList;
28  import java.util.Arrays;
29  import java.util.List;
30  
31  import org.hipparchus.exception.LocalizedCoreFormats;
32  import org.hipparchus.exception.MathIllegalArgumentException;
33  import org.hipparchus.exception.MathIllegalStateException;
34  import org.hipparchus.ode.ExpandableODE;
35  import org.hipparchus.ode.LocalizedODEFormats;
36  import org.hipparchus.ode.NamedParameterJacobianProvider;
37  import org.hipparchus.ode.ODEState;
38  import org.hipparchus.ode.OrdinaryDifferentialEquation;
39  import org.hipparchus.ode.ParameterConfiguration;
40  import org.hipparchus.ode.ParametersController;
41  import org.hipparchus.ode.SecondaryODE;
42  
43  /**
44   * This class defines a set of {@link SecondaryODE secondary equations} to
45   * compute the Jacobian matrices with respect to the initial state vector and, if
46   * any, to some parameters of the primary ODE set.
47   * <p>
48   * It is intended to be packed into an {@link ExpandableODE}
49   * in conjunction with a primary set of ODE, which may be:
50   * <ul>
51   * <li>a {@link FirstOrderDifferentialEquations}</li>
52   * <li>a {@link MainStateJacobianProvider}</li>
53   * </ul>
54   * In order to compute Jacobian matrices with respect to some parameters of the
55   * primary ODE set, the following parameter Jacobian providers may be set:
56   * <ul>
57   * <li>a {@link ParametersController}</li>
58   * </ul>
59   * </p>
60   *
61   * @see ExpandableODE
62   * @see FirstOrderDifferentialEquations
63   * @see MainStateJacobianProvider
64   * @see NamedParameterJacobianProvider
65   * @see ParametersController
66   * @deprecated as of 1.0, replaced with {@link org.hipparchus.ode.VariationalEquation}
67   */
68  @Deprecated
69  public class JacobianMatrices {
70  
71      /** Expandable first order differential equation. */
72      private ExpandableODE efode;
73  
74      /** Index of the instance in the expandable set. */
75      private int index;
76  
77      /** FODE with exact primary Jacobian computation skill. */
78      private MainStateJacobianProvider jode;
79  
80      /** FODE without exact parameter Jacobian computation skill. */
81      private ParametersController parametersController;
82  
83      /** Primary state vector dimension. */
84      private int stateDim;
85  
86      /** Selected parameters for parameter Jacobian computation. */
87      private MutableParameterConfiguration[] selectedParameters;
88  
89      /** FODE with exact parameter Jacobian computation skill. */
90      private List<NamedParameterJacobianProvider> jacobianProviders;
91  
92      /** Parameters dimension. */
93      private int paramDim;
94  
95      /** Boolean for selected parameters consistency. */
96      private boolean dirtyParameter;
97  
98      /** State and parameters Jacobian matrices in a row. */
99      private double[] matricesData;
100 
101     /** Simple constructor for a secondary equations set computing Jacobian matrices.
102      * <p>
103      * Parameters must belong to the supported ones given by {@link
104      * org.hipparchus.ode.Parameterizable#getParametersNames()}, so the primary set of differential
105      * equations must be {@link org.hipparchus.ode.Parameterizable}.
106      * </p>
107      * <p>Note that each selection clears the previous selected parameters.</p>
108      *
109      * @param fode the primary first order differential equations set to extend
110      * @param hY step used for finite difference computation with respect to state vector
111      * @param parameters parameters to consider for Jacobian matrices processing
112      * (may be null if parameters Jacobians is not desired)
113      * @exception MathIllegalArgumentException if there is a dimension mismatch between
114      * the steps array {@code hY} and the equation dimension
115      */
116     public JacobianMatrices(final OrdinaryDifferentialEquation fode, final double[] hY,
117                             final String... parameters)
118         throws MathIllegalArgumentException {
119         this(new MainStateJacobianWrapper(fode, hY), parameters);
120     }
121 
122     /** Simple constructor for a secondary equations set computing Jacobian matrices.
123      * <p>
124      * Parameters must belong to the supported ones given by {@link
125      * org.hipparchus.ode.Parameterizable#getParametersNames()}, so the primary set of differential
126      * equations must be {@link org.hipparchus.ode.Parameterizable}.
127      * </p>
128      * <p>Note that each selection clears the previous selected parameters.</p>
129      *
130      * @param jode the primary first order differential equations set to extend
131      * @param parameters parameters to consider for Jacobian matrices processing
132      * (may be null if parameters Jacobians is not desired)
133      */
134     public JacobianMatrices(final MainStateJacobianProvider jode,
135                             final String... parameters) {
136 
137         this.efode = null;
138         this.index = -1;
139 
140         this.jode = jode;
141         this.parametersController = null;
142 
143         this.stateDim = jode.getDimension();
144 
145         if (parameters == null) {
146             selectedParameters = null;
147             paramDim = 0;
148         } else {
149             this.selectedParameters = new MutableParameterConfiguration[parameters.length];
150             for (int i = 0; i < parameters.length; ++i) {
151                 selectedParameters[i] = new MutableParameterConfiguration(parameters[i], Double.NaN);
152             }
153             paramDim = parameters.length;
154         }
155         this.dirtyParameter = false;
156 
157         this.jacobianProviders = new ArrayList<>();
158 
159         // set the default initial state Jacobian to the identity
160         // and the default initial parameters Jacobian to the null matrix
161         matricesData = new double[(stateDim + paramDim) * stateDim];
162         for (int i = 0; i < stateDim; ++i) {
163             matricesData[i * (stateDim + 1)] = 1.0;
164         }
165 
166     }
167 
168     /** Register the variational equations for the Jacobians matrices to the expandable set.
169      * <p>
170      * This method must be called <em>before {@link #setUpInitialState(ODEState)}</em>
171      * </p>
172      * @param expandable expandable set into which variational equations should be registered
173      * @throws MathIllegalArgumentException if the dimension of the partial state does not
174      * match the selected equations set dimension
175      * @exception MismatchedEquations if the primary set of the expandable set does
176      * not match the one used to build the instance
177      * @see ExpandableODE#addSecondaryEquations(SecondaryODE)
178      * @see #setUpInitialState(ODEState)
179      */
180     public void registerVariationalEquations(final ExpandableODE expandable)
181         throws MathIllegalArgumentException, MismatchedEquations {
182 
183         // safety checks
184         final OrdinaryDifferentialEquation ode = (jode instanceof MainStateJacobianWrapper) ?
185                                                  ((MainStateJacobianWrapper) jode).ode :
186                                                  jode;
187         if (expandable.getPrimary() != ode) {
188             throw new MismatchedEquations();
189         }
190 
191         efode = expandable;
192         index = efode.addSecondaryEquations(new JacobiansSecondaryODE());
193 
194     }
195 
196     /** Set up initial state.
197      * <p>
198      * This method inserts the initial Jacobian matrices data into
199      * an {@link ODEState ODE state} by overriding the additional
200      * state components corresponding to the instance. It must be
201      * called prior to integrate the equations.
202      * </p>
203      * <p>
204      * This method must be called <em>after</em> {@link
205      * #registerVariationalEquations(ExpandableODE)},
206      * {@link #setInitialMainStateJacobian(double[][])} and
207      * {@link #setInitialParameterJacobian(String, double[])}.
208      * </p>
209      * @param initialState initial state, without the initial Jacobians
210      * matrices
211      * @return a new instance of initial state, with the initial Jacobians
212      * matrices properly initialized
213      */
214     public ODEState setUpInitialState(final ODEState initialState) { // NOPMD - PMD false positive
215 
216         // insert the matrices data into secondary states
217         final double[][] secondary = new double[efode.getMapper().getNumberOfEquations() - 1][];
218         for (int i = 0; i < initialState.getNumberOfSecondaryStates(); ++i) {
219             if (i + 1 != index) {
220                 secondary[i] = initialState.getSecondaryState(i + 1);
221             }
222         }
223         secondary[index - 1] = matricesData;
224 
225         // create an updated initial state
226         return new ODEState(initialState.getTime(), initialState.getPrimaryState(), secondary);
227 
228     }
229 
230     /** Add a parameter Jacobian provider.
231      * @param provider the parameter Jacobian provider to compute exactly the parameter Jacobian matrix
232      */
233     public void addParameterJacobianProvider(final NamedParameterJacobianProvider provider) {
234         jacobianProviders.add(provider);
235     }
236 
237     /** Set a parameter Jacobian provider.
238      * @param pc the controller to compute the parameter Jacobian matrix using finite differences
239      * @deprecated as of 1.0, replaced with {@link #setParametersController(ParametersController)}
240      */
241     @Deprecated
242     public void setParameterizedODE(final ParametersController pc) {
243         setParametersController(pc);
244     }
245 
246     /** Set a parameter Jacobian provider.
247      * @param parametersController the controller to compute the parameter Jacobian matrix using finite differences
248      */
249     public void setParametersController(final ParametersController parametersController) {
250         this.parametersController = parametersController;
251         dirtyParameter = true;
252     }
253 
254     /** Set the step associated to a parameter in order to compute by finite
255      *  difference the Jacobian matrix.
256      * <p>
257      * Needed if and only if the primary ODE set is a {@link ParametersController}.
258      * </p>
259      * <p>
260      * Given a non zero parameter value pval for the parameter, a reasonable value
261      * for such a step is {@code pval * FastMath.sqrt(Precision.EPSILON)}.
262      * </p>
263      * <p>
264      * A zero value for such a step doesn't enable to compute the parameter Jacobian matrix.
265      * </p>
266      * @param parameter parameter to consider for Jacobian processing
267      * @param hP step for Jacobian finite difference computation w.r.t. the specified parameter
268      * @see ParametersController
269      * @exception MathIllegalArgumentException if the parameter is not supported
270      */
271     public void setParameterStep(final String parameter, final double hP)
272         throws MathIllegalArgumentException {
273 
274         for (MutableParameterConfiguration param: selectedParameters) {
275             if (parameter.equals(param.getParameterName())) {
276                 param.setHP(hP);
277                 dirtyParameter = true;
278                 return;
279             }
280         }
281 
282         throw new MathIllegalArgumentException(LocalizedODEFormats.UNKNOWN_PARAMETER, parameter);
283 
284     }
285 
286     /** Set the initial value of the Jacobian matrix with respect to state.
287      * <p>
288      * If this method is not called, the initial value of the Jacobian
289      * matrix with respect to state is set to identity.
290      * </p>
291      * <p>
292      * This method must be called <em>before {@link #setUpInitialState(ODEState)}</em>
293      * </p>
294      * @param dYdY0 initial Jacobian matrix w.r.t. state
295      * @exception MathIllegalArgumentException if matrix dimensions are incorrect
296      */
297     public void setInitialMainStateJacobian(final double[][] dYdY0)
298         throws MathIllegalArgumentException {
299 
300         // Check dimensions
301         checkDimension(stateDim, dYdY0);
302         checkDimension(stateDim, dYdY0[0]);
303 
304         // store the matrix in row major order as a single dimension array
305         int i = 0;
306         for (final double[] row : dYdY0) {
307             System.arraycopy(row, 0, matricesData, i, stateDim);
308             i += stateDim;
309         }
310 
311     }
312 
313     /** Set the initial value of a column of the Jacobian matrix with respect to one parameter.
314      * <p>
315      * If this method is not called for some parameter, the initial value of
316      * the column of the Jacobian matrix with respect to this parameter is set to zero.
317      * </p>
318      * <p>
319      * This method must be called <em>before {@link #setUpInitialState(ODEState)}</em>
320      * </p>
321      * @param pName parameter name
322      * @param dYdP initial Jacobian column vector with respect to the parameter
323      * @exception MathIllegalArgumentException if a parameter is not supported
324      * @throws MathIllegalArgumentException if the column vector does not match state dimension
325      */
326     public void setInitialParameterJacobian(final String pName, final double[] dYdP)
327         throws MathIllegalArgumentException {
328 
329         // Check dimensions
330         checkDimension(stateDim, dYdP);
331 
332         // store the column in a global single dimension array
333         int i = stateDim * stateDim;
334         for (MutableParameterConfiguration param: selectedParameters) {
335             if (pName.equals(param.getParameterName())) {
336                 System.arraycopy(dYdP, 0, matricesData, i, stateDim);
337                 return;
338             }
339             i += stateDim;
340         }
341 
342         throw new MathIllegalArgumentException(LocalizedODEFormats.UNKNOWN_PARAMETER, pName);
343 
344     }
345 
346     /** Extract the Jacobian matrix with respect to state.
347      * @param state state from which to extract Jacobian matrix
348      * @return Jacobian matrix dY/dY0 with respect to state.
349      */
350     public double[][] extractMainSetJacobian(final ODEState state) {
351 
352         // get current state for this set of equations from the expandable fode
353         final double[] p = state.getSecondaryState(index);
354 
355         final double[][] dYdY0 = new double[stateDim][stateDim];
356         int j = 0;
357         for (int i = 0; i < stateDim; i++) {
358             System.arraycopy(p, j, dYdY0[i], 0, stateDim);
359             j += stateDim;
360         }
361 
362         return dYdY0;
363 
364     }
365 
366     /** Extract the Jacobian matrix with respect to one parameter.
367      * @param state state from which to extract Jacobian matrix
368      * @param pName name of the parameter for the computed Jacobian matrix
369      * @return Jacobian matrix dY/dP with respect to the named parameter
370      */
371     public double[] extractParameterJacobian(final ODEState state, final String pName) {
372 
373         // get current state for this set of equations from the expandable fode
374         final double[] p = state.getSecondaryState(index);
375 
376         final double[] dYdP = new double[stateDim];
377         int i = stateDim * stateDim;
378         for (MutableParameterConfiguration param: selectedParameters) {
379             if (param.getParameterName().equals(pName)) {
380                 System.arraycopy(p, i, dYdP, 0, stateDim);
381                 break;
382             }
383             i += stateDim;
384         }
385 
386         return dYdP;
387 
388     }
389 
390     /** Check array dimensions.
391      * @param expected expected dimension
392      * @param array (may be null if expected is 0)
393      * @throws MathIllegalArgumentException if the array dimension does not match the expected one
394      */
395     private void checkDimension(final int expected, final Object array)
396         throws MathIllegalArgumentException {
397         int arrayDimension = (array == null) ? 0 : Array.getLength(array);
398         if (arrayDimension != expected) {
399             throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
400                                                    arrayDimension, expected);
401         }
402     }
403 
404     /** Local implementation of secondary equations.
405      * <p>
406      * This class is an inner class to ensure proper scheduling of calls
407      * by forcing the use of {@link JacobianMatrices#registerVariationalEquations(ExpandableODE)}.
408      * </p>
409      */
410     private class JacobiansSecondaryODE implements SecondaryODE {
411 
412         /** {@inheritDoc} */
413         @Override
414         public int getDimension() {
415             return stateDim * (stateDim + paramDim);
416         }
417 
418         /** {@inheritDoc} */
419         @Override
420         public double[] computeDerivatives(final double t, final double[] y, final double[] yDot,
421                                            final double[] z)
422             throws MathIllegalArgumentException, MathIllegalStateException {
423 
424             try {
425 
426                 // Lazy initialization
427                 Constructor<ParameterConfiguration> configCtr =
428                                 ParameterConfiguration.class.getDeclaredConstructor(String.class, Double.TYPE);
429                 configCtr.setAccessible(true);
430                 @SuppressWarnings("unchecked")
431                 Constructor<NamedParameterJacobianProvider> providerCtr =
432                 (Constructor<NamedParameterJacobianProvider>)
433                 Class.forName("org.hipparchus.ode.ParameterJacobianWrapper").getDeclaredConstructor(OrdinaryDifferentialEquation.class,
434                                                                                                     double[].class,
435                                                                                                     ParametersController.class,
436                                                                                                     ParameterConfiguration[].class);
437                 providerCtr.setAccessible(true);
438                 if (dirtyParameter && (paramDim != 0)) {
439                     ParameterConfiguration [] immutable = new ParameterConfiguration[selectedParameters.length];
440                     for (int i = 0; i < selectedParameters.length; ++i) {
441                         immutable[i] = configCtr.newInstance(selectedParameters[i].getParameterName(),
442                                                              selectedParameters[i].getHP());
443                     }
444                     jacobianProviders.add(providerCtr.newInstance(jode, new double[jode.getDimension()],
445                                                                   parametersController, immutable));
446                     dirtyParameter = false;
447                 }
448 
449             } catch (InstantiationException | IllegalAccessException | IllegalArgumentException |
450                      InvocationTargetException | NoSuchMethodException | SecurityException | ClassNotFoundException e) {
451                 throw new MathIllegalStateException(e, LocalizedCoreFormats.SIMPLE_MESSAGE, e.getLocalizedMessage());
452             }
453 
454             // variational equations:
455             // from d[dy/dt]/dy0 and d[dy/dt]/dp to d[dy/dy0]/dt and d[dy/dp]/dt
456 
457             // compute Jacobian matrix with respect to primary state
458             double[][] dFdY = jode.computeMainStateJacobian(t, y, yDot);
459 
460             // Dispatch Jacobian matrix in the compound secondary state vector
461             final double[] zDot = new double[z.length];
462             for (int i = 0; i < stateDim; ++i) {
463                 final double[] dFdYi = dFdY[i];
464                 for (int j = 0; j < stateDim; ++j) {
465                     double s = 0;
466                     final int startIndex = j;
467                     int zIndex = startIndex;
468                     for (int l = 0; l < stateDim; ++l) {
469                         s += dFdYi[l] * z[zIndex];
470                         zIndex += stateDim;
471                     }
472                     zDot[startIndex + i * stateDim] = s;
473                 }
474             }
475 
476             if (paramDim != 0) {
477                 // compute Jacobian matrices with respect to parameters
478                 int startIndex = stateDim * stateDim;
479                 for (MutableParameterConfiguration param: selectedParameters) {
480                     boolean found = false;
481                     for (int k = 0 ; (!found) && (k < jacobianProviders.size()); ++k) {
482                         final NamedParameterJacobianProvider provider = jacobianProviders.get(k);
483                         if (provider.isSupported(param.getParameterName())) {
484                             final double[] dFdP =
485                                             provider.computeParameterJacobian(t, y, yDot,
486                                                                               param.getParameterName());
487                             for (int i = 0; i < stateDim; ++i) {
488                                 final double[] dFdYi = dFdY[i];
489                                 int zIndex = startIndex;
490                                 double s = dFdP[i];
491                                 for (int l = 0; l < stateDim; ++l) {
492                                     s += dFdYi[l] * z[zIndex];
493                                     zIndex++;
494                                 }
495                                 zDot[startIndex + i] = s;
496                             }
497                             found = true;
498                         }
499                     }
500                     if (! found) {
501                         Arrays.fill(zDot, startIndex, startIndex + stateDim, 0.0);
502                     }
503                     startIndex += stateDim;
504                 }
505             }
506 
507             return zDot;
508 
509         }
510     }
511 
512     /** Wrapper class to compute jacobian matrices by finite differences for ODE
513      *  which do not compute them by themselves.
514      */
515     private static class MainStateJacobianWrapper implements MainStateJacobianProvider {
516 
517         /** Raw ODE without jacobians computation skill to be wrapped into a MainStateJacobianProvider. */
518         private final OrdinaryDifferentialEquation ode;
519 
520         /** Steps for finite difference computation of the jacobian df/dy w.r.t. state. */
521         private final double[] hY;
522 
523         /** Wrap a {@link FirstOrderDifferentialEquations} into a {@link MainStateJacobianProvider}.
524          * @param ode original ODE problem, without jacobians computation skill
525          * @param hY step sizes to compute the jacobian df/dy
526          * @exception MathIllegalArgumentException if there is a dimension mismatch between
527          * the steps array {@code hY} and the equation dimension
528          */
529         MainStateJacobianWrapper(final OrdinaryDifferentialEquation ode,
530                                  final double[] hY)
531             throws MathIllegalArgumentException {
532             this.ode = ode;
533             this.hY = hY.clone();
534             if (hY.length != ode.getDimension()) {
535                 throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
536                                                        ode.getDimension(), hY.length);
537             }
538         }
539 
540         /** {@inheritDoc} */
541         @Override
542         public int getDimension() {
543             return ode.getDimension();
544         }
545 
546         /** {@inheritDoc} */
547         @Override
548         public double[] computeDerivatives(double t, double[] y)
549             throws MathIllegalArgumentException, MathIllegalStateException {
550             return ode.computeDerivatives(t, y);
551         }
552 
553         /** {@inheritDoc} */
554         @Override
555         public double[][] computeMainStateJacobian(double t, double[] y, double[] yDot)
556             throws MathIllegalArgumentException, MathIllegalStateException {
557 
558             final int n = ode.getDimension();
559             final double[][] dFdY = new double[n][n];
560             for (int j = 0; j < n; ++j) {
561                 final double savedYj = y[j];
562                 y[j] += hY[j];
563                 final double[] tmpDot = ode.computeDerivatives(t, y);
564                 for (int i = 0; i < n; ++i) {
565                     dFdY[i][j] = (tmpDot[i] - yDot[i]) / hY[j];
566                 }
567                 y[j] = savedYj;
568             }
569             return dFdY;
570         }
571 
572     }
573 
574     /**
575      * Special exception for equations mismatch.
576      */
577     public static class MismatchedEquations extends MathIllegalArgumentException {
578 
579         /** Serializable UID. */
580         private static final long serialVersionUID = 20120902L;
581 
582         /** Simple constructor. */
583         public MismatchedEquations() {
584             super(LocalizedODEFormats.UNMATCHED_ODE_IN_EXPANDED_SET);
585         }
586 
587     }
588 
589     /** Selected parameter for parameter Jacobian computation. */
590     private static class MutableParameterConfiguration {
591 
592         /** Parameter name. */
593         private String parameterName;
594 
595         /** Parameter step for finite difference computation. */
596         private double hP;
597 
598         /** Parameter name and step pair constructor.
599          * @param parameterName parameter name
600          * @param hP parameter step
601          */
602         MutableParameterConfiguration(final String parameterName, final double hP) {
603             this.parameterName = parameterName;
604             this.hP = hP;
605         }
606 
607         /** Get parameter name.
608          * @return parameterName parameter name
609          */
610         public String getParameterName() {
611             return parameterName;
612         }
613 
614         /** Get parameter step.
615          * @return hP parameter step
616          */
617         public double getHP() {
618             return hP;
619         }
620 
621         /** Set parameter step.
622          * @param hParam parameter step
623          */
624         public void setHP(final double hParam) {
625             this.hP = hParam;
626         }
627 
628     }
629 
630 }
631