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