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