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.math4.legacy.ode;
18  
19  import java.util.Arrays;
20  import java.util.Collection;
21  import java.util.HashMap;
22  import java.util.Map;
23  
24  import org.apache.commons.math4.legacy.exception.DimensionMismatchException;
25  import org.apache.commons.math4.legacy.exception.MaxCountExceededException;
26  
27  /** Wrapper class to compute Jacobian matrices by finite differences for ODE
28   *  which do not compute them by themselves.
29   *
30   * @since 3.0
31   */
32  class ParameterJacobianWrapper implements ParameterJacobianProvider {
33  
34      /** Main ODE set. */
35      private final FirstOrderDifferentialEquations fode;
36  
37      /** Raw ODE without Jacobian computation skill to be wrapped into a ParameterJacobianProvider. */
38      private final ParameterizedODE pode;
39  
40      /** Steps for finite difference computation of the Jacobian df/dp w.r.t. parameters. */
41      private final Map<String, Double> hParam;
42  
43      /** Wrap a {@link ParameterizedODE} into a {@link ParameterJacobianProvider}.
44       * @param fode main first order differential equations set
45       * @param pode secondary problem, without parameter Jacobian computation skill
46       * @param paramsAndSteps parameters and steps to compute the Jacobians df/dp
47       * @see JacobianMatrices#setParameterStep(String, double)
48       */
49      ParameterJacobianWrapper(final FirstOrderDifferentialEquations fode,
50                               final ParameterizedODE pode,
51                               final ParameterConfiguration[] paramsAndSteps) {
52          this.fode = fode;
53          this.pode = pode;
54          this.hParam = new HashMap<>();
55  
56          // set up parameters for jacobian computation
57          for (final ParameterConfiguration param : paramsAndSteps) {
58              final String name = param.getParameterName();
59              if (pode.isSupported(name)) {
60                  hParam.put(name, param.getHP());
61              }
62          }
63      }
64  
65      /** {@inheritDoc} */
66      @Override
67      public Collection<String> getParametersNames() {
68          return pode.getParametersNames();
69      }
70  
71      /** {@inheritDoc} */
72      @Override
73      public boolean isSupported(String name) {
74          return pode.isSupported(name);
75      }
76  
77      /** {@inheritDoc} */
78      @Override
79      public void computeParameterJacobian(double t, double[] y, double[] yDot,
80                                           String paramName, double[] dFdP)
81          throws DimensionMismatchException, MaxCountExceededException {
82  
83          final int n = fode.getDimension();
84          if (pode.isSupported(paramName)) {
85              final double[] tmpDot = new double[n];
86  
87              // compute the jacobian df/dp w.r.t. parameter
88              final double p  = pode.getParameter(paramName);
89              final double hP = hParam.get(paramName);
90              pode.setParameter(paramName, p + hP);
91              fode.computeDerivatives(t, y, tmpDot);
92              for (int i = 0; i < n; ++i) {
93                  dFdP[i] = (tmpDot[i] - yDot[i]) / hP;
94              }
95              pode.setParameter(paramName, p);
96          } else {
97              Arrays.fill(dFdP, 0, n, 0.0);
98          }
99      }
100 }