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    }