| /* |
| * 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; |
| |
| import java.util.ArrayList; |
| import java.util.List; |
| import java.io.Serializable; |
| |
| import org.apache.commons.math.MathRuntimeException; |
| import org.apache.commons.math.ode.DerivativeException; |
| import org.apache.commons.math.exception.util.LocalizedFormats; |
| import org.apache.commons.math.ode.sampling.StepHandler; |
| import org.apache.commons.math.ode.sampling.StepInterpolator; |
| import org.apache.commons.math.util.FastMath; |
| |
| /** |
| * This class stores all information provided by an ODE integrator |
| * during the integration process and build a continuous model of the |
| * solution from this. |
| * |
| * <p>This class act as a step handler from the integrator point of |
| * view. It is called iteratively during the integration process and |
| * stores a copy of all steps information in a sorted collection for |
| * later use. Once the integration process is over, the user can use |
| * the {@link #setInterpolatedTime setInterpolatedTime} and {@link |
| * #getInterpolatedState getInterpolatedState} to retrieve this |
| * information at any time. It is important to wait for the |
| * integration to be over before attempting to call {@link |
| * #setInterpolatedTime setInterpolatedTime} because some internal |
| * variables are set only once the last step has been handled.</p> |
| * |
| * <p>This is useful for example if the main loop of the user |
| * application should remain independent from the integration process |
| * or if one needs to mimic the behaviour of an analytical model |
| * despite a numerical model is used (i.e. one needs the ability to |
| * get the model value at any time or to navigate through the |
| * data).</p> |
| * |
| * <p>If problem modeling is done with several separate |
| * integration phases for contiguous intervals, the same |
| * ContinuousOutputModel can be used as step handler for all |
| * integration phases as long as they are performed in order and in |
| * the same direction. As an example, one can extrapolate the |
| * trajectory of a satellite with one model (i.e. one set of |
| * differential equations) up to the beginning of a maneuver, use |
| * another more complex model including thrusters modeling and |
| * accurate attitude control during the maneuver, and revert to the |
| * first model after the end of the maneuver. If the same continuous |
| * output model handles the steps of all integration phases, the user |
| * do not need to bother when the maneuver begins or ends, he has all |
| * the data available in a transparent manner.</p> |
| * |
| * <p>An important feature of this class is that it implements the |
| * <code>Serializable</code> interface. This means that the result of |
| * an integration can be serialized and reused later (if stored into a |
| * persistent medium like a filesystem or a database) or elsewhere (if |
| * sent to another application). Only the result of the integration is |
| * stored, there is no reference to the integrated problem by |
| * itself.</p> |
| * |
| * <p>One should be aware that the amount of data stored in a |
| * ContinuousOutputModel instance can be important if the state vector |
| * is large, if the integration interval is long or if the steps are |
| * small (which can result from small tolerance settings in {@link |
| * org.apache.commons.math.ode.nonstiff.AdaptiveStepsizeIntegrator adaptive |
| * step size integrators}).</p> |
| * |
| * @see StepHandler |
| * @see StepInterpolator |
| * @version $Revision: 1073158 $ $Date: 2011-02-21 22:46:52 +0100 (lun. 21 févr. 2011) $ |
| * @since 1.2 |
| */ |
| |
| public class ContinuousOutputModel |
| implements StepHandler, Serializable { |
| |
| /** Serializable version identifier */ |
| private static final long serialVersionUID = -1417964919405031606L; |
| |
| /** Initial integration time. */ |
| private double initialTime; |
| |
| /** Final integration time. */ |
| private double finalTime; |
| |
| /** Integration direction indicator. */ |
| private boolean forward; |
| |
| /** Current interpolator index. */ |
| private int index; |
| |
| /** Steps table. */ |
| private List<StepInterpolator> steps; |
| |
| /** Simple constructor. |
| * Build an empty continuous output model. |
| */ |
| public ContinuousOutputModel() { |
| steps = new ArrayList<StepInterpolator>(); |
| reset(); |
| } |
| |
| /** Append another model at the end of the instance. |
| * @param model model to add at the end of the instance |
| * @exception DerivativeException if user code called from step interpolator |
| * finalization triggers one |
| * @exception IllegalArgumentException if the model to append is not |
| * compatible with the instance (dimension of the state vector, |
| * propagation direction, hole between the dates) |
| */ |
| public void append(final ContinuousOutputModel model) |
| throws DerivativeException { |
| |
| if (model.steps.size() == 0) { |
| return; |
| } |
| |
| if (steps.size() == 0) { |
| initialTime = model.initialTime; |
| forward = model.forward; |
| } else { |
| |
| if (getInterpolatedState().length != model.getInterpolatedState().length) { |
| throw MathRuntimeException.createIllegalArgumentException( |
| LocalizedFormats.DIMENSIONS_MISMATCH_SIMPLE, |
| getInterpolatedState().length, model.getInterpolatedState().length); |
| } |
| |
| if (forward ^ model.forward) { |
| throw MathRuntimeException.createIllegalArgumentException( |
| LocalizedFormats.PROPAGATION_DIRECTION_MISMATCH); |
| } |
| |
| final StepInterpolator lastInterpolator = steps.get(index); |
| final double current = lastInterpolator.getCurrentTime(); |
| final double previous = lastInterpolator.getPreviousTime(); |
| final double step = current - previous; |
| final double gap = model.getInitialTime() - current; |
| if (FastMath.abs(gap) > 1.0e-3 * FastMath.abs(step)) { |
| throw MathRuntimeException.createIllegalArgumentException( |
| LocalizedFormats.HOLE_BETWEEN_MODELS_TIME_RANGES, FastMath.abs(gap)); |
| } |
| |
| } |
| |
| for (StepInterpolator interpolator : model.steps) { |
| steps.add(interpolator.copy()); |
| } |
| |
| index = steps.size() - 1; |
| finalTime = (steps.get(index)).getCurrentTime(); |
| |
| } |
| |
| /** Determines whether this handler needs dense output. |
| * <p>The essence of this class is to provide dense output over all |
| * steps, hence it requires the internal steps to provide themselves |
| * dense output. The method therefore returns always true.</p> |
| * @return always true |
| */ |
| public boolean requiresDenseOutput() { |
| return true; |
| } |
| |
| /** Reset the step handler. |
| * Initialize the internal data as required before the first step is |
| * handled. |
| */ |
| public void reset() { |
| initialTime = Double.NaN; |
| finalTime = Double.NaN; |
| forward = true; |
| index = 0; |
| steps.clear(); |
| } |
| |
| /** Handle the last accepted step. |
| * A copy of the information provided by the last step is stored in |
| * the instance for later use. |
| * @param interpolator interpolator for the last accepted step. |
| * @param isLast true if the step is the last one |
| * @exception DerivativeException if user code called from step interpolator |
| * finalization triggers one |
| */ |
| public void handleStep(final StepInterpolator interpolator, final boolean isLast) |
| throws DerivativeException { |
| |
| if (steps.size() == 0) { |
| initialTime = interpolator.getPreviousTime(); |
| forward = interpolator.isForward(); |
| } |
| |
| steps.add(interpolator.copy()); |
| |
| if (isLast) { |
| finalTime = interpolator.getCurrentTime(); |
| index = steps.size() - 1; |
| } |
| |
| } |
| |
| /** |
| * Get the initial integration time. |
| * @return initial integration time |
| */ |
| public double getInitialTime() { |
| return initialTime; |
| } |
| |
| /** |
| * Get the final integration time. |
| * @return final integration time |
| */ |
| public double getFinalTime() { |
| return finalTime; |
| } |
| |
| /** |
| * Get the time of the interpolated point. |
| * If {@link #setInterpolatedTime} has not been called, it returns |
| * the final integration time. |
| * @return interpolation point time |
| */ |
| public double getInterpolatedTime() { |
| return steps.get(index).getInterpolatedTime(); |
| } |
| |
| /** Set the time of the interpolated point. |
| * <p>This method should <strong>not</strong> be called before the |
| * integration is over because some internal variables are set only |
| * once the last step has been handled.</p> |
| * <p>Setting the time outside of the integration interval is now |
| * allowed (it was not allowed up to version 5.9 of Mantissa), but |
| * should be used with care since the accuracy of the interpolator |
| * will probably be very poor far from this interval. This allowance |
| * has been added to simplify implementation of search algorithms |
| * near the interval endpoints.</p> |
| * @param time time of the interpolated point |
| */ |
| public void setInterpolatedTime(final double time) { |
| |
| // initialize the search with the complete steps table |
| int iMin = 0; |
| final StepInterpolator sMin = steps.get(iMin); |
| double tMin = 0.5 * (sMin.getPreviousTime() + sMin.getCurrentTime()); |
| |
| int iMax = steps.size() - 1; |
| final StepInterpolator sMax = steps.get(iMax); |
| double tMax = 0.5 * (sMax.getPreviousTime() + sMax.getCurrentTime()); |
| |
| // handle points outside of the integration interval |
| // or in the first and last step |
| if (locatePoint(time, sMin) <= 0) { |
| index = iMin; |
| sMin.setInterpolatedTime(time); |
| return; |
| } |
| if (locatePoint(time, sMax) >= 0) { |
| index = iMax; |
| sMax.setInterpolatedTime(time); |
| return; |
| } |
| |
| // reduction of the table slice size |
| while (iMax - iMin > 5) { |
| |
| // use the last estimated index as the splitting index |
| final StepInterpolator si = steps.get(index); |
| final int location = locatePoint(time, si); |
| if (location < 0) { |
| iMax = index; |
| tMax = 0.5 * (si.getPreviousTime() + si.getCurrentTime()); |
| } else if (location > 0) { |
| iMin = index; |
| tMin = 0.5 * (si.getPreviousTime() + si.getCurrentTime()); |
| } else { |
| // we have found the target step, no need to continue searching |
| si.setInterpolatedTime(time); |
| return; |
| } |
| |
| // compute a new estimate of the index in the reduced table slice |
| final int iMed = (iMin + iMax) / 2; |
| final StepInterpolator sMed = steps.get(iMed); |
| final double tMed = 0.5 * (sMed.getPreviousTime() + sMed.getCurrentTime()); |
| |
| if ((FastMath.abs(tMed - tMin) < 1e-6) || (FastMath.abs(tMax - tMed) < 1e-6)) { |
| // too close to the bounds, we estimate using a simple dichotomy |
| index = iMed; |
| } else { |
| // estimate the index using a reverse quadratic polynom |
| // (reverse means we have i = P(t), thus allowing to simply |
| // compute index = P(time) rather than solving a quadratic equation) |
| final double d12 = tMax - tMed; |
| final double d23 = tMed - tMin; |
| final double d13 = tMax - tMin; |
| final double dt1 = time - tMax; |
| final double dt2 = time - tMed; |
| final double dt3 = time - tMin; |
| final double iLagrange = ((dt2 * dt3 * d23) * iMax - |
| (dt1 * dt3 * d13) * iMed + |
| (dt1 * dt2 * d12) * iMin) / |
| (d12 * d23 * d13); |
| index = (int) FastMath.rint(iLagrange); |
| } |
| |
| // force the next size reduction to be at least one tenth |
| final int low = FastMath.max(iMin + 1, (9 * iMin + iMax) / 10); |
| final int high = FastMath.min(iMax - 1, (iMin + 9 * iMax) / 10); |
| if (index < low) { |
| index = low; |
| } else if (index > high) { |
| index = high; |
| } |
| |
| } |
| |
| // now the table slice is very small, we perform an iterative search |
| index = iMin; |
| while ((index <= iMax) && (locatePoint(time, steps.get(index)) > 0)) { |
| ++index; |
| } |
| |
| steps.get(index).setInterpolatedTime(time); |
| |
| } |
| |
| /** |
| * Get the state vector of the interpolated point. |
| * @return state vector at time {@link #getInterpolatedTime} |
| * @exception DerivativeException if user code called from step interpolator |
| * finalization triggers one |
| */ |
| public double[] getInterpolatedState() throws DerivativeException { |
| return steps.get(index).getInterpolatedState(); |
| } |
| |
| /** Compare a step interval and a double. |
| * @param time point to locate |
| * @param interval step interval |
| * @return -1 if the double is before the interval, 0 if it is in |
| * the interval, and +1 if it is after the interval, according to |
| * the interval direction |
| */ |
| private int locatePoint(final double time, final StepInterpolator interval) { |
| if (forward) { |
| if (time < interval.getPreviousTime()) { |
| return -1; |
| } else if (time > interval.getCurrentTime()) { |
| return +1; |
| } else { |
| return 0; |
| } |
| } |
| if (time > interval.getPreviousTime()) { |
| return -1; |
| } else if (time < interval.getCurrentTime()) { |
| return +1; |
| } else { |
| return 0; |
| } |
| } |
| |
| } |