diff options
Diffstat (limited to 'src/main/java/org/apache/commons/math/ode/jacobians/FirstOrderIntegratorWithJacobians.java')
| -rw-r--r-- | src/main/java/org/apache/commons/math/ode/jacobians/FirstOrderIntegratorWithJacobians.java | 899 |
1 files changed, 899 insertions, 0 deletions
diff --git a/src/main/java/org/apache/commons/math/ode/jacobians/FirstOrderIntegratorWithJacobians.java b/src/main/java/org/apache/commons/math/ode/jacobians/FirstOrderIntegratorWithJacobians.java new file mode 100644 index 0000000..dc5ca18 --- /dev/null +++ b/src/main/java/org/apache/commons/math/ode/jacobians/FirstOrderIntegratorWithJacobians.java @@ -0,0 +1,899 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +package org.apache.commons.math.ode.jacobians; + +import java.io.IOException; +import java.io.ObjectInput; +import java.io.ObjectOutput; +import java.lang.reflect.Array; +import java.util.ArrayList; +import java.util.Collection; + +import org.apache.commons.math.MathRuntimeException; +import org.apache.commons.math.MaxEvaluationsExceededException; +import org.apache.commons.math.ode.DerivativeException; +import org.apache.commons.math.exception.util.LocalizedFormats; +import org.apache.commons.math.ode.ExtendedFirstOrderDifferentialEquations; +import org.apache.commons.math.ode.FirstOrderIntegrator; +import org.apache.commons.math.ode.IntegratorException; +import org.apache.commons.math.ode.events.EventException; +import org.apache.commons.math.ode.events.EventHandler; +import org.apache.commons.math.ode.sampling.StepHandler; +import org.apache.commons.math.ode.sampling.StepInterpolator; + +/** This class enhances a first order integrator for differential equations to + * compute also partial derivatives of the solution with respect to initial state + * and parameters. + * <p>In order to compute both the state and its derivatives, the ODE problem + * is extended with jacobians of the raw ODE and the variational equations are + * added to form a new compound problem of higher dimension. If the original ODE + * problem has dimension n and there are p parameters, the compound problem will + * have dimension n × (1 + n + p).</p> + * @see ParameterizedODE + * @see ODEWithJacobians + * @version $Revision: 1073158 $ $Date: 2011-02-21 22:46:52 +0100 (lun. 21 févr. 2011) $ + * @since 2.1 + * @deprecated as of 2.2 the complete package is deprecated, it will be replaced + * in 3.0 by a completely rewritten implementation + */ +@Deprecated +public class FirstOrderIntegratorWithJacobians { + + /** Underlying integrator for compound problem. */ + private final FirstOrderIntegrator integrator; + + /** Raw equations to integrate. */ + private final ODEWithJacobians ode; + + /** Maximal number of evaluations allowed. */ + private int maxEvaluations; + + /** Number of evaluations already performed. */ + private int evaluations; + + /** Build an enhanced integrator using internal differentiation to compute jacobians. + * @param integrator underlying integrator to solve the compound problem + * @param ode original problem (f in the equation y' = f(t, y)) + * @param p parameters array (may be null if {@link + * ParameterizedODE#getParametersDimension() + * getParametersDimension()} from original problem is zero) + * @param hY step sizes to use for computing the jacobian df/dy, must have the + * same dimension as the original problem + * @param hP step sizes to use for computing the jacobian df/dp, must have the + * same dimension as the original problem parameters dimension + * @see #FirstOrderIntegratorWithJacobians(FirstOrderIntegrator, + * ODEWithJacobians) + */ + public FirstOrderIntegratorWithJacobians(final FirstOrderIntegrator integrator, + final ParameterizedODE ode, + final double[] p, final double[] hY, final double[] hP) { + checkDimension(ode.getDimension(), hY); + checkDimension(ode.getParametersDimension(), p); + checkDimension(ode.getParametersDimension(), hP); + this.integrator = integrator; + this.ode = new FiniteDifferencesWrapper(ode, p, hY, hP); + setMaxEvaluations(-1); + } + + /** Build an enhanced integrator using ODE builtin jacobian computation features. + * @param integrator underlying integrator to solve the compound problem + * @param ode original problem, which can compute the jacobians by itself + * @see #FirstOrderIntegratorWithJacobians(FirstOrderIntegrator, + * ParameterizedODE, double[], double[], double[]) + */ + public FirstOrderIntegratorWithJacobians(final FirstOrderIntegrator integrator, + final ODEWithJacobians ode) { + this.integrator = integrator; + this.ode = ode; + setMaxEvaluations(-1); + } + + /** Add a step handler to this integrator. + * <p>The handler will be called by the integrator for each accepted + * step.</p> + * @param handler handler for the accepted steps + * @see #getStepHandlers() + * @see #clearStepHandlers() + */ + public void addStepHandler(StepHandlerWithJacobians handler) { + final int n = ode.getDimension(); + final int k = ode.getParametersDimension(); + integrator.addStepHandler(new StepHandlerWrapper(handler, n, k)); + } + + /** Get all the step handlers that have been added to the integrator. + * @return an unmodifiable collection of the added events handlers + * @see #addStepHandler(StepHandlerWithJacobians) + * @see #clearStepHandlers() + */ + public Collection<StepHandlerWithJacobians> getStepHandlers() { + final Collection<StepHandlerWithJacobians> handlers = + new ArrayList<StepHandlerWithJacobians>(); + for (final StepHandler handler : integrator.getStepHandlers()) { + if (handler instanceof StepHandlerWrapper) { + handlers.add(((StepHandlerWrapper) handler).getHandler()); + } + } + return handlers; + } + + /** Remove all the step handlers that have been added to the integrator. + * @see #addStepHandler(StepHandlerWithJacobians) + * @see #getStepHandlers() + */ + public void clearStepHandlers() { + integrator.clearStepHandlers(); + } + + /** Add an event handler to the integrator. + * @param handler event handler + * @param maxCheckInterval maximal time interval between switching + * function checks (this interval prevents missing sign changes in + * case the integration steps becomes very large) + * @param convergence convergence threshold in the event time search + * @param maxIterationCount upper limit of the iteration count in + * the event time search + * @see #getEventHandlers() + * @see #clearEventHandlers() + */ + public void addEventHandler(EventHandlerWithJacobians handler, + double maxCheckInterval, + double convergence, + int maxIterationCount) { + final int n = ode.getDimension(); + final int k = ode.getParametersDimension(); + integrator.addEventHandler(new EventHandlerWrapper(handler, n, k), + maxCheckInterval, convergence, maxIterationCount); + } + + /** Get all the event handlers that have been added to the integrator. + * @return an unmodifiable collection of the added events handlers + * @see #addEventHandler(EventHandlerWithJacobians, double, double, int) + * @see #clearEventHandlers() + */ + public Collection<EventHandlerWithJacobians> getEventHandlers() { + final Collection<EventHandlerWithJacobians> handlers = + new ArrayList<EventHandlerWithJacobians>(); + for (final EventHandler handler : integrator.getEventHandlers()) { + if (handler instanceof EventHandlerWrapper) { + handlers.add(((EventHandlerWrapper) handler).getHandler()); + } + } + return handlers; + } + + /** Remove all the event handlers that have been added to the integrator. + * @see #addEventHandler(EventHandlerWithJacobians, double, double, int) + * @see #getEventHandlers() + */ + public void clearEventHandlers() { + integrator.clearEventHandlers(); + } + + /** Integrate the differential equations and the variational equations up to the given time. + * <p>This method solves an Initial Value Problem (IVP) and also computes the derivatives + * of the solution with respect to initial state and parameters. This can be used as + * a basis to solve Boundary Value Problems (BVP).</p> + * <p>Since this method stores some internal state variables made + * available in its public interface during integration ({@link + * #getCurrentSignedStepsize()}), it is <em>not</em> thread-safe.</p> + * @param t0 initial time + * @param y0 initial value of the state vector at t0 + * @param dY0dP initial value of the state vector derivative with respect to the + * parameters at t0 + * @param t target time for the integration + * (can be set to a value smaller than <code>t0</code> for backward integration) + * @param y placeholder where to put the state vector at each successful + * step (and hence at the end of integration), can be the same object as y0 + * @param dYdY0 placeholder where to put the state vector derivative with respect + * to the initial state (dy[i]/dy0[j] is in element array dYdY0[i][j]) at each successful + * step (and hence at the end of integration) + * @param dYdP placeholder where to put the state vector derivative with respect + * to the parameters (dy[i]/dp[j] is in element array dYdP[i][j]) at each successful + * step (and hence at the end of integration) + * @return stop time, will be the same as target time if integration reached its + * target, but may be different if some event handler stops it at some point. + * @throws IntegratorException if the integrator cannot perform integration + * @throws DerivativeException this exception is propagated to the caller if + * the underlying user function triggers one + */ + public double integrate(final double t0, final double[] y0, final double[][] dY0dP, + final double t, final double[] y, + final double[][] dYdY0, final double[][] dYdP) + throws DerivativeException, IntegratorException { + + final int n = ode.getDimension(); + final int k = ode.getParametersDimension(); + checkDimension(n, y0); + checkDimension(n, y); + checkDimension(n, dYdY0); + checkDimension(n, dYdY0[0]); + if (k != 0) { + checkDimension(n, dY0dP); + checkDimension(k, dY0dP[0]); + checkDimension(n, dYdP); + checkDimension(k, dYdP[0]); + } + + // set up initial state, including partial derivatives + // the compound state z contains the raw state y and its derivatives + // with respect to initial state y0 and to parameters p + // y[i] is stored in z[i] + // dy[i]/dy0[j] is stored in z[n + i * n + j] + // dy[i]/dp[j] is stored in z[n * (n + 1) + i * k + j] + final double[] z = new double[n * (1 + n + k)]; + System.arraycopy(y0, 0, z, 0, n); + for (int i = 0; i < n; ++i) { + + // set diagonal element of dy/dy0 to 1.0 at t = t0 + z[i * (1 + n) + n] = 1.0; + + // set initial derivatives with respect to parameters + System.arraycopy(dY0dP[i], 0, z, n * (n + 1) + i * k, k); + + } + + // integrate the compound state variational equations + evaluations = 0; + final double stopTime = integrator.integrate(new MappingWrapper(), t0, z, t, z); + + // dispatch the final compound state into the state and partial derivatives arrays + dispatchCompoundState(z, y, dYdY0, dYdP); + + return stopTime; + + } + + /** Dispatch a compound state array into state and jacobians arrays. + * @param z compound state + * @param y raw state array to fill + * @param dydy0 jacobian array to fill + * @param dydp jacobian array to fill + */ + private static void dispatchCompoundState(final double[] z, final double[] y, + final double[][] dydy0, final double[][] dydp) { + + final int n = y.length; + final int k = dydp[0].length; + + // state + System.arraycopy(z, 0, y, 0, n); + + // jacobian with respect to initial state + for (int i = 0; i < n; ++i) { + System.arraycopy(z, n * (i + 1), dydy0[i], 0, n); + } + + // jacobian with respect to parameters + for (int i = 0; i < n; ++i) { + System.arraycopy(z, n * (n + 1) + i * k, dydp[i], 0, k); + } + + } + + /** Get the current value of the step start time t<sub>i</sub>. + * <p>This method can be called during integration (typically by + * the object implementing the {@link org.apache.commons.math.ode.FirstOrderDifferentialEquations + * differential equations} problem) if the value of the current step that + * is attempted is needed.</p> + * <p>The result is undefined if the method is called outside of + * calls to <code>integrate</code>.</p> + * @return current value of the step start time t<sub>i</sub> + */ + public double getCurrentStepStart() { + return integrator.getCurrentStepStart(); + } + + /** Get the current signed value of the integration stepsize. + * <p>This method can be called during integration (typically by + * the object implementing the {@link org.apache.commons.math.ode.FirstOrderDifferentialEquations + * differential equations} problem) if the signed value of the current stepsize + * that is tried is needed.</p> + * <p>The result is undefined if the method is called outside of + * calls to <code>integrate</code>.</p> + * @return current signed value of the stepsize + */ + public double getCurrentSignedStepsize() { + return integrator.getCurrentSignedStepsize(); + } + + /** Set the maximal number of differential equations function evaluations. + * <p>The purpose of this method is to avoid infinite loops which can occur + * for example when stringent error constraints are set or when lots of + * discrete events are triggered, thus leading to many rejected steps.</p> + * @param maxEvaluations maximal number of function evaluations (negative + * values are silently converted to maximal integer value, thus representing + * almost unlimited evaluations) + */ + public void setMaxEvaluations(int maxEvaluations) { + this.maxEvaluations = (maxEvaluations < 0) ? Integer.MAX_VALUE : maxEvaluations; + } + + /** Get the maximal number of functions evaluations. + * @return maximal number of functions evaluations + */ + public int getMaxEvaluations() { + return maxEvaluations; + } + + /** Get the number of evaluations of the differential equations function. + * <p> + * The number of evaluations corresponds to the last call to the + * <code>integrate</code> method. It is 0 if the method has not been called yet. + * </p> + * @return number of evaluations of the differential equations function + */ + public int getEvaluations() { + return evaluations; + } + + /** Check array dimensions. + * @param expected expected dimension + * @param array (may be null if expected is 0) + * @throws IllegalArgumentException if the array dimension does not match the expected one + */ + private void checkDimension(final int expected, final Object array) + throws IllegalArgumentException { + int arrayDimension = (array == null) ? 0 : Array.getLength(array); + if (arrayDimension != expected) { + throw MathRuntimeException.createIllegalArgumentException( + LocalizedFormats.DIMENSIONS_MISMATCH_SIMPLE, arrayDimension, expected); + } + } + + /** Wrapper class used to map state and jacobians into compound state. */ + private class MappingWrapper implements ExtendedFirstOrderDifferentialEquations { + + /** Current state. */ + private final double[] y; + + /** Time derivative of the current state. */ + private final double[] yDot; + + /** Derivatives of yDot with respect to state. */ + private final double[][] dFdY; + + /** Derivatives of yDot with respect to parameters. */ + private final double[][] dFdP; + + /** Simple constructor. + */ + public MappingWrapper() { + + final int n = ode.getDimension(); + final int k = ode.getParametersDimension(); + y = new double[n]; + yDot = new double[n]; + dFdY = new double[n][n]; + dFdP = new double[n][k]; + + } + + /** {@inheritDoc} */ + public int getDimension() { + final int n = y.length; + final int k = dFdP[0].length; + return n * (1 + n + k); + } + + /** {@inheritDoc} */ + public int getMainSetDimension() { + return ode.getDimension(); + } + + /** {@inheritDoc} */ + public void computeDerivatives(final double t, final double[] z, final double[] zDot) + throws DerivativeException { + + final int n = y.length; + final int k = dFdP[0].length; + + // compute raw ODE and its jacobians: dy/dt, d[dy/dt]/dy0 and d[dy/dt]/dp + System.arraycopy(z, 0, y, 0, n); + if (++evaluations > maxEvaluations) { + throw new DerivativeException(new MaxEvaluationsExceededException(maxEvaluations)); + } + ode.computeDerivatives(t, y, yDot); + ode.computeJacobians(t, y, yDot, dFdY, dFdP); + + // state part of the compound equations + System.arraycopy(yDot, 0, zDot, 0, n); + + // variational equations: from d[dy/dt]/dy0 to d[dy/dy0]/dt + for (int i = 0; i < n; ++i) { + final double[] dFdYi = dFdY[i]; + for (int j = 0; j < n; ++j) { + double s = 0; + final int startIndex = n + j; + int zIndex = startIndex; + for (int l = 0; l < n; ++l) { + s += dFdYi[l] * z[zIndex]; + zIndex += n; + } + zDot[startIndex + i * n] = s; + } + } + + // variational equations: from d[dy/dt]/dy0 and d[dy/dt]/dp to d[dy/dp]/dt + for (int i = 0; i < n; ++i) { + final double[] dFdYi = dFdY[i]; + final double[] dFdPi = dFdP[i]; + for (int j = 0; j < k; ++j) { + double s = dFdPi[j]; + final int startIndex = n * (n + 1) + j; + int zIndex = startIndex; + for (int l = 0; l < n; ++l) { + s += dFdYi[l] * z[zIndex]; + zIndex += k; + } + zDot[startIndex + i * k] = s; + } + } + + } + + } + + /** Wrapper class to compute jacobians by finite differences for ODE which do not compute them themselves. */ + private class FiniteDifferencesWrapper implements ODEWithJacobians { + + /** Raw ODE without jacobians computation. */ + private final ParameterizedODE ode; + + /** Parameters array (may be null if parameters dimension from original problem is zero) */ + private final double[] p; + + /** Step sizes to use for computing the jacobian df/dy. */ + private final double[] hY; + + /** Step sizes to use for computing the jacobian df/dp. */ + private final double[] hP; + + /** Temporary array for state derivatives used to compute jacobians. */ + private final double[] tmpDot; + + /** Simple constructor. + * @param ode original ODE problem, without jacobians computations + * @param p parameters array (may be null if parameters dimension from original problem is zero) + * @param hY step sizes to use for computing the jacobian df/dy + * @param hP step sizes to use for computing the jacobian df/dp + */ + public FiniteDifferencesWrapper(final ParameterizedODE ode, + final double[] p, final double[] hY, final double[] hP) { + this.ode = ode; + this.p = p.clone(); + this.hY = hY.clone(); + this.hP = hP.clone(); + tmpDot = new double[ode.getDimension()]; + } + + /** {@inheritDoc} */ + public int getDimension() { + return ode.getDimension(); + } + + /** {@inheritDoc} */ + public void computeDerivatives(double t, double[] y, double[] yDot) throws DerivativeException { + // this call to computeDerivatives has already been counted, + // we must not increment the counter again + ode.computeDerivatives(t, y, yDot); + } + + /** {@inheritDoc} */ + public int getParametersDimension() { + return ode.getParametersDimension(); + } + + /** {@inheritDoc} */ + public void computeJacobians(double t, double[] y, double[] yDot, + double[][] dFdY, double[][] dFdP) + throws DerivativeException { + + final int n = hY.length; + final int k = hP.length; + + evaluations += n + k; + if (evaluations > maxEvaluations) { + throw new DerivativeException(new MaxEvaluationsExceededException(maxEvaluations)); + } + + // compute df/dy where f is the ODE and y is the state array + for (int j = 0; j < n; ++j) { + final double savedYj = y[j]; + y[j] += hY[j]; + ode.computeDerivatives(t, y, tmpDot); + for (int i = 0; i < n; ++i) { + dFdY[i][j] = (tmpDot[i] - yDot[i]) / hY[j]; + } + y[j] = savedYj; + } + + // compute df/dp where f is the ODE and p is the parameters array + for (int j = 0; j < k; ++j) { + ode.setParameter(j, p[j] + hP[j]); + ode.computeDerivatives(t, y, tmpDot); + for (int i = 0; i < n; ++i) { + dFdP[i][j] = (tmpDot[i] - yDot[i]) / hP[j]; + } + ode.setParameter(j, p[j]); + } + + } + + } + + /** Wrapper for step handlers. */ + private static class StepHandlerWrapper implements StepHandler { + + /** Underlying step handler with jacobians. */ + private final StepHandlerWithJacobians handler; + + /** Dimension of the original ODE. */ + private final int n; + + /** Number of parameters. */ + private final int k; + + /** Simple constructor. + * @param handler underlying step handler with jacobians + * @param n dimension of the original ODE + * @param k number of parameters + */ + public StepHandlerWrapper(final StepHandlerWithJacobians handler, + final int n, final int k) { + this.handler = handler; + this.n = n; + this.k = k; + } + + /** Get the underlying step handler with jacobians. + * @return underlying step handler with jacobians + */ + public StepHandlerWithJacobians getHandler() { + return handler; + } + + /** {@inheritDoc} */ + public void handleStep(StepInterpolator interpolator, boolean isLast) + throws DerivativeException { + handler.handleStep(new StepInterpolatorWrapper(interpolator, n, k), isLast); + } + + /** {@inheritDoc} */ + public boolean requiresDenseOutput() { + return handler.requiresDenseOutput(); + } + + /** {@inheritDoc} */ + public void reset() { + handler.reset(); + } + + } + + /** Wrapper for step interpolators. */ + private static class StepInterpolatorWrapper + implements StepInterpolatorWithJacobians { + + /** Wrapped interpolator. */ + private StepInterpolator interpolator; + + /** State array. */ + private double[] y; + + /** Jacobian with respect to initial state dy/dy0. */ + private double[][] dydy0; + + /** Jacobian with respect to parameters dy/dp. */ + private double[][] dydp; + + /** Time derivative of the state array. */ + private double[] yDot; + + /** Time derivative of the sacobian with respect to initial state dy/dy0. */ + private double[][] dydy0Dot; + + /** Time derivative of the jacobian with respect to parameters dy/dp. */ + private double[][] dydpDot; + + /** Simple constructor. + * <p>This constructor is used only for externalization. It does nothing.</p> + */ + @SuppressWarnings("unused") + public StepInterpolatorWrapper() { + } + + /** Simple constructor. + * @param interpolator wrapped interpolator + * @param n dimension of the original ODE + * @param k number of parameters + */ + public StepInterpolatorWrapper(final StepInterpolator interpolator, + final int n, final int k) { + this.interpolator = interpolator; + y = new double[n]; + dydy0 = new double[n][n]; + dydp = new double[n][k]; + yDot = new double[n]; + dydy0Dot = new double[n][n]; + dydpDot = new double[n][k]; + } + + /** {@inheritDoc} */ + public void setInterpolatedTime(double time) { + interpolator.setInterpolatedTime(time); + } + + /** {@inheritDoc} */ + public boolean isForward() { + return interpolator.isForward(); + } + + /** {@inheritDoc} */ + public double getPreviousTime() { + return interpolator.getPreviousTime(); + } + + /** {@inheritDoc} */ + public double getInterpolatedTime() { + return interpolator.getInterpolatedTime(); + } + + /** {@inheritDoc} */ + public double[] getInterpolatedY() throws DerivativeException { + double[] extendedState = interpolator.getInterpolatedState(); + System.arraycopy(extendedState, 0, y, 0, y.length); + return y; + } + + /** {@inheritDoc} */ + public double[][] getInterpolatedDyDy0() throws DerivativeException { + double[] extendedState = interpolator.getInterpolatedState(); + final int n = y.length; + int start = n; + for (int i = 0; i < n; ++i) { + System.arraycopy(extendedState, start, dydy0[i], 0, n); + start += n; + } + return dydy0; + } + + /** {@inheritDoc} */ + public double[][] getInterpolatedDyDp() throws DerivativeException { + double[] extendedState = interpolator.getInterpolatedState(); + final int n = y.length; + final int k = dydp[0].length; + int start = n * (n + 1); + for (int i = 0; i < n; ++i) { + System.arraycopy(extendedState, start, dydp[i], 0, k); + start += k; + } + return dydp; + } + + /** {@inheritDoc} */ + public double[] getInterpolatedYDot() throws DerivativeException { + double[] extendedDerivatives = interpolator.getInterpolatedDerivatives(); + System.arraycopy(extendedDerivatives, 0, yDot, 0, yDot.length); + return yDot; + } + + /** {@inheritDoc} */ + public double[][] getInterpolatedDyDy0Dot() throws DerivativeException { + double[] extendedDerivatives = interpolator.getInterpolatedDerivatives(); + final int n = y.length; + int start = n; + for (int i = 0; i < n; ++i) { + System.arraycopy(extendedDerivatives, start, dydy0Dot[i], 0, n); + start += n; + } + return dydy0Dot; + } + + /** {@inheritDoc} */ + public double[][] getInterpolatedDyDpDot() throws DerivativeException { + double[] extendedDerivatives = interpolator.getInterpolatedDerivatives(); + final int n = y.length; + final int k = dydpDot[0].length; + int start = n * (n + 1); + for (int i = 0; i < n; ++i) { + System.arraycopy(extendedDerivatives, start, dydpDot[i], 0, k); + start += k; + } + return dydpDot; + } + + /** {@inheritDoc} */ + public double getCurrentTime() { + return interpolator.getCurrentTime(); + } + + /** {@inheritDoc} */ + public StepInterpolatorWithJacobians copy() throws DerivativeException { + final int n = y.length; + final int k = dydp[0].length; + StepInterpolatorWrapper copied = + new StepInterpolatorWrapper(interpolator.copy(), n, k); + copyArray(y, copied.y); + copyArray(dydy0, copied.dydy0); + copyArray(dydp, copied.dydp); + copyArray(yDot, copied.yDot); + copyArray(dydy0Dot, copied.dydy0Dot); + copyArray(dydpDot, copied.dydpDot); + return copied; + } + + /** {@inheritDoc} */ + public void writeExternal(ObjectOutput out) throws IOException { + out.writeObject(interpolator); + out.writeInt(y.length); + out.writeInt(dydp[0].length); + writeArray(out, y); + writeArray(out, dydy0); + writeArray(out, dydp); + writeArray(out, yDot); + writeArray(out, dydy0Dot); + writeArray(out, dydpDot); + } + + /** {@inheritDoc} */ + public void readExternal(ObjectInput in) throws IOException, ClassNotFoundException { + interpolator = (StepInterpolator) in.readObject(); + final int n = in.readInt(); + final int k = in.readInt(); + y = new double[n]; + dydy0 = new double[n][n]; + dydp = new double[n][k]; + yDot = new double[n]; + dydy0Dot = new double[n][n]; + dydpDot = new double[n][k]; + readArray(in, y); + readArray(in, dydy0); + readArray(in, dydp); + readArray(in, yDot); + readArray(in, dydy0Dot); + readArray(in, dydpDot); + } + + /** Copy an array. + * @param src source array + * @param dest destination array + */ + private static void copyArray(final double[] src, final double[] dest) { + System.arraycopy(src, 0, dest, 0, src.length); + } + + /** Copy an array. + * @param src source array + * @param dest destination array + */ + private static void copyArray(final double[][] src, final double[][] dest) { + for (int i = 0; i < src.length; ++i) { + copyArray(src[i], dest[i]); + } + } + + /** Write an array. + * @param out output stream + * @param array array to write + * @exception IOException if array cannot be read + */ + private static void writeArray(final ObjectOutput out, final double[] array) + throws IOException { + for (int i = 0; i < array.length; ++i) { + out.writeDouble(array[i]); + } + } + + /** Write an array. + * @param out output stream + * @param array array to write + * @exception IOException if array cannot be read + */ + private static void writeArray(final ObjectOutput out, final double[][] array) + throws IOException { + for (int i = 0; i < array.length; ++i) { + writeArray(out, array[i]); + } + } + + /** Read an array. + * @param in input stream + * @param array array to read + * @exception IOException if array cannot be read + */ + private static void readArray(final ObjectInput in, final double[] array) + throws IOException { + for (int i = 0; i < array.length; ++i) { + array[i] = in.readDouble(); + } + } + + /** Read an array. + * @param in input stream + * @param array array to read + * @exception IOException if array cannot be read + */ + private static void readArray(final ObjectInput in, final double[][] array) + throws IOException { + for (int i = 0; i < array.length; ++i) { + readArray(in, array[i]); + } + } + + } + + /** Wrapper for event handlers. */ + private static class EventHandlerWrapper implements EventHandler { + + /** Underlying event handler with jacobians. */ + private final EventHandlerWithJacobians handler; + + /** State array. */ + private double[] y; + + /** Jacobian with respect to initial state dy/dy0. */ + private double[][] dydy0; + + /** Jacobian with respect to parameters dy/dp. */ + private double[][] dydp; + + /** Simple constructor. + * @param handler underlying event handler with jacobians + * @param n dimension of the original ODE + * @param k number of parameters + */ + public EventHandlerWrapper(final EventHandlerWithJacobians handler, + final int n, final int k) { + this.handler = handler; + y = new double[n]; + dydy0 = new double[n][n]; + dydp = new double[n][k]; + } + + /** Get the underlying event handler with jacobians. + * @return underlying event handler with jacobians + */ + public EventHandlerWithJacobians getHandler() { + return handler; + } + + /** {@inheritDoc} */ + public int eventOccurred(double t, double[] z, boolean increasing) + throws EventException { + dispatchCompoundState(z, y, dydy0, dydp); + return handler.eventOccurred(t, y, dydy0, dydp, increasing); + } + + /** {@inheritDoc} */ + public double g(double t, double[] z) + throws EventException { + dispatchCompoundState(z, y, dydy0, dydp); + return handler.g(t, y, dydy0, dydp); + } + + /** {@inheritDoc} */ + public void resetState(double t, double[] z) + throws EventException { + dispatchCompoundState(z, y, dydy0, dydp); + handler.resetState(t, y, dydy0, dydp); + } + + } + +} |