JacobianMatrices.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;

  19. import java.lang.reflect.Array;
  20. import java.util.ArrayList;
  21. import java.util.Arrays;
  22. import java.util.List;

  24. import org.apache.commons.math4.legacy.exception.DimensionMismatchException;
  25. import org.apache.commons.math4.legacy.exception.MathIllegalArgumentException;
  26. import org.apache.commons.math4.legacy.exception.MaxCountExceededException;
  27. import org.apache.commons.math4.legacy.exception.util.LocalizedFormats;

  29. /**
  30.  * This class defines a set of {@link SecondaryEquations secondary equations} to
  31.  * compute the Jacobian matrices with respect to the initial state vector and, if
  32.  * any, to some parameters of the primary ODE set.
  33.  * <p>
  34.  * It is intended to be packed into an {@link ExpandableStatefulODE}
  35.  * in conjunction with a primary set of ODE, which may be:
  36.  * <ul>
  37.  * <li>a {@link FirstOrderDifferentialEquations}</li>
  38.  * <li>a {@link MainStateJacobianProvider}</li>
  39.  * </ul>
  40.  * In order to compute Jacobian matrices with respect to some parameters of the
  41.  * primary ODE set, the following parameter Jacobian providers may be set:
  42.  * <ul>
  43.  * <li>a {@link ParameterJacobianProvider}</li>
  44.  * <li>a {@link ParameterizedODE}</li>
  45.  * </ul>
  46.  *
  47.  * @see ExpandableStatefulODE
  48.  * @see FirstOrderDifferentialEquations
  49.  * @see MainStateJacobianProvider
  50.  * @see ParameterJacobianProvider
  51.  * @see ParameterizedODE
  52.  *
  53.  * @since 3.0
  54.  */
  55. public class JacobianMatrices {

  57.     /** Expandable first order differential equation. */
  58.     private ExpandableStatefulODE efode;

  60.     /** Index of the instance in the expandable set. */
  61.     private int index;

  63.     /** FODE with exact primary Jacobian computation skill. */
  64.     private MainStateJacobianProvider jode;

  66.     /** FODE without exact parameter Jacobian computation skill. */
  67.     private ParameterizedODE pode;

  69.     /** Main state vector dimension. */
  70.     private int stateDim;

  72.     /** Selected parameters for parameter Jacobian computation. */
  73.     private ParameterConfiguration[] selectedParameters;

  75.     /** FODE with exact parameter Jacobian computation skill. */
  76.     private List<ParameterJacobianProvider> jacobianProviders;

  78.     /** Parameters dimension. */
  79.     private int paramDim;

  81.     /** Boolean for selected parameters consistency. */
  82.     private boolean dirtyParameter;

  84.     /** State and parameters Jacobian matrices in a row. */
  85.     private double[] matricesData;

  87.     /** Simple constructor for a secondary equations set computing Jacobian matrices.
  88.      * <p>
  89.      * Parameters must belong to the supported ones given by {@link
  90.      * Parameterizable#getParametersNames()}, so the primary set of differential
  91.      * equations must be {@link Parameterizable}.
  92.      * </p>
  93.      * <p>Note that each selection clears the previous selected parameters.</p>
  94.      *
  95.      * @param fode the primary first order differential equations set to extend
  96.      * @param hY step used for finite difference computation with respect to state vector
  97.      * @param parameters parameters to consider for Jacobian matrices processing
  98.      * (may be null if parameters Jacobians is not desired)
  99.      * @exception DimensionMismatchException if there is a dimension mismatch between
  100.      * the steps array {@code hY} and the equation dimension
  101.      */
  102.     public JacobianMatrices(final FirstOrderDifferentialEquations fode, final double[] hY,
  103.                             final String... parameters)
  104.         throws DimensionMismatchException {
  105.         this(new MainStateJacobianWrapper(fode, hY), parameters);
  106.     }

  108.     /** Simple constructor for a secondary equations set computing Jacobian matrices.
  109.      * <p>
  110.      * Parameters must belong to the supported ones given by {@link
  111.      * Parameterizable#getParametersNames()}, so the primary set of differential
  112.      * equations must be {@link Parameterizable}.
  113.      * </p>
  114.      * <p>Note that each selection clears the previous selected parameters.</p>
  115.      *
  116.      * @param jode the primary first order differential equations set to extend
  117.      * @param parameters parameters to consider for Jacobian matrices processing
  118.      * (may be null if parameters Jacobians is not desired)
  119.      */
  120.     public JacobianMatrices(final MainStateJacobianProvider jode,
  121.                             final String... parameters) {

  123.         this.efode = null;
  124.         this.index = -1;

  126.         this.jode = jode;
  127.         this.pode = null;

  129.         this.stateDim = jode.getDimension();

  131.         if (parameters == null) {
  132.             selectedParameters = null;
  133.             paramDim = 0;
  134.         } else {
  135.             this.selectedParameters = new ParameterConfiguration[parameters.length];
  136.             for (int i = 0; i < parameters.length; ++i) {
  137.                 selectedParameters[i] = new ParameterConfiguration(parameters[i], Double.NaN);
  138.             }
  139.             paramDim = parameters.length;
  140.         }
  141.         this.dirtyParameter = false;

  143.         this.jacobianProviders = new ArrayList<>();

  145.         // set the default initial state Jacobian to the identity
  146.         // and the default initial parameters Jacobian to the null matrix
  147.         matricesData = new double[(stateDim + paramDim) * stateDim];
  148.         for (int i = 0; i < stateDim; ++i) {
  149.             matricesData[i * (stateDim + 1)] = 1.0;
  150.         }
  151.     }

  153.     /** Register the variational equations for the Jacobians matrices to the expandable set.
  154.      * @param expandable expandable set into which variational equations should be registered
  155.      * @throws DimensionMismatchException if the dimension of the partial state does not
  156.      * match the selected equations set dimension
  157.      * @exception MismatchedEquations if the primary set of the expandable set does
  158.      * not match the one used to build the instance
  159.      * @see ExpandableStatefulODE#addSecondaryEquations(SecondaryEquations)
  160.      */
  161.     public void registerVariationalEquations(final ExpandableStatefulODE expandable)
  162.         throws DimensionMismatchException, MismatchedEquations {

  164.         // safety checks
  165.         final FirstOrderDifferentialEquations ode = (jode instanceof MainStateJacobianWrapper) ?
  166.                                                     ((MainStateJacobianWrapper) jode).ode :
  167.                                                     jode;
  168.         if (expandable.getPrimary() != ode) {
  169.             throw new MismatchedEquations();
  170.         }

  172.         efode = expandable;
  173.         index = efode.addSecondaryEquations(new JacobiansSecondaryEquations());
  174.         efode.setSecondaryState(index, matricesData);
  175.     }

  177.     /** Add a parameter Jacobian provider.
  178.      * @param provider the parameter Jacobian provider to compute exactly the parameter Jacobian matrix
  179.      */
  180.     public void addParameterJacobianProvider(final ParameterJacobianProvider provider) {
  181.         jacobianProviders.add(provider);
  182.     }

  184.     /** Set a parameter Jacobian provider.
  185.      * @param parameterizedOde the parameterized ODE to compute the parameter Jacobian matrix using finite differences
  186.      */
  187.     public void setParameterizedODE(final ParameterizedODE parameterizedOde) {
  188.         this.pode = parameterizedOde;
  189.         dirtyParameter = true;
  190.     }

  192.     /** Set the step associated to a parameter in order to compute by finite
  193.      *  difference the Jacobian matrix.
  194.      * <p>
  195.      * Needed if and only if the primary ODE set is a {@link ParameterizedODE}.
  196.      * </p>
  197.      * <p>
  198.      * Given a non zero parameter value pval for the parameter, a reasonable value
  199.      * for such a step is {@code pval * JdkMath.sqrt(Precision.EPSILON)}.
  200.      * </p>
  201.      * <p>
  202.      * A zero value for such a step doesn't enable to compute the parameter Jacobian matrix.
  203.      * </p>
  204.      * @param parameter parameter to consider for Jacobian processing
  205.      * @param hP step for Jacobian finite difference computation w.r.t. the specified parameter
  206.      * @see ParameterizedODE
  207.      * @exception UnknownParameterException if the parameter is not supported
  208.      */
  209.     public void setParameterStep(final String parameter, final double hP)
  210.         throws UnknownParameterException {

  212.         for (ParameterConfiguration param: selectedParameters) {
  213.             if (parameter.equals(param.getParameterName())) {
  214.                 param.setHP(hP);
  215.                 dirtyParameter = true;
  216.                 return;
  217.             }
  218.         }

  220.         throw new UnknownParameterException(parameter);
  221.     }

  223.     /** Set the initial value of the Jacobian matrix with respect to state.
  224.      * <p>
  225.      * If this method is not called, the initial value of the Jacobian
  226.      * matrix with respect to state is set to identity.
  227.      * </p>
  228.      * @param dYdY0 initial Jacobian matrix w.r.t. state
  229.      * @exception DimensionMismatchException if matrix dimensions are incorrect
  230.      */
  231.     public void setInitialMainStateJacobian(final double[][] dYdY0)
  232.         throws DimensionMismatchException {

  234.         // Check dimensions
  235.         checkDimension(stateDim, dYdY0);
  236.         checkDimension(stateDim, dYdY0[0]);

  238.         // store the matrix in row major order as a single dimension array
  239.         int i = 0;
  240.         for (final double[] row : dYdY0) {
  241.             System.arraycopy(row, 0, matricesData, i, stateDim);
  242.             i += stateDim;
  243.         }

  245.         if (efode != null) {
  246.             efode.setSecondaryState(index, matricesData);
  247.         }
  248.     }

  250.     /** Set the initial value of a column of the Jacobian matrix with respect to one parameter.
  251.      * <p>
  252.      * If this method is not called for some parameter, the initial value of
  253.      * the column of the Jacobian matrix with respect to this parameter is set to zero.
  254.      * </p>
  255.      * @param pName parameter name
  256.      * @param dYdP initial Jacobian column vector with respect to the parameter
  257.      * @exception UnknownParameterException if a parameter is not supported
  258.      * @throws DimensionMismatchException if the column vector does not match state dimension
  259.      */
  260.     public void setInitialParameterJacobian(final String pName, final double[] dYdP)
  261.         throws UnknownParameterException, DimensionMismatchException {

  263.         // Check dimensions
  264.         checkDimension(stateDim, dYdP);

  266.         // store the column in a global single dimension array
  267.         int i = stateDim * stateDim;
  268.         for (ParameterConfiguration param: selectedParameters) {
  269.             if (pName.equals(param.getParameterName())) {
  270.                 System.arraycopy(dYdP, 0, matricesData, i, stateDim);
  271.                 if (efode != null) {
  272.                     efode.setSecondaryState(index, matricesData);
  273.                 }
  274.                 return;
  275.             }
  276.             i += stateDim;
  277.         }

  279.         throw new UnknownParameterException(pName);
  280.     }

  282.     /** Get the current value of the Jacobian matrix with respect to state.
  283.      * @param dYdY0 current Jacobian matrix with respect to state.
  284.      */
  285.     public void getCurrentMainSetJacobian(final double[][] dYdY0) {

  287.         // get current state for this set of equations from the expandable fode
  288.         double[] p = efode.getSecondaryState(index);

  290.         int j = 0;
  291.         for (int i = 0; i < stateDim; i++) {
  292.             System.arraycopy(p, j, dYdY0[i], 0, stateDim);
  293.             j += stateDim;
  294.         }
  295.     }

  297.     /** Get the current value of the Jacobian matrix with respect to one parameter.
  298.      * @param pName name of the parameter for the computed Jacobian matrix
  299.      * @param dYdP current Jacobian matrix with respect to the named parameter
  300.      */
  301.     public void getCurrentParameterJacobian(String pName, final double[] dYdP) {

  303.         // get current state for this set of equations from the expandable fode
  304.         double[] p = efode.getSecondaryState(index);

  306.         int i = stateDim * stateDim;
  307.         for (ParameterConfiguration param: selectedParameters) {
  308.             if (param.getParameterName().equals(pName)) {
  309.                 System.arraycopy(p, i, dYdP, 0, stateDim);
  310.                 return;
  311.             }
  312.             i += stateDim;
  313.         }
  314.     }

  316.     /** Check array dimensions.
  317.      * @param expected expected dimension
  318.      * @param array (may be null if expected is 0)
  319.      * @throws DimensionMismatchException if the array dimension does not match the expected one
  320.      */
  321.     private void checkDimension(final int expected, final Object array)
  322.         throws DimensionMismatchException {
  323.         int arrayDimension = (array == null) ? 0 : Array.getLength(array);
  324.         if (arrayDimension != expected) {
  325.             throw new DimensionMismatchException(arrayDimension, expected);
  326.         }
  327.     }

  329.     /** Local implementation of secondary equations.
  330.      * <p>
  331.      * This class is an inner class to ensure proper scheduling of calls
  332.      * by forcing the use of {@link JacobianMatrices#registerVariationalEquations(ExpandableStatefulODE)}.
  333.      * </p>
  334.      */
  335.     private final class JacobiansSecondaryEquations implements SecondaryEquations {

  337.         /** {@inheritDoc} */
  338.         @Override
  339.         public int getDimension() {
  340.             return stateDim * (stateDim + paramDim);
  341.         }

  343.         /** {@inheritDoc} */
  344.         @Override
  345.         public void computeDerivatives(final double t, final double[] y, final double[] yDot,
  346.                                        final double[] z, final double[] zDot)
  347.             throws MaxCountExceededException, DimensionMismatchException {

  349.             // Lazy initialization
  350.             if (dirtyParameter && paramDim != 0) {
  351.                 jacobianProviders.add(new ParameterJacobianWrapper(jode, pode, selectedParameters));
  352.                 dirtyParameter = false;
  353.             }

  355.             // variational equations:
  356.             // from d[dy/dt]/dy0 and d[dy/dt]/dp to d[dy/dy0]/dt and d[dy/dp]/dt

  358.             // compute Jacobian matrix with respect to primary state
  359.             double[][] dFdY = new double[stateDim][stateDim];
  360.             jode.computeMainStateJacobian(t, y, yDot, dFdY);

  362.             // Dispatch Jacobian matrix in the compound secondary state vector
  363.             for (int i = 0; i < stateDim; ++i) {
  364.                 final double[] dFdYi = dFdY[i];
  365.                 for (int j = 0; j < stateDim; ++j) {
  366.                     double s = 0;
  367.                     final int startIndex = j;
  368.                     int zIndex = startIndex;
  369.                     for (int l = 0; l < stateDim; ++l) {
  370.                         s += dFdYi[l] * z[zIndex];
  371.                         zIndex += stateDim;
  372.                     }
  373.                     zDot[startIndex + i * stateDim] = s;
  374.                 }
  375.             }

  377.             if (paramDim != 0) {
  378.                 // compute Jacobian matrices with respect to parameters
  379.                 double[] dFdP = new double[stateDim];
  380.                 int startIndex = stateDim * stateDim;
  381.                 for (ParameterConfiguration param: selectedParameters) {
  382.                     boolean found = false;
  383.                     for (int k = 0 ; !found && k < jacobianProviders.size(); ++k) {
  384.                         final ParameterJacobianProvider provider = jacobianProviders.get(k);
  385.                         if (provider.isSupported(param.getParameterName())) {
  386.                             provider.computeParameterJacobian(t, y, yDot,
  387.                                                               param.getParameterName(), dFdP);
  388.                             for (int i = 0; i < stateDim; ++i) {
  389.                                 final double[] dFdYi = dFdY[i];
  390.                                 int zIndex = startIndex;
  391.                                 double s = dFdP[i];
  392.                                 for (int l = 0; l < stateDim; ++l) {
  393.                                     s += dFdYi[l] * z[zIndex];
  394.                                     zIndex++;
  395.                                 }
  396.                                 zDot[startIndex + i] = s;
  397.                             }
  398.                             found = true;
  399.                         }
  400.                     }
  401.                     if (! found) {
  402.                         Arrays.fill(zDot, startIndex, startIndex + stateDim, 0.0);
  403.                     }
  404.                     startIndex += stateDim;
  405.                 }
  406.             }
  407.         }
  408.     }

  410.     /** Wrapper class to compute jacobian matrices by finite differences for ODE
  411.      *  which do not compute them by themselves.
  412.      */
  413.     private static final class MainStateJacobianWrapper implements MainStateJacobianProvider {

  415.         /** Raw ODE without jacobians computation skill to be wrapped into a MainStateJacobianProvider. */
  416.         private final FirstOrderDifferentialEquations ode;

  418.         /** Steps for finite difference computation of the jacobian df/dy w.r.t. state. */
  419.         private final double[] hY;

  421.         /** Wrap a {@link FirstOrderDifferentialEquations} into a {@link MainStateJacobianProvider}.
  422.          * @param ode original ODE problem, without jacobians computation skill
  423.          * @param hY step sizes to compute the jacobian df/dy
  424.          * @exception DimensionMismatchException if there is a dimension mismatch between
  425.          * the steps array {@code hY} and the equation dimension
  426.          */
  427.         MainStateJacobianWrapper(final FirstOrderDifferentialEquations ode,
  428.                                  final double[] hY)
  429.             throws DimensionMismatchException {
  430.             this.ode = ode;
  431.             this.hY = hY.clone();
  432.             if (hY.length != ode.getDimension()) {
  433.                 throw new DimensionMismatchException(ode.getDimension(), hY.length);
  434.             }
  435.         }

  437.         /** {@inheritDoc} */
  438.         @Override
  439.         public int getDimension() {
  440.             return ode.getDimension();
  441.         }

  443.         /** {@inheritDoc} */
  444.         @Override
  445.         public void computeDerivatives(double t, double[] y, double[] yDot)
  446.             throws MaxCountExceededException, DimensionMismatchException {
  447.             ode.computeDerivatives(t, y, yDot);
  448.         }

  450.         /** {@inheritDoc} */
  451.         @Override
  452.         public void computeMainStateJacobian(double t, double[] y, double[] yDot, double[][] dFdY)
  453.             throws MaxCountExceededException, DimensionMismatchException {

  455.             final int n = ode.getDimension();
  456.             final double[] tmpDot = new double[n];

  458.             for (int j = 0; j < n; ++j) {
  459.                 final double savedYj = y[j];
  460.                 y[j] += hY[j];
  461.                 ode.computeDerivatives(t, y, tmpDot);
  462.                 for (int i = 0; i < n; ++i) {
  463.                     dFdY[i][j] = (tmpDot[i] - yDot[i]) / hY[j];
  464.                 }
  465.                 y[j] = savedYj;
  466.             }
  467.         }
  468.     }

  470.     /**
  471.      * Special exception for equations mismatch.
  472.      * @since 3.1
  473.      */
  474.     public static class MismatchedEquations extends MathIllegalArgumentException {

  476.         /** Serializable UID. */
  477.         private static final long serialVersionUID = 20120902L;

  479.         /** Simple constructor. */
  480.         public MismatchedEquations() {
  481.             super(LocalizedFormats.UNMATCHED_ODE_IN_EXPANDED_SET);
  482.         }
  483.     }
  484. }
Created with JaCoCo 0.8.10.202304240956Code Coverage Report for Apache Commons Math 4.0-SNAPSHOT