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 * https://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
18 /*
19 * This is not the original file distributed by the Apache Software Foundation
20 * It has been modified by the Hipparchus project
21 */
22 package org.hipparchus.migration.ode;
23
24 import java.lang.reflect.Array;
25 import java.lang.reflect.Constructor;
26 import java.lang.reflect.InvocationTargetException;
27 import java.util.ArrayList;
28 import java.util.Arrays;
29 import java.util.List;
30
31 import org.hipparchus.exception.LocalizedCoreFormats;
32 import org.hipparchus.exception.MathIllegalArgumentException;
33 import org.hipparchus.exception.MathIllegalStateException;
34 import org.hipparchus.ode.ExpandableODE;
35 import org.hipparchus.ode.LocalizedODEFormats;
36 import org.hipparchus.ode.NamedParameterJacobianProvider;
37 import org.hipparchus.ode.ODEState;
38 import org.hipparchus.ode.OrdinaryDifferentialEquation;
39 import org.hipparchus.ode.ParameterConfiguration;
40 import org.hipparchus.ode.ParametersController;
41 import org.hipparchus.ode.SecondaryODE;
42
43 /**
44 * This class defines a set of {@link SecondaryODE secondary equations} to
45 * compute the Jacobian matrices with respect to the initial state vector and, if
46 * any, to some parameters of the primary ODE set.
47 * <p>
48 * It is intended to be packed into an {@link ExpandableODE}
49 * in conjunction with a primary set of ODE, which may be:</p>
50 * <ul>
51 * <li>a {@link FirstOrderDifferentialEquations}</li>
52 * <li>a {@link MainStateJacobianProvider}</li>
53 * </ul>
54 * <p>In order to compute Jacobian matrices with respect to some parameters of the
55 * primary ODE set, the following parameter Jacobian providers may be set:</p>
56 * <ul>
57 * <li>a {@link ParametersController}</li>
58 * </ul>
59 *
60 * @see ExpandableODE
61 * @see FirstOrderDifferentialEquations
62 * @see MainStateJacobianProvider
63 * @see NamedParameterJacobianProvider
64 * @see ParametersController
65 * @deprecated as of 1.0, replaced with {@link org.hipparchus.ode.VariationalEquation}
66 */
67 @Deprecated
68 public class JacobianMatrices {
69
70 /** Expandable first order differential equation. */
71 private ExpandableODE efode;
72
73 /** Index of the instance in the expandable set. */
74 private int index;
75
76 /** FODE with exact primary Jacobian computation skill. */
77 private MainStateJacobianProvider jode;
78
79 /** FODE without exact parameter Jacobian computation skill. */
80 private ParametersController parametersController;
81
82 /** Primary state vector dimension. */
83 private int stateDim;
84
85 /** Selected parameters for parameter Jacobian computation. */
86 private MutableParameterConfiguration[] selectedParameters;
87
88 /** FODE with exact parameter Jacobian computation skill. */
89 private List<NamedParameterJacobianProvider> jacobianProviders;
90
91 /** Parameters dimension. */
92 private int paramDim;
93
94 /** Boolean for selected parameters consistency. */
95 private boolean dirtyParameter;
96
97 /** State and parameters Jacobian matrices in a row. */
98 private double[] matricesData;
99
100 /** Simple constructor for a secondary equations set computing Jacobian matrices.
101 * <p>
102 * Parameters must belong to the supported ones given by {@link
103 * org.hipparchus.ode.Parameterizable#getParametersNames()}, so the primary set of differential
104 * equations must be {@link org.hipparchus.ode.Parameterizable}.
105 * </p>
106 * <p>Note that each selection clears the previous selected parameters.</p>
107 *
108 * @param fode the primary first order differential equations set to extend
109 * @param hY step used for finite difference computation with respect to state vector
110 * @param parameters parameters to consider for Jacobian matrices processing
111 * (may be null if parameters Jacobians is not desired)
112 * @exception MathIllegalArgumentException if there is a dimension mismatch between
113 * the steps array {@code hY} and the equation dimension
114 */
115 public JacobianMatrices(final OrdinaryDifferentialEquation fode, final double[] hY,
116 final String... parameters)
117 throws MathIllegalArgumentException {
118 this(new MainStateJacobianWrapper(fode, hY), parameters);
119 }
120
121 /** Simple constructor for a secondary equations set computing Jacobian matrices.
122 * <p>
123 * Parameters must belong to the supported ones given by {@link
124 * org.hipparchus.ode.Parameterizable#getParametersNames()}, so the primary set of differential
125 * equations must be {@link org.hipparchus.ode.Parameterizable}.
126 * </p>
127 * <p>Note that each selection clears the previous selected parameters.</p>
128 *
129 * @param jode the primary first order differential equations set to extend
130 * @param parameters parameters to consider for Jacobian matrices processing
131 * (may be null if parameters Jacobians is not desired)
132 */
133 public JacobianMatrices(final MainStateJacobianProvider jode,
134 final String... parameters) {
135
136 this.efode = null;
137 this.index = -1;
138
139 this.jode = jode;
140 this.parametersController = null;
141
142 this.stateDim = jode.getDimension();
143
144 if (parameters == null) {
145 selectedParameters = null;
146 paramDim = 0;
147 } else {
148 this.selectedParameters = new MutableParameterConfiguration[parameters.length];
149 for (int i = 0; i < parameters.length; ++i) {
150 selectedParameters[i] = new MutableParameterConfiguration(parameters[i], Double.NaN);
151 }
152 paramDim = parameters.length;
153 }
154 this.dirtyParameter = false;
155
156 this.jacobianProviders = new ArrayList<>();
157
158 // set the default initial state Jacobian to the identity
159 // and the default initial parameters Jacobian to the null matrix
160 matricesData = new double[(stateDim + paramDim) * stateDim];
161 for (int i = 0; i < stateDim; ++i) {
162 matricesData[i * (stateDim + 1)] = 1.0;
163 }
164
165 }
166
167 /** Register the variational equations for the Jacobians matrices to the expandable set.
168 * <p>
169 * This method must be called <em>before {@link #setUpInitialState(ODEState)}</em>
170 * </p>
171 * @param expandable expandable set into which variational equations should be registered
172 * @throws MathIllegalArgumentException if the dimension of the partial state does not
173 * match the selected equations set dimension
174 * @exception MismatchedEquations if the primary set of the expandable set does
175 * not match the one used to build the instance
176 * @see ExpandableODE#addSecondaryEquations(SecondaryODE)
177 * @see #setUpInitialState(ODEState)
178 */
179 public void registerVariationalEquations(final ExpandableODE expandable)
180 throws MathIllegalArgumentException, MismatchedEquations {
181
182 // safety checks
183 final OrdinaryDifferentialEquation ode = (jode instanceof MainStateJacobianWrapper) ?
184 ((MainStateJacobianWrapper) jode).ode :
185 jode;
186 if (expandable.getPrimary() != ode) {
187 throw new MismatchedEquations();
188 }
189
190 efode = expandable;
191 index = efode.addSecondaryEquations(new JacobiansSecondaryODE());
192
193 }
194
195 /** Set up initial state.
196 * <p>
197 * This method inserts the initial Jacobian matrices data into
198 * an {@link ODEState ODE state} by overriding the additional
199 * state components corresponding to the instance. It must be
200 * called prior to integrate the equations.
201 * </p>
202 * <p>
203 * This method must be called <em>after</em> {@link
204 * #registerVariationalEquations(ExpandableODE)},
205 * {@link #setInitialMainStateJacobian(double[][])} and
206 * {@link #setInitialParameterJacobian(String, double[])}.
207 * </p>
208 * @param initialState initial state, without the initial Jacobians
209 * matrices
210 * @return a new instance of initial state, with the initial Jacobians
211 * matrices properly initialized
212 */
213 public ODEState setUpInitialState(final ODEState initialState) { // NOPMD - PMD false positive
214
215 // insert the matrices data into secondary states
216 final double[][] secondary = new double[efode.getMapper().getNumberOfEquations() - 1][];
217 for (int i = 0; i < initialState.getNumberOfSecondaryStates(); ++i) {
218 if (i + 1 != index) {
219 secondary[i] = initialState.getSecondaryState(i + 1);
220 }
221 }
222 secondary[index - 1] = matricesData;
223
224 // create an updated initial state
225 return new ODEState(initialState.getTime(), initialState.getPrimaryState(), secondary);
226
227 }
228
229 /** Add a parameter Jacobian provider.
230 * @param provider the parameter Jacobian provider to compute exactly the parameter Jacobian matrix
231 */
232 public void addParameterJacobianProvider(final NamedParameterJacobianProvider provider) {
233 jacobianProviders.add(provider);
234 }
235
236 /** Set a parameter Jacobian provider.
237 * @param pc the controller to compute the parameter Jacobian matrix using finite differences
238 * @deprecated as of 1.0, replaced with {@link #setParametersController(ParametersController)}
239 */
240 @Deprecated
241 public void setParameterizedODE(final ParametersController pc) {
242 setParametersController(pc);
243 }
244
245 /** Set a parameter Jacobian provider.
246 * @param parametersController the controller to compute the parameter Jacobian matrix using finite differences
247 */
248 public void setParametersController(final ParametersController parametersController) {
249 this.parametersController = parametersController;
250 dirtyParameter = true;
251 }
252
253 /** Set the step associated to a parameter in order to compute by finite
254 * difference the Jacobian matrix.
255 * <p>
256 * Needed if and only if the primary ODE set is a {@link ParametersController}.
257 * </p>
258 * <p>
259 * Given a non zero parameter value pval for the parameter, a reasonable value
260 * for such a step is {@code pval * FastMath.sqrt(Precision.EPSILON)}.
261 * </p>
262 * <p>
263 * A zero value for such a step doesn't enable to compute the parameter Jacobian matrix.
264 * </p>
265 * @param parameter parameter to consider for Jacobian processing
266 * @param hP step for Jacobian finite difference computation w.r.t. the specified parameter
267 * @see ParametersController
268 * @exception MathIllegalArgumentException if the parameter is not supported
269 */
270 public void setParameterStep(final String parameter, final double hP)
271 throws MathIllegalArgumentException {
272
273 for (MutableParameterConfiguration param: selectedParameters) {
274 if (parameter.equals(param.getParameterName())) {
275 param.setHP(hP);
276 dirtyParameter = true;
277 return;
278 }
279 }
280
281 throw new MathIllegalArgumentException(LocalizedODEFormats.UNKNOWN_PARAMETER, parameter);
282
283 }
284
285 /** Set the initial value of the Jacobian matrix with respect to state.
286 * <p>
287 * If this method is not called, the initial value of the Jacobian
288 * matrix with respect to state is set to identity.
289 * </p>
290 * <p>
291 * This method must be called <em>before {@link #setUpInitialState(ODEState)}</em>
292 * </p>
293 * @param dYdY0 initial Jacobian matrix w.r.t. state
294 * @exception MathIllegalArgumentException if matrix dimensions are incorrect
295 */
296 public void setInitialMainStateJacobian(final double[][] dYdY0)
297 throws MathIllegalArgumentException {
298
299 // Check dimensions
300 checkDimension(stateDim, dYdY0);
301 checkDimension(stateDim, dYdY0[0]);
302
303 // store the matrix in row major order as a single dimension array
304 int i = 0;
305 for (final double[] row : dYdY0) {
306 System.arraycopy(row, 0, matricesData, i, stateDim);
307 i += stateDim;
308 }
309
310 }
311
312 /** Set the initial value of a column of the Jacobian matrix with respect to one parameter.
313 * <p>
314 * If this method is not called for some parameter, the initial value of
315 * the column of the Jacobian matrix with respect to this parameter is set to zero.
316 * </p>
317 * <p>
318 * This method must be called <em>before {@link #setUpInitialState(ODEState)}</em>
319 * </p>
320 * @param pName parameter name
321 * @param dYdP initial Jacobian column vector with respect to the parameter
322 * @exception MathIllegalArgumentException if a parameter is not supported
323 * @throws MathIllegalArgumentException if the column vector does not match state dimension
324 */
325 public void setInitialParameterJacobian(final String pName, final double[] dYdP)
326 throws MathIllegalArgumentException {
327
328 // Check dimensions
329 checkDimension(stateDim, dYdP);
330
331 // store the column in a global single dimension array
332 int i = stateDim * stateDim;
333 for (MutableParameterConfiguration param: selectedParameters) {
334 if (pName.equals(param.getParameterName())) {
335 System.arraycopy(dYdP, 0, matricesData, i, stateDim);
336 return;
337 }
338 i += stateDim;
339 }
340
341 throw new MathIllegalArgumentException(LocalizedODEFormats.UNKNOWN_PARAMETER, pName);
342
343 }
344
345 /** Extract the Jacobian matrix with respect to state.
346 * @param state state from which to extract Jacobian matrix
347 * @return Jacobian matrix dY/dY0 with respect to state.
348 */
349 public double[][] extractMainSetJacobian(final ODEState state) {
350
351 // get current state for this set of equations from the expandable fode
352 final double[] p = state.getSecondaryState(index);
353
354 final double[][] dYdY0 = new double[stateDim][stateDim];
355 int j = 0;
356 for (int i = 0; i < stateDim; i++) {
357 System.arraycopy(p, j, dYdY0[i], 0, stateDim);
358 j += stateDim;
359 }
360
361 return dYdY0;
362
363 }
364
365 /** Extract the Jacobian matrix with respect to one parameter.
366 * @param state state from which to extract Jacobian matrix
367 * @param pName name of the parameter for the computed Jacobian matrix
368 * @return Jacobian matrix dY/dP with respect to the named parameter
369 */
370 public double[] extractParameterJacobian(final ODEState state, final String pName) {
371
372 // get current state for this set of equations from the expandable fode
373 final double[] p = state.getSecondaryState(index);
374
375 final double[] dYdP = new double[stateDim];
376 int i = stateDim * stateDim;
377 for (MutableParameterConfiguration param: selectedParameters) {
378 if (param.getParameterName().equals(pName)) {
379 System.arraycopy(p, i, dYdP, 0, stateDim);
380 break;
381 }
382 i += stateDim;
383 }
384
385 return dYdP;
386
387 }
388
389 /** Check array dimensions.
390 * @param expected expected dimension
391 * @param array (may be null if expected is 0)
392 * @throws MathIllegalArgumentException if the array dimension does not match the expected one
393 */
394 private void checkDimension(final int expected, final Object array)
395 throws MathIllegalArgumentException {
396 int arrayDimension = (array == null) ? 0 : Array.getLength(array);
397 if (arrayDimension != expected) {
398 throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
399 arrayDimension, expected);
400 }
401 }
402
403 /** Local implementation of secondary equations.
404 * <p>
405 * This class is an inner class to ensure proper scheduling of calls
406 * by forcing the use of {@link JacobianMatrices#registerVariationalEquations(ExpandableODE)}.
407 * </p>
408 */
409 private class JacobiansSecondaryODE implements SecondaryODE {
410
411 /** {@inheritDoc} */
412 @Override
413 public int getDimension() {
414 return stateDim * (stateDim + paramDim);
415 }
416
417 /** {@inheritDoc} */
418 @Override
419 public double[] computeDerivatives(final double t, final double[] y, final double[] yDot,
420 final double[] z)
421 throws MathIllegalArgumentException, MathIllegalStateException {
422
423 try {
424
425 // Lazy initialization
426 Constructor<ParameterConfiguration> configCtr =
427 ParameterConfiguration.class.getDeclaredConstructor(String.class, Double.TYPE);
428 configCtr.setAccessible(true); // NOPMD
429 @SuppressWarnings("unchecked")
430 Constructor<NamedParameterJacobianProvider> providerCtr =
431 (Constructor<NamedParameterJacobianProvider>)
432 Class.forName("org.hipparchus.ode.ParameterJacobianWrapper").getDeclaredConstructor(OrdinaryDifferentialEquation.class,
433 double[].class,
434 ParametersController.class,
435 ParameterConfiguration[].class);
436 providerCtr.setAccessible(true); // NOPMD
437 if (dirtyParameter && (paramDim != 0)) {
438 ParameterConfiguration [] immutable = new ParameterConfiguration[selectedParameters.length];
439 for (int i = 0; i < selectedParameters.length; ++i) {
440 immutable[i] = configCtr.newInstance(selectedParameters[i].getParameterName(),
441 selectedParameters[i].getHP());
442 }
443 jacobianProviders.add(providerCtr.newInstance(jode, new double[jode.getDimension()],
444 parametersController, immutable));
445 dirtyParameter = false;
446 }
447
448 } catch (InstantiationException | IllegalAccessException | IllegalArgumentException |
449 InvocationTargetException | NoSuchMethodException | SecurityException | ClassNotFoundException e) {
450 throw new MathIllegalStateException(e, LocalizedCoreFormats.SIMPLE_MESSAGE, e.getLocalizedMessage());
451 }
452
453 // variational equations:
454 // from d[dy/dt]/dy0 and d[dy/dt]/dp to d[dy/dy0]/dt and d[dy/dp]/dt
455
456 // compute Jacobian matrix with respect to primary state
457 double[][] dFdY = jode.computeMainStateJacobian(t, y, yDot);
458
459 // Dispatch Jacobian matrix in the compound secondary state vector
460 final double[] zDot = new double[z.length];
461 for (int i = 0; i < stateDim; ++i) {
462 final double[] dFdYi = dFdY[i];
463 for (int j = 0; j < stateDim; ++j) {
464 double s = 0;
465 final int startIndex = j;
466 int zIndex = startIndex;
467 for (int l = 0; l < stateDim; ++l) {
468 s += dFdYi[l] * z[zIndex];
469 zIndex += stateDim;
470 }
471 zDot[startIndex + i * stateDim] = s;
472 }
473 }
474
475 if (paramDim != 0) {
476 // compute Jacobian matrices with respect to parameters
477 int startIndex = stateDim * stateDim;
478 for (MutableParameterConfiguration param: selectedParameters) {
479 boolean found = false;
480 for (int k = 0 ; (!found) && (k < jacobianProviders.size()); ++k) {
481 final NamedParameterJacobianProvider provider = jacobianProviders.get(k);
482 if (provider.isSupported(param.getParameterName())) {
483 final double[] dFdP =
484 provider.computeParameterJacobian(t, y, yDot,
485 param.getParameterName());
486 for (int i = 0; i < stateDim; ++i) {
487 final double[] dFdYi = dFdY[i];
488 int zIndex = startIndex;
489 double s = dFdP[i];
490 for (int l = 0; l < stateDim; ++l) {
491 s += dFdYi[l] * z[zIndex];
492 zIndex++;
493 }
494 zDot[startIndex + i] = s;
495 }
496 found = true;
497 }
498 }
499 if (! found) {
500 Arrays.fill(zDot, startIndex, startIndex + stateDim, 0.0);
501 }
502 startIndex += stateDim;
503 }
504 }
505
506 return zDot;
507
508 }
509 }
510
511 /** Wrapper class to compute jacobian matrices by finite differences for ODE
512 * which do not compute them by themselves.
513 */
514 private static class MainStateJacobianWrapper implements MainStateJacobianProvider {
515
516 /** Raw ODE without jacobians computation skill to be wrapped into a MainStateJacobianProvider. */
517 private final OrdinaryDifferentialEquation ode;
518
519 /** Steps for finite difference computation of the jacobian df/dy w.r.t. state. */
520 private final double[] hY;
521
522 /** Wrap a {@link FirstOrderDifferentialEquations} into a {@link MainStateJacobianProvider}.
523 * @param ode original ODE problem, without jacobians computation skill
524 * @param hY step sizes to compute the jacobian df/dy
525 * @exception MathIllegalArgumentException if there is a dimension mismatch between
526 * the steps array {@code hY} and the equation dimension
527 */
528 MainStateJacobianWrapper(final OrdinaryDifferentialEquation ode,
529 final double[] hY)
530 throws MathIllegalArgumentException {
531 this.ode = ode;
532 this.hY = hY.clone();
533 if (hY.length != ode.getDimension()) {
534 throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
535 ode.getDimension(), hY.length);
536 }
537 }
538
539 /** {@inheritDoc} */
540 @Override
541 public int getDimension() {
542 return ode.getDimension();
543 }
544
545 /** {@inheritDoc} */
546 @Override
547 public double[] computeDerivatives(double t, double[] y)
548 throws MathIllegalArgumentException, MathIllegalStateException {
549 return ode.computeDerivatives(t, y);
550 }
551
552 /** {@inheritDoc} */
553 @Override
554 public double[][] computeMainStateJacobian(double t, double[] y, double[] yDot)
555 throws MathIllegalArgumentException, MathIllegalStateException {
556
557 final int n = ode.getDimension();
558 final double[][] dFdY = new double[n][n];
559 for (int j = 0; j < n; ++j) {
560 final double savedYj = y[j];
561 y[j] += hY[j];
562 final double[] tmpDot = ode.computeDerivatives(t, y);
563 for (int i = 0; i < n; ++i) {
564 dFdY[i][j] = (tmpDot[i] - yDot[i]) / hY[j];
565 }
566 y[j] = savedYj;
567 }
568 return dFdY;
569 }
570
571 }
572
573 /**
574 * Special exception for equations mismatch.
575 */
576 public static class MismatchedEquations extends MathIllegalArgumentException {
577
578 /** Serializable UID. */
579 private static final long serialVersionUID = 20120902L;
580
581 /** Simple constructor. */
582 public MismatchedEquations() {
583 super(LocalizedODEFormats.UNMATCHED_ODE_IN_EXPANDED_SET);
584 }
585
586 }
587
588 /** Selected parameter for parameter Jacobian computation. */
589 private static class MutableParameterConfiguration {
590
591 /** Parameter name. */
592 private String parameterName;
593
594 /** Parameter step for finite difference computation. */
595 private double hP;
596
597 /** Parameter name and step pair constructor.
598 * @param parameterName parameter name
599 * @param hP parameter step
600 */
601 MutableParameterConfiguration(final String parameterName, final double hP) {
602 this.parameterName = parameterName;
603 this.hP = hP;
604 }
605
606 /** Get parameter name.
607 * @return parameterName parameter name
608 */
609 public String getParameterName() {
610 return parameterName;
611 }
612
613 /** Get parameter step.
614 * @return hP parameter step
615 */
616 public double getHP() {
617 return hP;
618 }
619
620 /** Set parameter step.
621 * @param hParam parameter step
622 */
623 public void setHP(final double hParam) {
624 this.hP = hParam;
625 }
626
627 }
628
629 }
630