| /* |
| * 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. |
| * <p>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).</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. |
| * <p>The handler will be called by the integrator for each accepted |
| * step.</p> |
| * @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 Collection<StepHandlerWithJacobians> getStepHandlers() { |
| final Collection<StepHandlerWithJacobians> handlers = |
| new ArrayList<StepHandlerWithJacobians>(); |
| for (final StepHandler handler : integrator.getStepHandlers()) { |
| if (handler instanceof StepHandlerWrapper) { |
| handlers.add(((StepHandlerWrapper) handler).getHandler()); |
| } |
| } |
| return handlers; |
| } |
| |
| /** Remove all the step handlers that have been added to the integrator. |
| * @see #addStepHandler(StepHandlerWithJacobians) |
| * @see #getStepHandlers() |
| */ |
| public void clearStepHandlers() { |
| integrator.clearStepHandlers(); |
| } |
| |
| /** Add an event handler to the integrator. |
| * @param handler event handler |
| * @param maxCheckInterval maximal time interval between switching |
| * function checks (this interval prevents missing sign changes in |
| * case the integration steps becomes very large) |
| * @param convergence convergence threshold in the event time search |
| * @param maxIterationCount upper limit of the iteration count in |
| * the event time search |
| * @see #getEventHandlers() |
| * @see #clearEventHandlers() |
| */ |
| public void addEventHandler(EventHandlerWithJacobians handler, |
| double maxCheckInterval, |
| double convergence, |
| int maxIterationCount) { |
| final int n = ode.getDimension(); |
| final int k = ode.getParametersDimension(); |
| integrator.addEventHandler(new EventHandlerWrapper(handler, n, k), |
| maxCheckInterval, convergence, maxIterationCount); |
| } |
| |
| /** Get all the event handlers that have been added to the integrator. |
| * @return an unmodifiable collection of the added events handlers |
| * @see #addEventHandler(EventHandlerWithJacobians, double, double, int) |
| * @see #clearEventHandlers() |
| */ |
| public Collection<EventHandlerWithJacobians> getEventHandlers() { |
| final Collection<EventHandlerWithJacobians> handlers = |
| new ArrayList<EventHandlerWithJacobians>(); |
| for (final EventHandler handler : integrator.getEventHandlers()) { |
| if (handler instanceof EventHandlerWrapper) { |
| handlers.add(((EventHandlerWrapper) handler).getHandler()); |
| } |
| } |
| return handlers; |
| } |
| |
| /** Remove all the event handlers that have been added to the integrator. |
| * @see #addEventHandler(EventHandlerWithJacobians, double, double, int) |
| * @see #getEventHandlers() |
| */ |
| public void clearEventHandlers() { |
| integrator.clearEventHandlers(); |
| } |
| |
| /** Integrate the differential equations and the variational equations up to the given time. |
| * <p>This 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).</p> |
| * <p>Since this method stores some internal state variables made |
| * available in its public interface during integration ({@link |
| * #getCurrentSignedStepsize()}), it is <em>not</em> thread-safe.</p> |
| * @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 than <code>t0</code> 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 t<sub>i</sub>. |
| * <p>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.</p> |
| * <p>The result is undefined if the method is called outside of |
| * calls to <code>integrate</code>.</p> |
| * @return current value of the step start time t<sub>i</sub> |
| */ |
| public double getCurrentStepStart() { |
| return integrator.getCurrentStepStart(); |
| } |
| |
| /** Get the current signed value of the integration stepsize. |
| * <p>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.</p> |
| * <p>The result is undefined if the method is called outside of |
| * calls to <code>integrate</code>.</p> |
| * @return current signed value of the stepsize |
| */ |
| public double getCurrentSignedStepsize() { |
| return integrator.getCurrentSignedStepsize(); |
| } |
| |
| /** Set the maximal number of differential equations function evaluations. |
| * <p>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.</p> |
| * @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. |
| * <p> |
| * The number of evaluations corresponds to the last call to the |
| * <code>integrate</code> method. It is 0 if the method has not been called yet. |
| * </p> |
| * @return number of evaluations of the differential equations function |
| */ |
| public int getEvaluations() { |
| return evaluations; |
| } |
| |
| /** Check array dimensions. |
| * @param expected expected dimension |
| * @param array (may be null if expected is 0) |
| * @throws IllegalArgumentException if the array dimension does not match the expected one |
| */ |
| private void checkDimension(final int expected, final Object array) |
| throws IllegalArgumentException { |
| int arrayDimension = (array == null) ? 0 : Array.getLength(array); |
| if (arrayDimension != expected) { |
| throw MathRuntimeException.createIllegalArgumentException( |
| LocalizedFormats.DIMENSIONS_MISMATCH_SIMPLE, arrayDimension, expected); |
| } |
| } |
| |
| /** Wrapper class used to map state and jacobians into compound state. */ |
| private class MappingWrapper implements ExtendedFirstOrderDifferentialEquations { |
| |
| /** Current state. */ |
| private final double[] y; |
| |
| /** Time derivative of the current state. */ |
| private final double[] yDot; |
| |
| /** Derivatives of yDot with respect to state. */ |
| private final double[][] dFdY; |
| |
| /** Derivatives of yDot with respect to parameters. */ |
| private final double[][] dFdP; |
| |
| /** Simple constructor. |
| */ |
| public MappingWrapper() { |
| |
| final int n = ode.getDimension(); |
| final int k = ode.getParametersDimension(); |
| y = new double[n]; |
| yDot = new double[n]; |
| dFdY = new double[n][n]; |
| dFdP = new double[n][k]; |
| |
| } |
| |
| /** {@inheritDoc} */ |
| public int getDimension() { |
| final int n = y.length; |
| final int k = dFdP[0].length; |
| return n * (1 + n + k); |
| } |
| |
| /** {@inheritDoc} */ |
| public int getMainSetDimension() { |
| return ode.getDimension(); |
| } |
| |
| /** {@inheritDoc} */ |
| public void computeDerivatives(final double t, final double[] z, final double[] zDot) |
| throws DerivativeException { |
| |
| final int n = y.length; |
| final int k = dFdP[0].length; |
| |
| // compute raw ODE and its jacobians: dy/dt, d[dy/dt]/dy0 and d[dy/dt]/dp |
| System.arraycopy(z, 0, y, 0, n); |
| if (++evaluations > maxEvaluations) { |
| throw new DerivativeException(new MaxEvaluationsExceededException(maxEvaluations)); |
| } |
| ode.computeDerivatives(t, y, yDot); |
| ode.computeJacobians(t, y, yDot, dFdY, dFdP); |
| |
| // state part of the compound equations |
| System.arraycopy(yDot, 0, zDot, 0, n); |
| |
| // variational equations: from d[dy/dt]/dy0 to d[dy/dy0]/dt |
| for (int i = 0; i < n; ++i) { |
| final double[] dFdYi = dFdY[i]; |
| for (int j = 0; j < n; ++j) { |
| double s = 0; |
| final int startIndex = n + j; |
| int zIndex = startIndex; |
| for (int l = 0; l < n; ++l) { |
| s += dFdYi[l] * z[zIndex]; |
| zIndex += n; |
| } |
| zDot[startIndex + i * n] = s; |
| } |
| } |
| |
| // variational equations: from d[dy/dt]/dy0 and d[dy/dt]/dp to d[dy/dp]/dt |
| for (int i = 0; i < n; ++i) { |
| final double[] dFdYi = dFdY[i]; |
| final double[] dFdPi = dFdP[i]; |
| for (int j = 0; j < k; ++j) { |
| double s = dFdPi[j]; |
| final int startIndex = n * (n + 1) + j; |
| int zIndex = startIndex; |
| for (int l = 0; l < n; ++l) { |
| s += dFdYi[l] * z[zIndex]; |
| zIndex += k; |
| } |
| zDot[startIndex + i * k] = s; |
| } |
| } |
| |
| } |
| |
| } |
| |
| /** Wrapper class to compute jacobians by finite differences for ODE which do not compute them themselves. */ |
| private class FiniteDifferencesWrapper implements ODEWithJacobians { |
| |
| /** Raw ODE without jacobians computation. */ |
| private final ParameterizedODE ode; |
| |
| /** Parameters array (may be null if parameters dimension from original problem is zero) */ |
| private final double[] p; |
| |
| /** Step sizes to use for computing the jacobian df/dy. */ |
| private final double[] hY; |
| |
| /** Step sizes to use for computing the jacobian df/dp. */ |
| private final double[] hP; |
| |
| /** Temporary array for state derivatives used to compute jacobians. */ |
| private final double[] tmpDot; |
| |
| /** Simple constructor. |
| * @param ode original ODE problem, without jacobians computations |
| * @param p parameters array (may be null if parameters dimension from original problem is zero) |
| * @param hY step sizes to use for computing the jacobian df/dy |
| * @param hP step sizes to use for computing the jacobian df/dp |
| */ |
| public FiniteDifferencesWrapper(final ParameterizedODE ode, |
| final double[] p, final double[] hY, final double[] hP) { |
| this.ode = ode; |
| this.p = p.clone(); |
| this.hY = hY.clone(); |
| this.hP = hP.clone(); |
| tmpDot = new double[ode.getDimension()]; |
| } |
| |
| /** {@inheritDoc} */ |
| public int getDimension() { |
| return ode.getDimension(); |
| } |
| |
| /** {@inheritDoc} */ |
| public void computeDerivatives(double t, double[] y, double[] yDot) throws DerivativeException { |
| // this call to computeDerivatives has already been counted, |
| // we must not increment the counter again |
| ode.computeDerivatives(t, y, yDot); |
| } |
| |
| /** {@inheritDoc} */ |
| public int getParametersDimension() { |
| return ode.getParametersDimension(); |
| } |
| |
| /** {@inheritDoc} */ |
| public void computeJacobians(double t, double[] y, double[] yDot, |
| double[][] dFdY, double[][] dFdP) |
| throws DerivativeException { |
| |
| final int n = hY.length; |
| final int k = hP.length; |
| |
| evaluations += n + k; |
| if (evaluations > maxEvaluations) { |
| throw new DerivativeException(new MaxEvaluationsExceededException(maxEvaluations)); |
| } |
| |
| // compute df/dy where f is the ODE and y is the state array |
| for (int j = 0; j < n; ++j) { |
| final double savedYj = y[j]; |
| y[j] += hY[j]; |
| ode.computeDerivatives(t, y, tmpDot); |
| for (int i = 0; i < n; ++i) { |
| dFdY[i][j] = (tmpDot[i] - yDot[i]) / hY[j]; |
| } |
| y[j] = savedYj; |
| } |
| |
| // compute df/dp where f is the ODE and p is the parameters array |
| for (int j = 0; j < k; ++j) { |
| ode.setParameter(j, p[j] + hP[j]); |
| ode.computeDerivatives(t, y, tmpDot); |
| for (int i = 0; i < n; ++i) { |
| dFdP[i][j] = (tmpDot[i] - yDot[i]) / hP[j]; |
| } |
| ode.setParameter(j, p[j]); |
| } |
| |
| } |
| |
| } |
| |
| /** Wrapper for step handlers. */ |
| private static class StepHandlerWrapper implements StepHandler { |
| |
| /** Underlying step handler with jacobians. */ |
| private final StepHandlerWithJacobians handler; |
| |
| /** Dimension of the original ODE. */ |
| private final int n; |
| |
| /** Number of parameters. */ |
| private final int k; |
| |
| /** Simple constructor. |
| * @param handler underlying step handler with jacobians |
| * @param n dimension of the original ODE |
| * @param k number of parameters |
| */ |
| public StepHandlerWrapper(final StepHandlerWithJacobians handler, |
| final int n, final int k) { |
| this.handler = handler; |
| this.n = n; |
| this.k = k; |
| } |
| |
| /** Get the underlying step handler with jacobians. |
| * @return underlying step handler with jacobians |
| */ |
| public StepHandlerWithJacobians getHandler() { |
| return handler; |
| } |
| |
| /** {@inheritDoc} */ |
| public void handleStep(StepInterpolator interpolator, boolean isLast) |
| throws DerivativeException { |
| handler.handleStep(new StepInterpolatorWrapper(interpolator, n, k), isLast); |
| } |
| |
| /** {@inheritDoc} */ |
| public boolean requiresDenseOutput() { |
| return handler.requiresDenseOutput(); |
| } |
| |
| /** {@inheritDoc} */ |
| public void reset() { |
| handler.reset(); |
| } |
| |
| } |
| |
| /** Wrapper for step interpolators. */ |
| private static class StepInterpolatorWrapper |
| implements StepInterpolatorWithJacobians { |
| |
| /** Wrapped interpolator. */ |
| private StepInterpolator interpolator; |
| |
| /** 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; |
| |
| /** Time derivative of the state array. */ |
| private double[] yDot; |
| |
| /** Time derivative of the sacobian with respect to initial state dy/dy0. */ |
| private double[][] dydy0Dot; |
| |
| /** Time derivative of the jacobian with respect to parameters dy/dp. */ |
| private double[][] dydpDot; |
| |
| /** Simple constructor. |
| * <p>This constructor is used only for externalization. It does nothing.</p> |
| */ |
| @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); |
| } |
| |
| } |
| |
| } |
