/* * Licensed to the Apache Software Foundation (ASF) under one or more * contributor license agreements. See the NOTICE file distributed with * this work for additional information regarding copyright ownership. * The ASF licenses this file to You under the Apache License, Version 2.0 * (the "License"); you may not use this file except in compliance with * the License. You may obtain a copy of the License at * * http://www.apache.org/licenses/LICENSE-2.0 * * Unless required by applicable law or agreed to in writing, software * distributed under the License is distributed on an "AS IS" BASIS, * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. * See the License for the specific language governing permissions and * limitations under the License. */ package org.apache.commons.math.ode.jacobians; import java.io.IOException; import java.io.ObjectInput; import java.io.ObjectOutput; import java.lang.reflect.Array; import java.util.ArrayList; import java.util.Collection; import org.apache.commons.math.MathRuntimeException; import org.apache.commons.math.MaxEvaluationsExceededException; import org.apache.commons.math.ode.DerivativeException; import org.apache.commons.math.exception.util.LocalizedFormats; import org.apache.commons.math.ode.ExtendedFirstOrderDifferentialEquations; import org.apache.commons.math.ode.FirstOrderIntegrator; import org.apache.commons.math.ode.IntegratorException; import org.apache.commons.math.ode.events.EventException; import org.apache.commons.math.ode.events.EventHandler; import org.apache.commons.math.ode.sampling.StepHandler; import org.apache.commons.math.ode.sampling.StepInterpolator; /** This class enhances a first order integrator for differential equations to * compute also partial derivatives of the solution with respect to initial state * and parameters. *
In order to compute both the state and its derivatives, the ODE problem * is extended with jacobians of the raw ODE and the variational equations are * added to form a new compound problem of higher dimension. If the original ODE * problem has dimension n and there are p parameters, the compound problem will * have dimension n × (1 + n + p).
* @see ParameterizedODE * @see ODEWithJacobians * @version $Revision: 1073158 $ $Date: 2011-02-21 22:46:52 +0100 (lun. 21 févr. 2011) $ * @since 2.1 * @deprecated as of 2.2 the complete package is deprecated, it will be replaced * in 3.0 by a completely rewritten implementation */ @Deprecated public class FirstOrderIntegratorWithJacobians { /** Underlying integrator for compound problem. */ private final FirstOrderIntegrator integrator; /** Raw equations to integrate. */ private final ODEWithJacobians ode; /** Maximal number of evaluations allowed. */ private int maxEvaluations; /** Number of evaluations already performed. */ private int evaluations; /** Build an enhanced integrator using internal differentiation to compute jacobians. * @param integrator underlying integrator to solve the compound problem * @param ode original problem (f in the equation y' = f(t, y)) * @param p parameters array (may be null if {@link * ParameterizedODE#getParametersDimension() * getParametersDimension()} from original problem is zero) * @param hY step sizes to use for computing the jacobian df/dy, must have the * same dimension as the original problem * @param hP step sizes to use for computing the jacobian df/dp, must have the * same dimension as the original problem parameters dimension * @see #FirstOrderIntegratorWithJacobians(FirstOrderIntegrator, * ODEWithJacobians) */ public FirstOrderIntegratorWithJacobians(final FirstOrderIntegrator integrator, final ParameterizedODE ode, final double[] p, final double[] hY, final double[] hP) { checkDimension(ode.getDimension(), hY); checkDimension(ode.getParametersDimension(), p); checkDimension(ode.getParametersDimension(), hP); this.integrator = integrator; this.ode = new FiniteDifferencesWrapper(ode, p, hY, hP); setMaxEvaluations(-1); } /** Build an enhanced integrator using ODE builtin jacobian computation features. * @param integrator underlying integrator to solve the compound problem * @param ode original problem, which can compute the jacobians by itself * @see #FirstOrderIntegratorWithJacobians(FirstOrderIntegrator, * ParameterizedODE, double[], double[], double[]) */ public FirstOrderIntegratorWithJacobians(final FirstOrderIntegrator integrator, final ODEWithJacobians ode) { this.integrator = integrator; this.ode = ode; setMaxEvaluations(-1); } /** Add a step handler to this integrator. *The handler will be called by the integrator for each accepted * step.
* @param handler handler for the accepted steps * @see #getStepHandlers() * @see #clearStepHandlers() */ public void addStepHandler(StepHandlerWithJacobians handler) { final int n = ode.getDimension(); final int k = ode.getParametersDimension(); integrator.addStepHandler(new StepHandlerWrapper(handler, n, k)); } /** Get all the step handlers that have been added to the integrator. * @return an unmodifiable collection of the added events handlers * @see #addStepHandler(StepHandlerWithJacobians) * @see #clearStepHandlers() */ public CollectionThis method solves an Initial Value Problem (IVP) and also computes the derivatives * of the solution with respect to initial state and parameters. This can be used as * a basis to solve Boundary Value Problems (BVP).
*Since this method stores some internal state variables made * available in its public interface during integration ({@link * #getCurrentSignedStepsize()}), it is not thread-safe.
* @param t0 initial time * @param y0 initial value of the state vector at t0 * @param dY0dP initial value of the state vector derivative with respect to the * parameters at t0 * @param t target time for the integration * (can be set to a value smaller thant0
for backward integration)
* @param y placeholder where to put the state vector at each successful
* step (and hence at the end of integration), can be the same object as y0
* @param dYdY0 placeholder where to put the state vector derivative with respect
* to the initial state (dy[i]/dy0[j] is in element array dYdY0[i][j]) at each successful
* step (and hence at the end of integration)
* @param dYdP placeholder where to put the state vector derivative with respect
* to the parameters (dy[i]/dp[j] is in element array dYdP[i][j]) at each successful
* step (and hence at the end of integration)
* @return stop time, will be the same as target time if integration reached its
* target, but may be different if some event handler stops it at some point.
* @throws IntegratorException if the integrator cannot perform integration
* @throws DerivativeException this exception is propagated to the caller if
* the underlying user function triggers one
*/
public double integrate(final double t0, final double[] y0, final double[][] dY0dP,
final double t, final double[] y,
final double[][] dYdY0, final double[][] dYdP)
throws DerivativeException, IntegratorException {
final int n = ode.getDimension();
final int k = ode.getParametersDimension();
checkDimension(n, y0);
checkDimension(n, y);
checkDimension(n, dYdY0);
checkDimension(n, dYdY0[0]);
if (k != 0) {
checkDimension(n, dY0dP);
checkDimension(k, dY0dP[0]);
checkDimension(n, dYdP);
checkDimension(k, dYdP[0]);
}
// set up initial state, including partial derivatives
// the compound state z contains the raw state y and its derivatives
// with respect to initial state y0 and to parameters p
// y[i] is stored in z[i]
// dy[i]/dy0[j] is stored in z[n + i * n + j]
// dy[i]/dp[j] is stored in z[n * (n + 1) + i * k + j]
final double[] z = new double[n * (1 + n + k)];
System.arraycopy(y0, 0, z, 0, n);
for (int i = 0; i < n; ++i) {
// set diagonal element of dy/dy0 to 1.0 at t = t0
z[i * (1 + n) + n] = 1.0;
// set initial derivatives with respect to parameters
System.arraycopy(dY0dP[i], 0, z, n * (n + 1) + i * k, k);
}
// integrate the compound state variational equations
evaluations = 0;
final double stopTime = integrator.integrate(new MappingWrapper(), t0, z, t, z);
// dispatch the final compound state into the state and partial derivatives arrays
dispatchCompoundState(z, y, dYdY0, dYdP);
return stopTime;
}
/** Dispatch a compound state array into state and jacobians arrays.
* @param z compound state
* @param y raw state array to fill
* @param dydy0 jacobian array to fill
* @param dydp jacobian array to fill
*/
private static void dispatchCompoundState(final double[] z, final double[] y,
final double[][] dydy0, final double[][] dydp) {
final int n = y.length;
final int k = dydp[0].length;
// state
System.arraycopy(z, 0, y, 0, n);
// jacobian with respect to initial state
for (int i = 0; i < n; ++i) {
System.arraycopy(z, n * (i + 1), dydy0[i], 0, n);
}
// jacobian with respect to parameters
for (int i = 0; i < n; ++i) {
System.arraycopy(z, n * (n + 1) + i * k, dydp[i], 0, k);
}
}
/** Get the current value of the step start time ti.
* This method can be called during integration (typically by * the object implementing the {@link org.apache.commons.math.ode.FirstOrderDifferentialEquations * differential equations} problem) if the value of the current step that * is attempted is needed.
*The result is undefined if the method is called outside of
* calls to integrate
.
This method can be called during integration (typically by * the object implementing the {@link org.apache.commons.math.ode.FirstOrderDifferentialEquations * differential equations} problem) if the signed value of the current stepsize * that is tried is needed.
*The result is undefined if the method is called outside of
* calls to integrate
.
The purpose of this method is to avoid infinite loops which can occur * for example when stringent error constraints are set or when lots of * discrete events are triggered, thus leading to many rejected steps.
* @param maxEvaluations maximal number of function evaluations (negative * values are silently converted to maximal integer value, thus representing * almost unlimited evaluations) */ public void setMaxEvaluations(int maxEvaluations) { this.maxEvaluations = (maxEvaluations < 0) ? Integer.MAX_VALUE : maxEvaluations; } /** Get the maximal number of functions evaluations. * @return maximal number of functions evaluations */ public int getMaxEvaluations() { return maxEvaluations; } /** Get the number of evaluations of the differential equations function. *
* The number of evaluations corresponds to the last call to the
* integrate
method. It is 0 if the method has not been called yet.
*
This constructor is used only for externalization. It does nothing.
*/ @SuppressWarnings("unused") public StepInterpolatorWrapper() { } /** Simple constructor. * @param interpolator wrapped interpolator * @param n dimension of the original ODE * @param k number of parameters */ public StepInterpolatorWrapper(final StepInterpolator interpolator, final int n, final int k) { this.interpolator = interpolator; y = new double[n]; dydy0 = new double[n][n]; dydp = new double[n][k]; yDot = new double[n]; dydy0Dot = new double[n][n]; dydpDot = new double[n][k]; } /** {@inheritDoc} */ public void setInterpolatedTime(double time) { interpolator.setInterpolatedTime(time); } /** {@inheritDoc} */ public boolean isForward() { return interpolator.isForward(); } /** {@inheritDoc} */ public double getPreviousTime() { return interpolator.getPreviousTime(); } /** {@inheritDoc} */ public double getInterpolatedTime() { return interpolator.getInterpolatedTime(); } /** {@inheritDoc} */ public double[] getInterpolatedY() throws DerivativeException { double[] extendedState = interpolator.getInterpolatedState(); System.arraycopy(extendedState, 0, y, 0, y.length); return y; } /** {@inheritDoc} */ public double[][] getInterpolatedDyDy0() throws DerivativeException { double[] extendedState = interpolator.getInterpolatedState(); final int n = y.length; int start = n; for (int i = 0; i < n; ++i) { System.arraycopy(extendedState, start, dydy0[i], 0, n); start += n; } return dydy0; } /** {@inheritDoc} */ public double[][] getInterpolatedDyDp() throws DerivativeException { double[] extendedState = interpolator.getInterpolatedState(); final int n = y.length; final int k = dydp[0].length; int start = n * (n + 1); for (int i = 0; i < n; ++i) { System.arraycopy(extendedState, start, dydp[i], 0, k); start += k; } return dydp; } /** {@inheritDoc} */ public double[] getInterpolatedYDot() throws DerivativeException { double[] extendedDerivatives = interpolator.getInterpolatedDerivatives(); System.arraycopy(extendedDerivatives, 0, yDot, 0, yDot.length); return yDot; } /** {@inheritDoc} */ public double[][] getInterpolatedDyDy0Dot() throws DerivativeException { double[] extendedDerivatives = interpolator.getInterpolatedDerivatives(); final int n = y.length; int start = n; for (int i = 0; i < n; ++i) { System.arraycopy(extendedDerivatives, start, dydy0Dot[i], 0, n); start += n; } return dydy0Dot; } /** {@inheritDoc} */ public double[][] getInterpolatedDyDpDot() throws DerivativeException { double[] extendedDerivatives = interpolator.getInterpolatedDerivatives(); final int n = y.length; final int k = dydpDot[0].length; int start = n * (n + 1); for (int i = 0; i < n; ++i) { System.arraycopy(extendedDerivatives, start, dydpDot[i], 0, k); start += k; } return dydpDot; } /** {@inheritDoc} */ public double getCurrentTime() { return interpolator.getCurrentTime(); } /** {@inheritDoc} */ public StepInterpolatorWithJacobians copy() throws DerivativeException { final int n = y.length; final int k = dydp[0].length; StepInterpolatorWrapper copied = new StepInterpolatorWrapper(interpolator.copy(), n, k); copyArray(y, copied.y); copyArray(dydy0, copied.dydy0); copyArray(dydp, copied.dydp); copyArray(yDot, copied.yDot); copyArray(dydy0Dot, copied.dydy0Dot); copyArray(dydpDot, copied.dydpDot); return copied; } /** {@inheritDoc} */ public void writeExternal(ObjectOutput out) throws IOException { out.writeObject(interpolator); out.writeInt(y.length); out.writeInt(dydp[0].length); writeArray(out, y); writeArray(out, dydy0); writeArray(out, dydp); writeArray(out, yDot); writeArray(out, dydy0Dot); writeArray(out, dydpDot); } /** {@inheritDoc} */ public void readExternal(ObjectInput in) throws IOException, ClassNotFoundException { interpolator = (StepInterpolator) in.readObject(); final int n = in.readInt(); final int k = in.readInt(); y = new double[n]; dydy0 = new double[n][n]; dydp = new double[n][k]; yDot = new double[n]; dydy0Dot = new double[n][n]; dydpDot = new double[n][k]; readArray(in, y); readArray(in, dydy0); readArray(in, dydp); readArray(in, yDot); readArray(in, dydy0Dot); readArray(in, dydpDot); } /** Copy an array. * @param src source array * @param dest destination array */ private static void copyArray(final double[] src, final double[] dest) { System.arraycopy(src, 0, dest, 0, src.length); } /** Copy an array. * @param src source array * @param dest destination array */ private static void copyArray(final double[][] src, final double[][] dest) { for (int i = 0; i < src.length; ++i) { copyArray(src[i], dest[i]); } } /** Write an array. * @param out output stream * @param array array to write * @exception IOException if array cannot be read */ private static void writeArray(final ObjectOutput out, final double[] array) throws IOException { for (int i = 0; i < array.length; ++i) { out.writeDouble(array[i]); } } /** Write an array. * @param out output stream * @param array array to write * @exception IOException if array cannot be read */ private static void writeArray(final ObjectOutput out, final double[][] array) throws IOException { for (int i = 0; i < array.length; ++i) { writeArray(out, array[i]); } } /** Read an array. * @param in input stream * @param array array to read * @exception IOException if array cannot be read */ private static void readArray(final ObjectInput in, final double[] array) throws IOException { for (int i = 0; i < array.length; ++i) { array[i] = in.readDouble(); } } /** Read an array. * @param in input stream * @param array array to read * @exception IOException if array cannot be read */ private static void readArray(final ObjectInput in, final double[][] array) throws IOException { for (int i = 0; i < array.length; ++i) { readArray(in, array[i]); } } } /** Wrapper for event handlers. */ private static class EventHandlerWrapper implements EventHandler { /** Underlying event handler with jacobians. */ private final EventHandlerWithJacobians handler; /** State array. */ private double[] y; /** Jacobian with respect to initial state dy/dy0. */ private double[][] dydy0; /** Jacobian with respect to parameters dy/dp. */ private double[][] dydp; /** Simple constructor. * @param handler underlying event handler with jacobians * @param n dimension of the original ODE * @param k number of parameters */ public EventHandlerWrapper(final EventHandlerWithJacobians handler, final int n, final int k) { this.handler = handler; y = new double[n]; dydy0 = new double[n][n]; dydp = new double[n][k]; } /** Get the underlying event handler with jacobians. * @return underlying event handler with jacobians */ public EventHandlerWithJacobians getHandler() { return handler; } /** {@inheritDoc} */ public int eventOccurred(double t, double[] z, boolean increasing) throws EventException { dispatchCompoundState(z, y, dydy0, dydp); return handler.eventOccurred(t, y, dydy0, dydp, increasing); } /** {@inheritDoc} */ public double g(double t, double[] z) throws EventException { dispatchCompoundState(z, y, dydy0, dydp); return handler.g(t, y, dydy0, dydp); } /** {@inheritDoc} */ public void resetState(double t, double[] z) throws EventException { dispatchCompoundState(z, y, dydy0, dydp); handler.resetState(t, y, dydy0, dydp); } } }