src/main/java/org/apache/commons/math/ode/jacobians/FirstOrderIntegratorWithJacobians.java

/*
 * 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 &times; (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);
        }

    }

}