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.math3.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.math3.exception.DimensionMismatchException;
25 import org.apache.commons.math3.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      */
ParameterJacobianWrapper(final FirstOrderDifferentialEquations fode, final ParameterizedODE pode, final ParameterConfiguration[] paramsAndSteps)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<String, Double>();
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} */
getParametersNames()66     public Collection<String> getParametersNames() {
67         return pode.getParametersNames();
68     }
69 
70     /** {@inheritDoc} */
isSupported(String name)71     public boolean isSupported(String name) {
72         return pode.isSupported(name);
73     }
74 
75     /** {@inheritDoc} */
computeParameterJacobian(double t, double[] y, double[] yDot, String paramName, double[] dFdP)76     public void computeParameterJacobian(double t, double[] y, double[] yDot,
77                                          String paramName, double[] dFdP)
78         throws DimensionMismatchException, MaxCountExceededException {
79 
80         final int n = fode.getDimension();
81         if (pode.isSupported(paramName)) {
82             final double[] tmpDot = new double[n];
83 
84             // compute the jacobian df/dp w.r.t. parameter
85             final double p  = pode.getParameter(paramName);
86             final double hP = hParam.get(paramName);
87             pode.setParameter(paramName, p + hP);
88             fode.computeDerivatives(t, y, tmpDot);
89             for (int i = 0; i < n; ++i) {
90                 dFdP[i] = (tmpDot[i] - yDot[i]) / hP;
91             }
92             pode.setParameter(paramName, p);
93         } else {
94             Arrays.fill(dFdP, 0, n, 0.0);
95         }
96 
97     }
98 
99 }
100