001 /*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements. See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License. You may obtain a copy of the License at
008 *
009 * http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017 package org.apache.commons.math3.ode;
018
019 import java.lang.reflect.Array;
020 import java.util.ArrayList;
021 import java.util.Arrays;
022 import java.util.List;
023
024 import org.apache.commons.math3.exception.DimensionMismatchException;
025 import org.apache.commons.math3.exception.MathIllegalArgumentException;
026 import org.apache.commons.math3.exception.MaxCountExceededException;
027 import org.apache.commons.math3.exception.util.LocalizedFormats;
028
029 /**
030 * This class defines a set of {@link SecondaryEquations secondary equations} to
031 * compute the Jacobian matrices with respect to the initial state vector and, if
032 * any, to some parameters of the primary ODE set.
033 * <p>
034 * It is intended to be packed into an {@link ExpandableStatefulODE}
035 * in conjunction with a primary set of ODE, which may be:
036 * <ul>
037 * <li>a {@link FirstOrderDifferentialEquations}</li>
038 * <li>a {@link MainStateJacobianProvider}</li>
039 * </ul>
040 * In order to compute Jacobian matrices with respect to some parameters of the
041 * primary ODE set, the following parameter Jacobian providers may be set:
042 * <ul>
043 * <li>a {@link ParameterJacobianProvider}</li>
044 * <li>a {@link ParameterizedODE}</li>
045 * </ul>
046 * </p>
047 *
048 * @see ExpandableStatefulODE
049 * @see FirstOrderDifferentialEquations
050 * @see MainStateJacobianProvider
051 * @see ParameterJacobianProvider
052 * @see ParameterizedODE
053 *
054 * @version $Id: JacobianMatrices.java 1422447 2012-12-16 01:38:40Z psteitz $
055 * @since 3.0
056 */
057 public class JacobianMatrices {
058
059 /** Expandable first order differential equation. */
060 private ExpandableStatefulODE efode;
061
062 /** Index of the instance in the expandable set. */
063 private int index;
064
065 /** FODE with exact primary Jacobian computation skill. */
066 private MainStateJacobianProvider jode;
067
068 /** FODE without exact parameter Jacobian computation skill. */
069 private ParameterizedODE pode;
070
071 /** Main state vector dimension. */
072 private int stateDim;
073
074 /** Selected parameters for parameter Jacobian computation. */
075 private ParameterConfiguration[] selectedParameters;
076
077 /** FODE with exact parameter Jacobian computation skill. */
078 private List<ParameterJacobianProvider> jacobianProviders;
079
080 /** Parameters dimension. */
081 private int paramDim;
082
083 /** Boolean for selected parameters consistency. */
084 private boolean dirtyParameter;
085
086 /** State and parameters Jacobian matrices in a row. */
087 private double[] matricesData;
088
089 /** Simple constructor for a secondary equations set computing Jacobian matrices.
090 * <p>
091 * Parameters must belong to the supported ones given by {@link
092 * Parameterizable#getParametersNames()}, so the primary set of differential
093 * equations must be {@link Parameterizable}.
094 * </p>
095 * <p>Note that each selection clears the previous selected parameters.</p>
096 *
097 * @param fode the primary first order differential equations set to extend
098 * @param hY step used for finite difference computation with respect to state vector
099 * @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