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  package org.apache.commons.math3.ode;
18  
19  import java.lang.reflect.Array;
20  import java.util.ArrayList;
21  import java.util.Arrays;
22  import java.util.List;
23  
24  import org.apache.commons.math3.exception.DimensionMismatchException;
25  import org.apache.commons.math3.exception.MathIllegalArgumentException;
26  import org.apache.commons.math3.exception.MaxCountExceededException;
27  import org.apache.commons.math3.exception.util.LocalizedFormats;
28  
29  /**
30   * This class defines a set of {@link SecondaryEquations secondary equations} to
31   * compute the Jacobian matrices with respect to the initial state vector and, if
32   * any, to some parameters of the primary ODE set.
33   * <p>
34   * It is intended to be packed into an {@link ExpandableStatefulODE}
35   * in conjunction with a primary set of ODE, which may be:
36   * <ul>
37   * <li>a {@link FirstOrderDifferentialEquations}</li>
38   * <li>a {@link MainStateJacobianProvider}</li>
39   * </ul>
40   * In order to compute Jacobian matrices with respect to some parameters of the
41   * primary ODE set, the following parameter Jacobian providers may be set:
42   * <ul>
43   * <li>a {@link ParameterJacobianProvider}</li>
44   * <li>a {@link ParameterizedODE}</li>
45   * </ul>
46   * </p>
47   *
48   * @see ExpandableStatefulODE
49   * @see FirstOrderDifferentialEquations
50   * @see MainStateJacobianProvider
51   * @see ParameterJacobianProvider
52   * @see ParameterizedODE
53   *
54   * @since 3.0
55   */
56  public class JacobianMatrices {
57  
58      /** Expandable first order differential equation. */
59      private ExpandableStatefulODE efode;
60  
61      /** Index of the instance in the expandable set. */
62      private int index;
63  
64      /** FODE with exact primary Jacobian computation skill. */
65      private MainStateJacobianProvider jode;
66  
67      /** FODE without exact parameter Jacobian computation skill. */
68      private ParameterizedODE pode;
69  
70      /** Main state vector dimension. */
71      private int stateDim;
72  
73      /** Selected parameters for parameter Jacobian computation. */
74      private ParameterConfiguration[] selectedParameters;
75  
76      /** FODE with exact parameter Jacobian computation skill. */
77      private List<ParameterJacobianProvider> jacobianProviders;
78  
79      /** Parameters dimension. */
80      private int paramDim;
81  
82      /** Boolean for selected parameters consistency. */
83      private boolean dirtyParameter;
84  
85      /** State and parameters Jacobian matrices in a row. */
86      private double[] matricesData;
87  
88      /** Simple constructor for a secondary equations set computing Jacobian matrices.
89       * <p>
90       * Parameters must belong to the supported ones given by {@link
91       * Parameterizable#getParametersNames()}, so the primary set of differential
92       * equations must be {@link Parameterizable}.
93       * </p>
94       * <p>Note that each selection clears the previous selected parameters.</p>
95       *
96       * @param fode the primary first order differential equations set to extend
97       * @param hY step used for finite difference computation with respect to state vector
98       * @param parameters parameters to consider for Jacobian matrices processing
99       * (may be null if parameters Jacobians is not desired)
100      * @exception DimensionMismatchException if there is a dimension mismatch between
101      * the steps array {@code hY} and the equation dimension
102      */
103     public JacobianMatrices(final FirstOrderDifferentialEquations fode, final double[] hY,
104                             final String... parameters)
105         throws DimensionMismatchException {
106         this(new MainStateJacobianWrapper(fode, hY), parameters);
107     }
108 
109     /** Simple constructor for a secondary equations set computing Jacobian matrices.
110      * <p>
111      * Parameters must belong to the supported ones given by {@link
112      * Parameterizable#getParametersNames()}, so the primary set of differential
113      * equations must be {@link Parameterizable}.
114      * </p>
115      * <p>Note that each selection clears the previous selected parameters.</p>
116      *
117      * @param jode the primary first order differential equations set to extend
118      * @param parameters parameters to consider for Jacobian matrices processing
119      * (may be null if parameters Jacobians is not desired)
120      */
121     public JacobianMatrices(final MainStateJacobianProvider jode,
122                             final String... parameters) {
123 
124         this.efode = null;
125         this.index = -1;
126 
127         this.jode = jode;
128         this.pode = null;
129 
130         this.stateDim = jode.getDimension();
131 
132         if (parameters == null) {
133             selectedParameters = null;
134             paramDim = 0;
135         } else {
136             this.selectedParameters = new ParameterConfiguration[parameters.length];
137             for (int i = 0; i < parameters.length; ++i) {
138                 selectedParameters[i] = new ParameterConfiguration(parameters[i], Double.NaN);
139             }
140             paramDim = parameters.length;
141         }
142         this.dirtyParameter = false;
143 
144         this.jacobianProviders = new ArrayList<ParameterJacobianProvider>();
145 
146         // set the default initial state Jacobian to the identity
147         // and the default initial parameters Jacobian to the null matrix
148         matricesData = new double[(stateDim + paramDim) * stateDim];
149         for (int i = 0; i < stateDim; ++i) {
150             matricesData[i * (stateDim + 1)] = 1.0;
151         }
152 
153     }
154 
155     /** Register the variational equations for the Jacobians matrices to the expandable set.
156      * @param expandable expandable set into which variational equations should be registered
157      * @throws DimensionMismatchException if the dimension of the partial state does not
158      * match the selected equations set dimension
159      * @exception MismatchedEquations if the primary set of the expandable set does
160      * not match the one used to build the instance
161      * @see ExpandableStatefulODE#addSecondaryEquations(SecondaryEquations)
162      */
163     public void registerVariationalEquations(final ExpandableStatefulODE expandable)
164         throws DimensionMismatchException, MismatchedEquations {
165 
166         // safety checks
167         final FirstOrderDifferentialEquations ode = (jode instanceof MainStateJacobianWrapper) ?
168                                                     ((MainStateJacobianWrapper) jode).ode :
169                                                     jode;
170         if (expandable.getPrimary() != ode) {
171             throw new MismatchedEquations();
172         }
173 
174         efode = expandable;
175         index = efode.addSecondaryEquations(new JacobiansSecondaryEquations());
176         efode.setSecondaryState(index, matricesData);
177 
178     }
179 
180     /** Add a parameter Jacobian provider.
181      * @param provider the parameter Jacobian provider to compute exactly the parameter Jacobian matrix
182      */
183     public void addParameterJacobianProvider(final ParameterJacobianProvider provider) {
184         jacobianProviders.add(provider);
185     }
186 
187     /** Set a parameter Jacobian provider.
188      * @param parameterizedOde the parameterized ODE to compute the parameter Jacobian matrix using finite differences
189      */
190     public void setParameterizedODE(final ParameterizedODE parameterizedOde) {
191         this.pode = parameterizedOde;
192         dirtyParameter = true;
193     }
194 
195     /** Set the step associated to a parameter in order to compute by finite
196      *  difference the Jacobian matrix.
197      * <p>
198      * Needed if and only if the primary ODE set is a {@link ParameterizedODE}.
199      * </p>
200      * <p>
201      * Given a non zero parameter value pval for the parameter, a reasonable value
202      * for such a step is {@code pval * FastMath.sqrt(Precision.EPSILON)}.
203      * </p>
204      * <p>
205      * A zero value for such a step doesn't enable to compute the parameter Jacobian matrix.
206      * </p>
207      * @param parameter parameter to consider for Jacobian processing
208      * @param hP step for Jacobian finite difference computation w.r.t. the specified parameter
209      * @see ParameterizedODE
210      * @exception UnknownParameterException if the parameter is not supported
211      */
212     public void setParameterStep(final String parameter, final double hP)
213         throws UnknownParameterException {
214 
215         for (ParameterConfiguration param: selectedParameters) {
216             if (parameter.equals(param.getParameterName())) {
217                 param.setHP(hP);
218                 dirtyParameter = true;
219                 return;
220             }
221         }
222 
223         throw new UnknownParameterException(parameter);
224 
225     }
226 
227     /** Set the initial value of the Jacobian matrix with respect to state.
228      * <p>
229      * If this method is not called, the initial value of the Jacobian
230      * matrix with respect to state is set to identity.
231      * </p>
232      * @param dYdY0 initial Jacobian matrix w.r.t. state
233      * @exception DimensionMismatchException if matrix dimensions are incorrect
234      */
235     public void setInitialMainStateJacobian(final double[][] dYdY0)
236         throws DimensionMismatchException {
237 
238         // Check dimensions
239         checkDimension(stateDim, dYdY0);
240         checkDimension(stateDim, dYdY0[0]);
241 
242         // store the matrix in row major order as a single dimension array
243         int i = 0;
244         for (final double[] row : dYdY0) {
245             System.arraycopy(row, 0, matricesData, i, stateDim);
246             i += stateDim;
247         }
248 
249         if (efode != null) {
250             efode.setSecondaryState(index, matricesData);
251         }
252 
253     }
254 
255     /** Set the initial value of a column of the Jacobian matrix with respect to one parameter.
256      * <p>
257      * If this method is not called for some parameter, the initial value of
258      * the column of the Jacobian matrix with respect to this parameter is set to zero.
259      * </p>
260      * @param pName parameter name
261      * @param dYdP initial Jacobian column vector with respect to the parameter
262      * @exception UnknownParameterException if a parameter is not supported
263      * @throws DimensionMismatchException if the column vector does not match state dimension
264      */
265     public void setInitialParameterJacobian(final String pName, final double[] dYdP)
266         throws UnknownParameterException, DimensionMismatchException {
267 
268         // Check dimensions
269         checkDimension(stateDim, dYdP);
270 
271         // store the column in a global single dimension array
272         int i = stateDim * stateDim;
273         for (ParameterConfiguration param: selectedParameters) {
274             if (pName.equals(param.getParameterName())) {
275                 System.arraycopy(dYdP, 0, matricesData, i, stateDim);
276                 if (efode != null) {
277                     efode.setSecondaryState(index, matricesData);
278                 }
279                 return;
280             }
281             i += stateDim;
282         }
283 
284         throw new UnknownParameterException(pName);
285 
286     }
287 
288     /** Get the current value of the Jacobian matrix with respect to state.
289      * @param dYdY0 current Jacobian matrix with respect to state.
290      */
291     public void getCurrentMainSetJacobian(final double[][] dYdY0) {
292 
293         // get current state for this set of equations from the expandable fode
294         double[] p = efode.getSecondaryState(index);
295 
296         int j = 0;
297         for (int i = 0; i < stateDim; i++) {
298             System.arraycopy(p, j, dYdY0[i], 0, stateDim);
299             j += stateDim;
300         }
301 
302     }
303 
304     /** Get the current value of the Jacobian matrix with respect to one parameter.
305      * @param pName name of the parameter for the computed Jacobian matrix
306      * @param dYdP current Jacobian matrix with respect to the named parameter
307      */
308     public void getCurrentParameterJacobian(String pName, final double[] dYdP) {
309 
310         // get current state for this set of equations from the expandable fode
311         double[] p = efode.getSecondaryState(index);
312 
313         int i = stateDim * stateDim;
314         for (ParameterConfiguration param: selectedParameters) {
315             if (param.getParameterName().equals(pName)) {
316                 System.arraycopy(p, i, dYdP, 0, stateDim);
317                 return;
318             }
319             i += stateDim;
320         }
321 
322     }
323 
324     /** Check array dimensions.
325      * @param expected expected dimension
326      * @param array (may be null if expected is 0)
327      * @throws DimensionMismatchException if the array dimension does not match the expected one
328      */
329     private void checkDimension(final int expected, final Object array)
330         throws DimensionMismatchException {
331         int arrayDimension = (array == null) ? 0 : Array.getLength(array);
332         if (arrayDimension != expected) {
333             throw new DimensionMismatchException(arrayDimension, expected);
334         }
335     }
336 
337     /** Local implementation of secondary equations.
338      * <p>
339      * This class is an inner class to ensure proper scheduling of calls
340      * by forcing the use of {@link JacobianMatrices#registerVariationalEquations(ExpandableStatefulODE)}.
341      * </p>
342      */
343     private class JacobiansSecondaryEquations implements SecondaryEquations {
344 
345         /** {@inheritDoc} */
346         public int getDimension() {
347             return stateDim * (stateDim + paramDim);
348         }
349 
350         /** {@inheritDoc} */
351         public void computeDerivatives(final double t, final double[] y, final double[] yDot,
352                                        final double[] z, final double[] zDot)
353             throws MaxCountExceededException, DimensionMismatchException {
354 
355             // Lazy initialization
356             if (dirtyParameter && (paramDim != 0)) {
357                 jacobianProviders.add(new ParameterJacobianWrapper(jode, pode, selectedParameters));
358                 dirtyParameter = false;
359             }
360 
361             // variational equations:
362             // from d[dy/dt]/dy0 and d[dy/dt]/dp to d[dy/dy0]/dt and d[dy/dp]/dt
363 
364             // compute Jacobian matrix with respect to primary state
365             double[][] dFdY = new double[stateDim][stateDim];
366             jode.computeMainStateJacobian(t, y, yDot, dFdY);
367 
368             // Dispatch Jacobian matrix in the compound secondary state vector
369             for (int i = 0; i < stateDim; ++i) {
370                 final double[] dFdYi = dFdY[i];
371                 for (int j = 0; j < stateDim; ++j) {
372                     double s = 0;
373                     final int startIndex = j;
374                     int zIndex = startIndex;
375                     for (int l = 0; l < stateDim; ++l) {
376                         s += dFdYi[l] * z[zIndex];
377                         zIndex += stateDim;
378                     }
379                     zDot[startIndex + i * stateDim] = s;
380                 }
381             }
382 
383             if (paramDim != 0) {
384                 // compute Jacobian matrices with respect to parameters
385                 double[] dFdP = new double[stateDim];
386                 int startIndex = stateDim * stateDim;
387                 for (ParameterConfiguration param: selectedParameters) {
388                     boolean found = false;
389                     for (int k = 0 ; (!found) && (k < jacobianProviders.size()); ++k) {
390                         final ParameterJacobianProvider provider = jacobianProviders.get(k);
391                         if (provider.isSupported(param.getParameterName())) {
392                             provider.computeParameterJacobian(t, y, yDot,
393                                                               param.getParameterName(), dFdP);
394                             for (int i = 0; i < stateDim; ++i) {
395                                 final double[] dFdYi = dFdY[i];
396                                 int zIndex = startIndex;
397                                 double s = dFdP[i];
398                                 for (int l = 0; l < stateDim; ++l) {
399                                     s += dFdYi[l] * z[zIndex];
400                                     zIndex++;
401                                 }
402                                 zDot[startIndex + i] = s;
403                             }
404                             found = true;
405                         }
406                     }
407                     if (! found) {
408                         Arrays.fill(zDot, startIndex, startIndex + stateDim, 0.0);
409                     }
410                     startIndex += stateDim;
411                 }
412             }
413 
414         }
415     }
416 
417     /** Wrapper class to compute jacobian matrices by finite differences for ODE
418      *  which do not compute them by themselves.
419      */
420     private static class MainStateJacobianWrapper implements MainStateJacobianProvider {
421 
422         /** Raw ODE without jacobians computation skill to be wrapped into a MainStateJacobianProvider. */
423         private final FirstOrderDifferentialEquations ode;
424 
425         /** Steps for finite difference computation of the jacobian df/dy w.r.t. state. */
426         private final double[] hY;
427 
428         /** Wrap a {@link FirstOrderDifferentialEquations} into a {@link MainStateJacobianProvider}.
429          * @param ode original ODE problem, without jacobians computation skill
430          * @param hY step sizes to compute the jacobian df/dy
431          * @exception DimensionMismatchException if there is a dimension mismatch between
432          * the steps array {@code hY} and the equation dimension
433          */
434         public MainStateJacobianWrapper(final FirstOrderDifferentialEquations ode,
435                                         final double[] hY)
436             throws DimensionMismatchException {
437             this.ode = ode;
438             this.hY = hY.clone();
439             if (hY.length != ode.getDimension()) {
440                 throw new DimensionMismatchException(ode.getDimension(), hY.length);
441             }
442         }
443 
444         /** {@inheritDoc} */
445         public int getDimension() {
446             return ode.getDimension();
447         }
448 
449         /** {@inheritDoc} */
450         public void computeDerivatives(double t, double[] y, double[] yDot)
451             throws MaxCountExceededException, DimensionMismatchException {
452             ode.computeDerivatives(t, y, yDot);
453         }
454 
455         /** {@inheritDoc} */
456         public void computeMainStateJacobian(double t, double[] y, double[] yDot, double[][] dFdY)
457             throws MaxCountExceededException, DimensionMismatchException {
458 
459             final int n = ode.getDimension();
460             final double[] tmpDot = new double[n];
461 
462             for (int j = 0; j < n; ++j) {
463                 final double savedYj = y[j];
464                 y[j] += hY[j];
465                 ode.computeDerivatives(t, y, tmpDot);
466                 for (int i = 0; i < n; ++i) {
467                     dFdY[i][j] = (tmpDot[i] - yDot[i]) / hY[j];
468                 }
469                 y[j] = savedYj;
470             }
471         }
472 
473     }
474 
475     /**
476      * Special exception for equations mismatch.
477      * @since 3.1
478      */
479     public static class MismatchedEquations extends MathIllegalArgumentException {
480 
481         /** Serializable UID. */
482         private static final long serialVersionUID = 20120902L;
483 
484         /** Simple constructor. */
485         public MismatchedEquations() {
486             super(LocalizedFormats.UNMATCHED_ODE_IN_EXPANDED_SET);
487         }
488 
489     }
490 
491 }
492