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