ParameterJacobianWrapper.java

  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. import java.util.Arrays;
  19. import java.util.Collection;
  20. import java.util.HashMap;
  21. import java.util.Map;

  22. import org.apache.commons.math4.legacy.exception.DimensionMismatchException;
  23. import org.apache.commons.math4.legacy.exception.MaxCountExceededException;

  24. /** Wrapper class to compute Jacobian matrices by finite differences for ODE
  25.  *  which do not compute them by themselves.
  26.  *
  27.  * @since 3.0
  28.  */
  29. class ParameterJacobianWrapper implements ParameterJacobianProvider {

  30.     /** Main ODE set. */
  31.     private final FirstOrderDifferentialEquations fode;

  32.     /** Raw ODE without Jacobian computation skill to be wrapped into a ParameterJacobianProvider. */
  33.     private final ParameterizedODE pode;

  34.     /** Steps for finite difference computation of the Jacobian df/dp w.r.t. parameters. */
  35.     private final Map<String, Double> hParam;

  36.     /** Wrap a {@link ParameterizedODE} into a {@link ParameterJacobianProvider}.
  37.      * @param fode main first order differential equations set
  38.      * @param pode secondary problem, without parameter Jacobian computation skill
  39.      * @param paramsAndSteps parameters and steps to compute the Jacobians df/dp
  40.      * @see JacobianMatrices#setParameterStep(String, double)
  41.      */
  42.     ParameterJacobianWrapper(final FirstOrderDifferentialEquations fode,
  43.                              final ParameterizedODE pode,
  44.                              final ParameterConfiguration[] paramsAndSteps) {
  45.         this.fode = fode;
  46.         this.pode = pode;
  47.         this.hParam = new HashMap<>();

  48.         // set up parameters for jacobian computation
  49.         for (final ParameterConfiguration param : paramsAndSteps) {
  50.             final String name = param.getParameterName();
  51.             if (pode.isSupported(name)) {
  52.                 hParam.put(name, param.getHP());
  53.             }
  54.         }
  55.     }

  56.     /** {@inheritDoc} */
  57.     @Override
  58.     public Collection<String> getParametersNames() {
  59.         return pode.getParametersNames();
  60.     }

  61.     /** {@inheritDoc} */
  62.     @Override
  63.     public boolean isSupported(String name) {
  64.         return pode.isSupported(name);
  65.     }

  66.     /** {@inheritDoc} */
  67.     @Override
  68.     public void computeParameterJacobian(double t, double[] y, double[] yDot,
  69.                                          String paramName, double[] dFdP)
  70.         throws DimensionMismatchException, MaxCountExceededException {

  71.         final int n = fode.getDimension();
  72.         if (pode.isSupported(paramName)) {
  73.             final double[] tmpDot = new double[n];

  74.             // compute the jacobian df/dp w.r.t. parameter
  75.             final double p  = pode.getParameter(paramName);
  76.             final double hP = hParam.get(paramName);
  77.             pode.setParameter(paramName, p + hP);
  78.             fode.computeDerivatives(t, y, tmpDot);
  79.             for (int i = 0; i < n; ++i) {
  80.                 dFdP[i] = (tmpDot[i] - yDot[i]) / hP;
  81.             }
  82.             pode.setParameter(paramName, p);
  83.         } else {
  84.             Arrays.fill(dFdP, 0, n, 0.0);
  85.         }
  86.     }
  87. }