    001    /*
    002     * Licensed to the Apache Software Foundation (ASF) under one or more
    003     * contributor license agreements.  See the NOTICE file distributed with
    004     * this work for additional information regarding copyright ownership.
    005     * The ASF licenses this file to You under the Apache License, Version 2.0
    006     * (the "License"); you may not use this file except in compliance with
    007     * the License.  You may obtain a copy of the License at
    008     *
    009     *      http://www.apache.org/licenses/LICENSE-2.0
    010     *
    011     * Unless required by applicable law or agreed to in writing, software
    012     * distributed under the License is distributed on an "AS IS" BASIS,
    013     * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
    014     * See the License for the specific language governing permissions and
    015     * limitations under the License.
    016     */
    017
    018    package org.apache.commons.math.ode.jacobians;
    019
    020    import java.io.IOException;
    021    import java.io.ObjectInput;
    022    import java.io.ObjectOutput;
    023    import java.lang.reflect.Array;
    024    import java.util.ArrayList;
    025    import java.util.Collection;
    026
    027    import org.apache.commons.math.MathRuntimeException;
    028    import org.apache.commons.math.MaxEvaluationsExceededException;
    029    import org.apache.commons.math.ode.DerivativeException;
    030    import org.apache.commons.math.exception.util.LocalizedFormats;
    031    import org.apache.commons.math.ode.ExtendedFirstOrderDifferentialEquations;
    032    import org.apache.commons.math.ode.FirstOrderIntegrator;
    033    import org.apache.commons.math.ode.IntegratorException;
    034    import org.apache.commons.math.ode.events.EventException;
    035    import org.apache.commons.math.ode.events.EventHandler;
    036    import org.apache.commons.math.ode.sampling.StepHandler;
    037    import org.apache.commons.math.ode.sampling.StepInterpolator;
    038
    039    /** This class enhances a first order integrator for differential equations to
    040     * compute also partial derivatives of the solution with respect to initial state
    041     * and parameters.
    042     * <p>In order to compute both the state and its derivatives, the ODE problem
    043     * is extended with jacobians of the raw ODE and the variational equations are
    044     * added to form a new compound problem of higher dimension. If the original ODE
    045     * problem has dimension n and there are p parameters, the compound problem will
    046     * have dimension n &times; (1 + n + p).</p>
    047     * @see ParameterizedODE
    048     * @see ODEWithJacobians
    049     * @version $Revision: 1073158 $ $Date: 2011-02-21 22:46:52 +0100 (lun. 21 f??vr. 2011) $
    050     * @since 2.1
    051     * @deprecated as of 2.2 the complete package is deprecated, it will be replaced
    052     * in 3.0 by a completely rewritten implementation
    053     */
    054    @Deprecated
    055    public class FirstOrderIntegratorWithJacobians {
    056
    057        /** Underlying integrator for compound problem. */
    058        private final FirstOrderIntegrator integrator;
    059
    060        /** Raw equations to integrate. */
    061        private final ODEWithJacobians ode;
    062
    063        /** Maximal number of evaluations allowed. */
    064        private int maxEvaluations;
    065
    066        /** Number of evaluations already performed. */
    067        private int evaluations;
    068
    069        /** Build an enhanced integrator using internal differentiation to compute jacobians.
    070         * @param integrator underlying integrator to solve the compound problem
    071         * @param ode original problem (f in the equation y' = f(t, y))
    072         * @param p parameters array (may be null if {@link
    073         * ParameterizedODE#getParametersDimension()
    074         * getParametersDimension()} from original problem is zero)
    075         * @param hY step sizes to use for computing the jacobian df/dy, must have the
    076         * same dimension as the original problem
    077         * @param hP step sizes to use for computing the jacobian df/dp, must have the
    078         * same dimension as the original problem parameters dimension
    079         * @see #FirstOrderIntegratorWithJacobians(FirstOrderIntegrator,
    080         * ODEWithJacobians)
    081         */
    082        public FirstOrderIntegratorWithJacobians(final FirstOrderIntegrator integrator,
    083                                                 final ParameterizedODE ode,
    084                                                 final double[] p, final double[] hY, final double[] hP) {
    085            checkDimension(ode.getDimension(), hY);
    086            checkDimension(ode.getParametersDimension(), p);
    087            checkDimension(ode.getParametersDimension(), hP);
    088            this.integrator = integrator;
    089            this.ode = new FiniteDifferencesWrapper(ode, p, hY, hP);
    090            setMaxEvaluations(-1);
    091        }
    092
    093        /** Build an enhanced integrator using ODE builtin jacobian computation features.
    094         * @param integrator underlying integrator to solve the compound problem
    095         * @param ode original problem, which can compute the jacobians by itself
    096         * @see #FirstOrderIntegratorWithJacobians(FirstOrderIntegrator,
    097         * ParameterizedODE, double[], double[], double[])
    098         */
    099        public FirstOrderIntegratorWithJacobians(final FirstOrderIntegrator integrator,
    100                                                 final ODEWithJacobians ode) {
    101            this.integrator = integrator;
    102            this.ode = ode;
    103            setMaxEvaluations(-1);
    104        }
    105
    106        /** Add a step handler to this integrator.
    107         * <p>The handler will be called by the integrator for each accepted
    108         * step.</p>
    109         * @param handler handler for the accepted steps
    110         * @see #getStepHandlers()
    111         * @see #clearStepHandlers()
    112         */
    113        public void addStepHandler(StepHandlerWithJacobians handler) {
    114            final int n = ode.getDimension();
    115            final int k = ode.getParametersDimension();
    116            integrator.addStepHandler(new StepHandlerWrapper(handler, n, k));
    117        }
    118
    119        /** Get all the step handlers that have been added to the integrator.
    120         * @return an unmodifiable collection of the added events handlers
    121         * @see #addStepHandler(StepHandlerWithJacobians)
    122         * @see #clearStepHandlers()
    123         */
    124        public Collection<StepHandlerWithJacobians> getStepHandlers() {
    125            final Collection<StepHandlerWithJacobians> handlers =
    126                new ArrayList<StepHandlerWithJacobians>();
    127            for (final StepHandler handler : integrator.getStepHandlers()) {
    128                if (handler instanceof StepHandlerWrapper) {
    129                    handlers.add(((StepHandlerWrapper) handler).getHandler());
    130                }
    131            }
    132            return handlers;
    133        }
    134
    135        /** Remove all the step handlers that have been added to the integrator.
    136         * @see #addStepHandler(StepHandlerWithJacobians)
    137         * @see #getStepHandlers()
    138         */
    139        public void clearStepHandlers() {
    140            integrator.clearStepHandlers();
    141        }
    142
    143        /** Add an event handler to the integrator.
    144         * @param handler event handler
    145         * @param maxCheckInterval maximal time interval between switching
    146         * function checks (this interval prevents missing sign changes in
    147         * case the integration steps becomes very large)
    148         * @param convergence convergence threshold in the event time search
    149         * @param maxIterationCount upper limit of the iteration count in
    150         * the event time search
    151         * @see #getEventHandlers()
    152         * @see #clearEventHandlers()
    153         */
    154        public void addEventHandler(EventHandlerWithJacobians handler,
    155                                    double maxCheckInterval,
    156                                    double convergence,
    157                                    int maxIterationCount) {
    158            final int n = ode.getDimension();
    159            final int k = ode.getParametersDimension();
    160            integrator.addEventHandler(new EventHandlerWrapper(handler, n, k),
    161                                       maxCheckInterval, convergence, maxIterationCount);
    162        }
    163
    164        /** Get all the event handlers that have been added to the integrator.
    165         * @return an unmodifiable collection of the added events handlers
    166         * @see #addEventHandler(EventHandlerWithJacobians, double, double, int)
    167         * @see #clearEventHandlers()
    168         */
    169        public Collection<EventHandlerWithJacobians> getEventHandlers() {
    170            final Collection<EventHandlerWithJacobians> handlers =
    171                new ArrayList<EventHandlerWithJacobians>();
    172            for (final EventHandler handler : integrator.getEventHandlers()) {
    173                if (handler instanceof EventHandlerWrapper) {
    174                    handlers.add(((EventHandlerWrapper) handler).getHandler());
    175                }
    176            }
    177            return handlers;
    178        }
    179
    180        /** Remove all the event handlers that have been added to the integrator.
    181         * @see #addEventHandler(EventHandlerWithJacobians, double, double, int)
    182         * @see #getEventHandlers()
    183         */
    184        public void clearEventHandlers() {
    185            integrator.clearEventHandlers();
    186        }
    187
    188        /** Integrate the differential equations and the variational equations up to the given time.
    189         * <p>This method solves an Initial Value Problem (IVP) and also computes the derivatives
    190         * of the solution with respect to initial state and parameters. This can be used as
    191         * a basis to solve Boundary Value Problems (BVP).</p>
    192         * <p>Since this method stores some internal state variables made
    193         * available in its public interface during integration ({@link
    194         * #getCurrentSignedStepsize()}), it is <em>not</em> thread-safe.</p>
    195         * @param t0 initial time
    196         * @param y0 initial value of the state vector at t0
    197         * @param dY0dP initial value of the state vector derivative with respect to the
    198         * parameters at t0
    199         * @param t target time for the integration
    200         * (can be set to a value smaller than <code>t0</code> for backward integration)
    201         * @param y placeholder where to put the state vector at each successful
    202         *  step (and hence at the end of integration), can be the same object as y0
    203         * @param dYdY0 placeholder where to put the state vector derivative with respect
    204         * to the initial state (dy[i]/dy0[j] is in element array dYdY0[i][j]) at each successful
    205         *  step (and hence at the end of integration)
    206         * @param dYdP placeholder where to put the state vector derivative with respect
    207         * to the parameters (dy[i]/dp[j] is in element array dYdP[i][j]) at each successful
    208         *  step (and hence at the end of integration)
    209         * @return stop time, will be the same as target time if integration reached its
    210         * target, but may be different if some event handler stops it at some point.
    211         * @throws IntegratorException if the integrator cannot perform integration
    212         * @throws DerivativeException this exception is propagated to the caller if
    213         * the underlying user function triggers one
    214         */
    215        public double integrate(final double t0, final double[] y0, final double[][] dY0dP,
    216                                final double t, final double[] y,
    217                                final double[][] dYdY0, final double[][] dYdP)
    218            throws DerivativeException, IntegratorException {
    219
    220            final int n = ode.getDimension();
    221            final int k = ode.getParametersDimension();
    222            checkDimension(n, y0);
    223            checkDimension(n, y);
    224            checkDimension(n, dYdY0);
    225            checkDimension(n, dYdY0[0]);
    226            if (k != 0) {
    227                checkDimension(n, dY0dP);
    228                checkDimension(k, dY0dP[0]);
    229                checkDimension(n, dYdP);
    230                checkDimension(k, dYdP[0]);
    231            }
    232
    233            // set up initial state, including partial derivatives
    234            // the compound state z contains the raw state y and its derivatives
    235            // with respect to initial state y0 and to parameters p
    236            //    y[i]         is stored in z[i]
    237            //    dy[i]/dy0[j] is stored in z[n + i * n + j]
    238            //    dy[i]/dp[j]  is stored in z[n * (n + 1) + i * k + j]
    239            final double[] z = new double[n * (1 + n + k)];
    240            System.arraycopy(y0, 0, z, 0, n);
    241            for (int i = 0; i < n; ++i) {
    242
    243                // set diagonal element of dy/dy0 to 1.0 at t = t0
    244                z[i * (1 + n) + n] = 1.0;
    245
    246                // set initial derivatives with respect to parameters
    247                System.arraycopy(dY0dP[i], 0, z, n * (n + 1) + i * k, k);
    248
    249            }
    250
    251            // integrate the compound state variational equations
    252            evaluations = 0;
    253            final double stopTime = integrator.integrate(new MappingWrapper(), t0, z, t, z);
    254
    255            // dispatch the final compound state into the state and partial derivatives arrays
    256            dispatchCompoundState(z, y, dYdY0, dYdP);
    257
    258            return stopTime;
    259
    260        }
    261
    262        /** Dispatch a compound state array into state and jacobians arrays.
    263         * @param z compound state
    264         * @param y raw state array to fill
    265         * @param dydy0 jacobian array to fill
    266         * @param dydp jacobian array to fill
    267         */
    268        private static void dispatchCompoundState(final double[] z, final double[] y,
    269                                                  final double[][] dydy0, final double[][] dydp) {
    270
    271            final int n = y.length;
    272            final int k = dydp[0].length;
    273
    274            // state
    275            System.arraycopy(z, 0, y, 0, n);
    276
    277            // jacobian with respect to initial state
    278            for (int i = 0; i < n; ++i) {
    279                System.arraycopy(z, n * (i + 1), dydy0[i], 0, n);
    280            }
    281
    282            // jacobian with respect to parameters
    283            for (int i = 0; i < n; ++i) {
    284                System.arraycopy(z, n * (n + 1) + i * k, dydp[i], 0, k);
    285            }
    286
    287        }
    288
    289        /** Get the current value of the step start time t<sub>i</sub>.
    290         * <p>This method can be called during integration (typically by
    291         * the object implementing the {@link org.apache.commons.math.ode.FirstOrderDifferentialEquations
    292         * differential equations} problem) if the value of the current step that
    293         * is attempted is needed.</p>
    294         * <p>The result is undefined if the method is called outside of
    295         * calls to <code>integrate</code>.</p>
    296         * @return current value of the step start time t<sub>i</sub>
    297         */
    298        public double getCurrentStepStart() {
    299            return integrator.getCurrentStepStart();
    300        }
    301
    302        /** Get the current signed value of the integration stepsize.
    303         * <p>This method can be called during integration (typically by
    304         * the object implementing the {@link org.apache.commons.math.ode.FirstOrderDifferentialEquations
    305         * differential equations} problem) if the signed value of the current stepsize
    306         * that is tried is needed.</p>
    307         * <p>The result is undefined if the method is called outside of
    308         * calls to <code>integrate</code>.</p>
    309         * @return current signed value of the stepsize
    310         */
    311        public double getCurrentSignedStepsize() {
    312            return integrator.getCurrentSignedStepsize();
    313        }
    314
    315        /** Set the maximal number of differential equations function evaluations.
    316         * <p>The purpose of this method is to avoid infinite loops which can occur
    317         * for example when stringent error constraints are set or when lots of
    318         * discrete events are triggered, thus leading to many rejected steps.</p>
    319         * @param maxEvaluations maximal number of function evaluations (negative
    320         * values are silently converted to maximal integer value, thus representing
    321         * almost unlimited evaluations)
    322         */
    323        public void setMaxEvaluations(int maxEvaluations) {
    324            this.maxEvaluations = (maxEvaluations < 0) ? Integer.MAX_VALUE : maxEvaluations;
    325        }
    326
    327        /** Get the maximal number of functions evaluations.
    328         * @return maximal number of functions evaluations
    329         */
    330        public int getMaxEvaluations() {
    331            return maxEvaluations;
    332        }
    333
    334        /** Get the number of evaluations of the differential equations function.
    335         * <p>
    336         * The number of evaluations corresponds to the last call to the
    337         * <code>integrate</code> method. It is 0 if the method has not been called yet.
    338         * </p>
    339         * @return number of evaluations of the differential equations function
    340         */
    341        public int getEvaluations() {
    342            return evaluations;
    343        }
    344
    345        /** Check array dimensions.
    346         * @param expected expected dimension
    347         * @param array (may be null if expected is 0)
    348         * @throws IllegalArgumentException if the array dimension does not match the expected one
    349         */
    350        private void checkDimension(final int expected, final Object array)
    351            throws IllegalArgumentException {
    352            int arrayDimension = (array == null) ? 0 : Array.getLength(array);
    353            if (arrayDimension != expected) {
    354                throw MathRuntimeException.createIllegalArgumentException(
    355                      LocalizedFormats.DIMENSIONS_MISMATCH_SIMPLE, arrayDimension, expected);
    356            }
    357        }
    358
    359        /** Wrapper class used to map state and jacobians into compound state. */
    360        private class MappingWrapper implements  ExtendedFirstOrderDifferentialEquations {
    361
    362            /** Current state. */
    363            private final double[]   y;
    364
    365            /** Time derivative of the current state. */
    366            private final double[]   yDot;
    367
    368            /** Derivatives of yDot with respect to state. */
    369            private final double[][] dFdY;
    370
    371            /** Derivatives of yDot with respect to parameters. */
    372            private final double[][] dFdP;
    373
    374            /** Simple constructor.
    375             */
    376            public MappingWrapper() {
    377
    378                final int n = ode.getDimension();
    379                final int k = ode.getParametersDimension();
    380                y    = new double[n];
    381                yDot = new double[n];
    382                dFdY = new double[n][n];
    383                dFdP = new double[n][k];
    384
    385            }
    386
    387            /** {@inheritDoc} */
    388            public int getDimension() {
    389                final int n = y.length;
    390                final int k = dFdP[0].length;
    391                return n * (1 + n + k);
    392            }
    393
    394            /** {@inheritDoc} */
    395            public int getMainSetDimension() {
    396                return ode.getDimension();
    397            }
    398
    399            /** {@inheritDoc} */
    400            public void computeDerivatives(final double t, final double[] z, final double[] zDot)
    401                throws DerivativeException {
    402
    403                final int n = y.length;
    404                final int k = dFdP[0].length;
    405
    406                // compute raw ODE and its jacobians: dy/dt, d[dy/dt]/dy0 and d[dy/dt]/dp
    407                System.arraycopy(z,    0, y,    0, n);
    408                if (++evaluations > maxEvaluations) {
    409                    throw new DerivativeException(new MaxEvaluationsExceededException(maxEvaluations));
    410                }
    411                ode.computeDerivatives(t, y, yDot);
    412                ode.computeJacobians(t, y, yDot, dFdY, dFdP);
    413
    414                // state part of the compound equations
    415                System.arraycopy(yDot, 0, zDot, 0, n);
    416
    417                // variational equations: from d[dy/dt]/dy0 to d[dy/dy0]/dt
    418                for (int i = 0; i < n; ++i) {
    419                    final double[] dFdYi = dFdY[i];
    420                    for (int j = 0; j < n; ++j) {
    421                        double s = 0;
    422                        final int startIndex = n + j;
    423                        int zIndex = startIndex;
    424                        for (int l = 0; l < n; ++l) {
    425                            s += dFdYi[l] * z[zIndex];
    426                            zIndex += n;
    427                        }
    428                        zDot[startIndex + i * n] = s;
    429                    }
    430                }
    431
    432                // variational equations: from d[dy/dt]/dy0 and d[dy/dt]/dp to d[dy/dp]/dt
    433                for (int i = 0; i < n; ++i) {
    434                    final double[] dFdYi = dFdY[i];
    435                    final double[] dFdPi = dFdP[i];
    436                    for (int j = 0; j < k; ++j) {
    437                        double s = dFdPi[j];
    438                        final int startIndex = n * (n + 1) + j;
    439                        int zIndex = startIndex;
    440                        for (int l = 0; l < n; ++l) {
    441                            s += dFdYi[l] * z[zIndex];
    442                            zIndex += k;
    443                        }
    444                        zDot[startIndex + i * k] = s;
    445                    }
    446                }
    447
    448            }
    449
    450        }
    451
    452        /** Wrapper class to compute jacobians by finite differences for ODE which do not compute them themselves. */
    453        private class FiniteDifferencesWrapper implements ODEWithJacobians {
    454
    455            /** Raw ODE without jacobians computation. */
    456            private final ParameterizedODE ode;
    457
    458            /** Parameters array (may be null if parameters dimension from original problem is zero) */
    459            private final double[] p;
    460
    461            /** Step sizes to use for computing the jacobian df/dy. */
    462            private final double[] hY;
    463
    464            /** Step sizes to use for computing the jacobian df/dp. */
    465            private final double[] hP;
    466
    467            /** Temporary array for state derivatives used to compute jacobians. */
    468            private final double[] tmpDot;
    469
    470            /** Simple constructor.
    471             * @param ode original ODE problem, without jacobians computations
    472             * @param p parameters array (may be null if parameters dimension from original problem is zero)
    473             * @param hY step sizes to use for computing the jacobian df/dy
    474             * @param hP step sizes to use for computing the jacobian df/dp
    475             */
    476            public FiniteDifferencesWrapper(final ParameterizedODE ode,
    477                                            final double[] p, final double[] hY, final double[] hP) {
    478                this.ode = ode;
    479                this.p  = p.clone();
    480                this.hY = hY.clone();
    481                this.hP = hP.clone();
    482                tmpDot = new double[ode.getDimension()];
    483            }
    484
    485            /** {@inheritDoc} */
    486            public int getDimension() {
    487                return ode.getDimension();
    488            }
    489
    490            /** {@inheritDoc} */
    491            public void computeDerivatives(double t, double[] y, double[] yDot) throws DerivativeException {
    492                // this call to computeDerivatives has already been counted,
    493                // we must not increment the counter again
    494                ode.computeDerivatives(t, y, yDot);
    495            }
    496
    497            /** {@inheritDoc} */
    498            public int getParametersDimension() {
    499                return ode.getParametersDimension();
    500            }
    501
    502            /** {@inheritDoc} */
    503            public void computeJacobians(double t, double[] y, double[] yDot,
    504                                         double[][] dFdY, double[][] dFdP)
    505                throws DerivativeException {
    506
    507                final int n = hY.length;
    508                final int k = hP.length;
    509
    510                evaluations += n + k;
    511                if (evaluations > maxEvaluations) {
    512                    throw new DerivativeException(new MaxEvaluationsExceededException(maxEvaluations));
    513                }
    514
    515                // compute df/dy where f is the ODE and y is the state array
    516                for (int j = 0; j < n; ++j) {
    517                    final double savedYj = y[j];
    518                    y[j] += hY[j];
    519                    ode.computeDerivatives(t, y, tmpDot);
    520                    for (int i = 0; i < n; ++i) {
    521                        dFdY[i][j] = (tmpDot[i] - yDot[i]) / hY[j];
    522                    }
    523                    y[j] = savedYj;
    524                }
    525
    526                // compute df/dp where f is the ODE and p is the parameters array
    527                for (int j = 0; j < k; ++j) {
    528                    ode.setParameter(j, p[j] +  hP[j]);
    529                    ode.computeDerivatives(t, y, tmpDot);
    530                    for (int i = 0; i < n; ++i) {
    531                        dFdP[i][j] = (tmpDot[i] - yDot[i]) / hP[j];
    532                    }
    533                    ode.setParameter(j, p[j]);
    534                }
    535
    536            }
    537
    538        }
    539
    540        /** Wrapper for step handlers. */
    541        private static class StepHandlerWrapper implements StepHandler {
    542
    543            /** Underlying step handler with jacobians. */
    544            private final StepHandlerWithJacobians handler;
    545
    546            /** Dimension of the original ODE. */
    547            private final int n;
    548
    549            /** Number of parameters. */
    550            private final int k;
    551
    552            /** Simple constructor.
    553             * @param handler underlying step handler with jacobians
    554             * @param n dimension of the original ODE
    555             * @param k number of parameters
    556             */
    557            public StepHandlerWrapper(final StepHandlerWithJacobians handler,
    558                                      final int n, final int k) {
    559                this.handler = handler;
    560                this.n       = n;
    561                this.k       = k;
    562            }
    563
    564            /** Get the underlying step handler with jacobians.
    565             * @return underlying step handler with jacobians
    566             */
    567            public StepHandlerWithJacobians getHandler() {
    568                return handler;
    569            }
    570
    571            /** {@inheritDoc} */
    572            public void handleStep(StepInterpolator interpolator, boolean isLast)
    573                throws DerivativeException {
    574                handler.handleStep(new StepInterpolatorWrapper(interpolator, n, k), isLast);
    575            }
    576
    577            /** {@inheritDoc} */
    578            public boolean requiresDenseOutput() {
    579                return handler.requiresDenseOutput();
    580            }
    581
    582            /** {@inheritDoc} */
    583            public void reset() {
    584                handler.reset();
    585            }
    586
    587        }
    588
    589        /** Wrapper for step interpolators. */
    590        private static class StepInterpolatorWrapper
    591            implements StepInterpolatorWithJacobians {
    592
    593            /** Wrapped interpolator. */
    594            private StepInterpolator interpolator;
    595
    596            /** State array. */
    597            private double[] y;
    598
    599            /** Jacobian with respect to initial state dy/dy0. */
    600            private double[][] dydy0;
    601
    602            /** Jacobian with respect to parameters dy/dp. */
    603            private double[][] dydp;
    604
    605            /** Time derivative of the state array. */
    606            private double[] yDot;
    607
    608            /** Time derivative of the sacobian with respect to initial state dy/dy0. */
    609            private double[][] dydy0Dot;
    610
    611            /** Time derivative of the jacobian with respect to parameters dy/dp. */
    612            private double[][] dydpDot;
    613
    614            /** Simple constructor.
    615             * <p>This constructor is used only for externalization. It does nothing.</p>
    616             */
    617            @SuppressWarnings("unused")
    618            public StepInterpolatorWrapper() {
    619            }
    620
    621            /** Simple constructor.
    622             * @param interpolator wrapped interpolator
    623             * @param n dimension of the original ODE
    624             * @param k number of parameters
    625             */
    626            public StepInterpolatorWrapper(final StepInterpolator interpolator,
    627                                           final int n, final int k) {
    628                this.interpolator = interpolator;
    629                y        = new double[n];
    630                dydy0    = new double[n][n];
    631                dydp     = new double[n][k];
    632                yDot     = new double[n];
    633                dydy0Dot = new double[n][n];
    634                dydpDot  = new double[n][k];
    635            }
    636
    637            /** {@inheritDoc} */
    638            public void setInterpolatedTime(double time) {
    639                interpolator.setInterpolatedTime(time);
    640            }
    641
    642            /** {@inheritDoc} */
    643            public boolean isForward() {
    644                return interpolator.isForward();
    645            }
    646
    647            /** {@inheritDoc} */
    648            public double getPreviousTime() {
    649                return interpolator.getPreviousTime();
    650            }
    651
    652            /** {@inheritDoc} */
    653            public double getInterpolatedTime() {
    654                return interpolator.getInterpolatedTime();
    655            }
    656
    657            /** {@inheritDoc} */
    658            public double[] getInterpolatedY() throws DerivativeException {
    659                double[] extendedState = interpolator.getInterpolatedState();
    660                System.arraycopy(extendedState, 0, y, 0, y.length);
    661                return y;
    662            }
    663
    664            /** {@inheritDoc} */
    665            public double[][] getInterpolatedDyDy0() throws DerivativeException {
    666                double[] extendedState = interpolator.getInterpolatedState();
    667                final int n = y.length;
    668                int start = n;
    669                for (int i = 0; i < n; ++i) {
    670                    System.arraycopy(extendedState, start, dydy0[i], 0, n);
    671                    start += n;
    672                }
    673                return dydy0;
    674            }
    675
    676            /** {@inheritDoc} */
    677            public double[][] getInterpolatedDyDp() throws DerivativeException {
    678                double[] extendedState = interpolator.getInterpolatedState();
    679                final int n = y.length;
    680                final int k = dydp[0].length;
    681                int start = n * (n + 1);
    682                for (int i = 0; i < n; ++i) {
    683                    System.arraycopy(extendedState, start, dydp[i], 0, k);
    684                    start += k;
    685                }
    686                return dydp;
    687            }
    688
    689            /** {@inheritDoc} */
    690            public double[] getInterpolatedYDot() throws DerivativeException {
    691                double[] extendedDerivatives = interpolator.getInterpolatedDerivatives();
    692                System.arraycopy(extendedDerivatives, 0, yDot, 0, yDot.length);
    693                return yDot;
    694            }
    695
    696            /** {@inheritDoc} */
    697            public double[][] getInterpolatedDyDy0Dot() throws DerivativeException {
    698                double[] extendedDerivatives = interpolator.getInterpolatedDerivatives();
    699                final int n = y.length;
    700                int start = n;
    701                for (int i = 0; i < n; ++i) {
    702                    System.arraycopy(extendedDerivatives, start, dydy0Dot[i], 0, n);
    703                    start += n;
    704                }
    705                return dydy0Dot;
    706            }
    707
    708            /** {@inheritDoc} */
    709            public double[][] getInterpolatedDyDpDot() throws DerivativeException {
    710                double[] extendedDerivatives = interpolator.getInterpolatedDerivatives();
    711                final int n = y.length;
    712                final int k = dydpDot[0].length;
    713                int start = n * (n + 1);
    714                for (int i = 0; i < n; ++i) {
    715                    System.arraycopy(extendedDerivatives, start, dydpDot[i], 0, k);
    716                    start += k;
    717                }
    718                return dydpDot;
    719            }
    720
    721            /** {@inheritDoc} */
    722            public double getCurrentTime() {
    723                return interpolator.getCurrentTime();
    724            }
    725
    726            /** {@inheritDoc} */
    727            public StepInterpolatorWithJacobians copy() throws DerivativeException {
    728                final int n = y.length;
    729                final int k = dydp[0].length;
    730                StepInterpolatorWrapper copied =
    731                    new StepInterpolatorWrapper(interpolator.copy(), n, k);
    732                copyArray(y,        copied.y);
    733                copyArray(dydy0,    copied.dydy0);
    734                copyArray(dydp,     copied.dydp);
    735                copyArray(yDot,     copied.yDot);
    736                copyArray(dydy0Dot, copied.dydy0Dot);
    737                copyArray(dydpDot,  copied.dydpDot);
    738                return copied;
    739            }
    740
    741            /** {@inheritDoc} */
    742            public void writeExternal(ObjectOutput out) throws IOException {
    743                out.writeObject(interpolator);
    744                out.writeInt(y.length);
    745                out.writeInt(dydp[0].length);
    746                writeArray(out, y);
    747                writeArray(out, dydy0);
    748                writeArray(out, dydp);
    749                writeArray(out, yDot);
    750                writeArray(out, dydy0Dot);
    751                writeArray(out, dydpDot);
    752            }
    753
    754            /** {@inheritDoc} */
    755            public void readExternal(ObjectInput in) throws IOException, ClassNotFoundException {
    756                interpolator = (StepInterpolator) in.readObject();
    757                final int n = in.readInt();
    758                final int k = in.readInt();
    759                y        = new double[n];
    760                dydy0    = new double[n][n];
    761                dydp     = new double[n][k];
    762                yDot     = new double[n];
    763                dydy0Dot = new double[n][n];
    764                dydpDot  = new double[n][k];
    765                readArray(in, y);
    766                readArray(in, dydy0);
    767                readArray(in, dydp);
    768                readArray(in, yDot);
    769                readArray(in, dydy0Dot);
    770                readArray(in, dydpDot);
    771            }
    772
    773            /** Copy an array.
    774             * @param src source array
    775             * @param dest destination array
    776             */
    777            private static void copyArray(final double[] src, final double[] dest) {
    778                System.arraycopy(src, 0, dest, 0, src.length);
    779            }
    780
    781            /** Copy an array.
    782             * @param src source array
    783             * @param dest destination array
    784             */
    785            private static void copyArray(final double[][] src, final double[][] dest) {
    786                for (int i = 0; i < src.length; ++i) {
    787                    copyArray(src[i], dest[i]);
    788                }
    789            }
    790
    791            /** Write an array.
    792             * @param out output stream
    793             * @param array array to write
    794             * @exception IOException if array cannot be read
    795             */
    796            private static void writeArray(final ObjectOutput out, final double[] array)
    797                throws IOException {
    798                for (int i = 0; i < array.length; ++i) {
    799                    out.writeDouble(array[i]);
    800                }
    801            }
    802
    803            /** Write an array.
    804             * @param out output stream
    805             * @param array array to write
    806             * @exception IOException if array cannot be read
    807             */
    808            private static void writeArray(final ObjectOutput out, final double[][] array)
    809                throws IOException {
    810                for (int i = 0; i < array.length; ++i) {
    811                    writeArray(out, array[i]);
    812                }
    813            }
    814
    815            /** Read an array.
    816             * @param in input stream
    817             * @param array array to read
    818             * @exception IOException if array cannot be read
    819             */
    820            private static void readArray(final ObjectInput in, final double[] array)
    821                throws IOException {
    822                for (int i = 0; i < array.length; ++i) {
    823                    array[i] = in.readDouble();
    824                }
    825            }
    826
    827            /** Read an array.
    828             * @param in input stream
    829             * @param array array to read
    830             * @exception IOException if array cannot be read
    831             */
    832            private static void readArray(final ObjectInput in, final double[][] array)
    833                throws IOException {
    834                for (int i = 0; i < array.length; ++i) {
    835                    readArray(in, array[i]);
    836                }
    837            }
    838
    839        }
    840
    841        /** Wrapper for event handlers. */
    842        private static class EventHandlerWrapper implements EventHandler {
    843
    844            /** Underlying event handler with jacobians. */
    845            private final EventHandlerWithJacobians handler;
    846
    847            /** State array. */
    848            private double[] y;
    849
    850            /** Jacobian with respect to initial state dy/dy0. */
    851            private double[][] dydy0;
    852
    853            /** Jacobian with respect to parameters dy/dp. */
    854            private double[][] dydp;
    855
    856            /** Simple constructor.
    857             * @param handler underlying event handler with jacobians
    858             * @param n dimension of the original ODE
    859             * @param k number of parameters
    860             */
    861            public EventHandlerWrapper(final EventHandlerWithJacobians handler,
    862                                       final int n, final int k) {
    863                this.handler = handler;
    864                y        = new double[n];
    865                dydy0    = new double[n][n];
    866                dydp     = new double[n][k];
    867            }
    868
    869            /** Get the underlying event handler with jacobians.
    870             * @return underlying event handler with jacobians
    871             */
    872            public EventHandlerWithJacobians getHandler() {
    873                return handler;
    874            }
    875
    876            /** {@inheritDoc} */
    877            public int eventOccurred(double t, double[] z, boolean increasing)
    878                throws EventException {
    879                dispatchCompoundState(z, y, dydy0, dydp);
    880                return handler.eventOccurred(t, y, dydy0, dydp, increasing);
    881            }
    882
    883            /** {@inheritDoc} */
    884            public double g(double t, double[] z)
    885                throws EventException {
    886                dispatchCompoundState(z, y, dydy0, dydp);
    887                return handler.g(t, y, dydy0, dydp);
    888            }
    889
    890            /** {@inheritDoc} */
    891            public void resetState(double t, double[] z)
    892                throws EventException {
    893                dispatchCompoundState(z, y, dydy0, dydp);
    894                handler.resetState(t, y, dydy0, dydp);
    895            }
    896
    897        }
    898
    899    }