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