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
018package org.apache.commons.math3.ode;
019
020import java.io.Serializable;
021import java.util.ArrayList;
022import java.util.List;
023
024import org.apache.commons.math3.exception.DimensionMismatchException;
025import org.apache.commons.math3.exception.MathIllegalArgumentException;
026import org.apache.commons.math3.exception.MaxCountExceededException;
027import org.apache.commons.math3.exception.util.LocalizedFormats;
028import org.apache.commons.math3.ode.sampling.StepHandler;
029import org.apache.commons.math3.ode.sampling.StepInterpolator;
030import org.apache.commons.math3.util.FastMath;
031
032/**
033 * This class stores all information provided by an ODE integrator
034 * during the integration process and build a continuous model of the
035 * solution from this.
036 *
037 * <p>This class act as a step handler from the integrator point of
038 * view. It is called iteratively during the integration process and
039 * stores a copy of all steps information in a sorted collection for
040 * later use. Once the integration process is over, the user can use
041 * the {@link #setInterpolatedTime setInterpolatedTime} and {@link
042 * #getInterpolatedState getInterpolatedState} to retrieve this
043 * information at any time. It is important to wait for the
044 * integration to be over before attempting to call {@link
045 * #setInterpolatedTime setInterpolatedTime} because some internal
046 * variables are set only once the last step has been handled.</p>
047 *
048 * <p>This is useful for example if the main loop of the user
049 * application should remain independent from the integration process
050 * or if one needs to mimic the behaviour of an analytical model
051 * despite a numerical model is used (i.e. one needs the ability to
052 * get the model value at any time or to navigate through the
053 * data).</p>
054 *
055 * <p>If problem modeling is done with several separate
056 * integration phases for contiguous intervals, the same
057 * ContinuousOutputModel can be used as step handler for all
058 * integration phases as long as they are performed in order and in
059 * the same direction. As an example, one can extrapolate the
060 * trajectory of a satellite with one model (i.e. one set of
061 * differential equations) up to the beginning of a maneuver, use
062 * another more complex model including thrusters modeling and
063 * accurate attitude control during the maneuver, and revert to the
064 * first model after the end of the maneuver. If the same continuous
065 * output model handles the steps of all integration phases, the user
066 * do not need to bother when the maneuver begins or ends, he has all
067 * the data available in a transparent manner.</p>
068 *
069 * <p>An important feature of this class is that it implements the
070 * <code>Serializable</code> interface. This means that the result of
071 * an integration can be serialized and reused later (if stored into a
072 * persistent medium like a filesystem or a database) or elsewhere (if
073 * sent to another application). Only the result of the integration is
074 * stored, there is no reference to the integrated problem by
075 * itself.</p>
076 *
077 * <p>One should be aware that the amount of data stored in a
078 * ContinuousOutputModel instance can be important if the state vector
079 * is large, if the integration interval is long or if the steps are
080 * small (which can result from small tolerance settings in {@link
081 * org.apache.commons.math3.ode.nonstiff.AdaptiveStepsizeIntegrator adaptive
082 * step size integrators}).</p>
083 *
084 * @see StepHandler
085 * @see StepInterpolator
086 * @version $Id: ContinuousOutputModel.java 1463684 2013-04-02 19:04:13Z luc $
087 * @since 1.2
088 */
089
090public class ContinuousOutputModel
091  implements StepHandler, Serializable {
092
093    /** Serializable version identifier */
094    private static final long serialVersionUID = -1417964919405031606L;
095
096    /** Initial integration time. */
097    private double initialTime;
098
099    /** Final integration time. */
100    private double finalTime;
101
102    /** Integration direction indicator. */
103    private boolean forward;
104
105    /** Current interpolator index. */
106    private int index;
107
108    /** Steps table. */
109    private List<StepInterpolator> steps;
110
111  /** Simple constructor.
112   * Build an empty continuous output model.
113   */
114  public ContinuousOutputModel() {
115    steps = new ArrayList<StepInterpolator>();
116    initialTime = Double.NaN;
117    finalTime   = Double.NaN;
118    forward     = true;
119    index       = 0;
120  }
121
122  /** Append another model at the end of the instance.
123   * @param model model to add at the end of the instance
124   * @exception MathIllegalArgumentException if the model to append is not
125   * compatible with the instance (dimension of the state vector,
126   * propagation direction, hole between the dates)
127   * @exception MaxCountExceededException if the number of functions evaluations is exceeded
128   * during step finalization
129   */
130  public void append(final ContinuousOutputModel model)
131    throws MathIllegalArgumentException, MaxCountExceededException {
132
133    if (model.steps.size() == 0) {
134      return;
135    }
136
137    if (steps.size() == 0) {
138      initialTime = model.initialTime;
139      forward     = model.forward;
140    } else {
141
142      if (getInterpolatedState().length != model.getInterpolatedState().length) {
143          throw new DimensionMismatchException(model.getInterpolatedState().length,
144                                               getInterpolatedState().length);
145      }
146
147      if (forward ^ model.forward) {
148          throw new MathIllegalArgumentException(LocalizedFormats.PROPAGATION_DIRECTION_MISMATCH);
149      }
150
151      final StepInterpolator lastInterpolator = steps.get(index);
152      final double current  = lastInterpolator.getCurrentTime();
153      final double previous = lastInterpolator.getPreviousTime();
154      final double step = current - previous;
155      final double gap = model.getInitialTime() - current;
156      if (FastMath.abs(gap) > 1.0e-3 * FastMath.abs(step)) {
157        throw new MathIllegalArgumentException(LocalizedFormats.HOLE_BETWEEN_MODELS_TIME_RANGES,
158                                               FastMath.abs(gap));
159      }
160
161    }
162
163    for (StepInterpolator interpolator : model.steps) {
164      steps.add(interpolator.copy());
165    }
166
167    index = steps.size() - 1;
168    finalTime = (steps.get(index)).getCurrentTime();
169
170  }
171
172  /** {@inheritDoc} */
173  public void init(double t0, double[] y0, double t) {
174    initialTime = Double.NaN;
175    finalTime   = Double.NaN;
176    forward     = true;
177    index       = 0;
178    steps.clear();
179  }
180
181  /** Handle the last accepted step.
182   * A copy of the information provided by the last step is stored in
183   * the instance for later use.
184   * @param interpolator interpolator for the last accepted step.
185   * @param isLast true if the step is the last one
186   * @exception MaxCountExceededException if the number of functions evaluations is exceeded
187   * during step finalization
188   */
189  public void handleStep(final StepInterpolator interpolator, final boolean isLast)
190      throws MaxCountExceededException {
191
192    if (steps.size() == 0) {
193      initialTime = interpolator.getPreviousTime();
194      forward     = interpolator.isForward();
195    }
196
197    steps.add(interpolator.copy());
198
199    if (isLast) {
200      finalTime = interpolator.getCurrentTime();
201      index     = steps.size() - 1;
202    }
203
204  }
205
206  /**
207   * Get the initial integration time.
208   * @return initial integration time
209   */
210  public double getInitialTime() {
211    return initialTime;
212  }
213
214  /**
215   * Get the final integration time.
216   * @return final integration time
217   */
218  public double getFinalTime() {
219    return finalTime;
220  }
221
222  /**
223   * Get the time of the interpolated point.
224   * If {@link #setInterpolatedTime} has not been called, it returns
225   * the final integration time.
226   * @return interpolation point time
227   */
228  public double getInterpolatedTime() {
229    return steps.get(index).getInterpolatedTime();
230  }
231
232  /** Set the time of the interpolated point.
233   * <p>This method should <strong>not</strong> be called before the
234   * integration is over because some internal variables are set only
235   * once the last step has been handled.</p>
236   * <p>Setting the time outside of the integration interval is now
237   * allowed (it was not allowed up to version 5.9 of Mantissa), but
238   * should be used with care since the accuracy of the interpolator
239   * will probably be very poor far from this interval. This allowance
240   * has been added to simplify implementation of search algorithms
241   * near the interval endpoints.</p>
242   * @param time time of the interpolated point
243   */
244  public void setInterpolatedTime(final double time) {
245
246      // initialize the search with the complete steps table
247      int iMin = 0;
248      final StepInterpolator sMin = steps.get(iMin);
249      double tMin = 0.5 * (sMin.getPreviousTime() + sMin.getCurrentTime());
250
251      int iMax = steps.size() - 1;
252      final StepInterpolator sMax = steps.get(iMax);
253      double tMax = 0.5 * (sMax.getPreviousTime() + sMax.getCurrentTime());
254
255      // handle points outside of the integration interval
256      // or in the first and last step
257      if (locatePoint(time, sMin) <= 0) {
258        index = iMin;
259        sMin.setInterpolatedTime(time);
260        return;
261      }
262      if (locatePoint(time, sMax) >= 0) {
263        index = iMax;
264        sMax.setInterpolatedTime(time);
265        return;
266      }
267
268      // reduction of the table slice size
269      while (iMax - iMin > 5) {
270
271        // use the last estimated index as the splitting index
272        final StepInterpolator si = steps.get(index);
273        final int location = locatePoint(time, si);
274        if (location < 0) {
275          iMax = index;
276          tMax = 0.5 * (si.getPreviousTime() + si.getCurrentTime());
277        } else if (location > 0) {
278          iMin = index;
279          tMin = 0.5 * (si.getPreviousTime() + si.getCurrentTime());
280        } else {
281          // we have found the target step, no need to continue searching
282          si.setInterpolatedTime(time);
283          return;
284        }
285
286        // compute a new estimate of the index in the reduced table slice
287        final int iMed = (iMin + iMax) / 2;
288        final StepInterpolator sMed = steps.get(iMed);
289        final double tMed = 0.5 * (sMed.getPreviousTime() + sMed.getCurrentTime());
290
291        if ((FastMath.abs(tMed - tMin) < 1e-6) || (FastMath.abs(tMax - tMed) < 1e-6)) {
292          // too close to the bounds, we estimate using a simple dichotomy
293          index = iMed;
294        } else {
295          // estimate the index using a reverse quadratic polynom
296          // (reverse means we have i = P(t), thus allowing to simply
297          // compute index = P(time) rather than solving a quadratic equation)
298          final double d12 = tMax - tMed;
299          final double d23 = tMed - tMin;
300          final double d13 = tMax - tMin;
301          final double dt1 = time - tMax;
302          final double dt2 = time - tMed;
303          final double dt3 = time - tMin;
304          final double iLagrange = ((dt2 * dt3 * d23) * iMax -
305                                    (dt1 * dt3 * d13) * iMed +
306                                    (dt1 * dt2 * d12) * iMin) /
307                                   (d12 * d23 * d13);
308          index = (int) FastMath.rint(iLagrange);
309        }
310
311        // force the next size reduction to be at least one tenth
312        final int low  = FastMath.max(iMin + 1, (9 * iMin + iMax) / 10);
313        final int high = FastMath.min(iMax - 1, (iMin + 9 * iMax) / 10);
314        if (index < low) {
315          index = low;
316        } else if (index > high) {
317          index = high;
318        }
319
320      }
321
322      // now the table slice is very small, we perform an iterative search
323      index = iMin;
324      while ((index <= iMax) && (locatePoint(time, steps.get(index)) > 0)) {
325        ++index;
326      }
327
328      steps.get(index).setInterpolatedTime(time);
329
330  }
331
332  /**
333   * Get the state vector of the interpolated point.
334   * @return state vector at time {@link #getInterpolatedTime}
335   * @exception MaxCountExceededException if the number of functions evaluations is exceeded
336   * @see #getInterpolatedSecondaryState(int)
337   */
338  public double[] getInterpolatedState() throws MaxCountExceededException {
339    return steps.get(index).getInterpolatedState();
340  }
341
342  /** Get the interpolated secondary state corresponding to the secondary equations.
343   * @param secondaryStateIndex index of the secondary set, as returned by {@link
344   * org.apache.commons.math3.ode.ExpandableStatefulODE#addSecondaryEquations(
345   * org.apache.commons.math3.ode.SecondaryEquations)
346   * ExpandableStatefulODE.addSecondaryEquations(SecondaryEquations)}
347   * @return interpolated secondary state at the current interpolation date
348   * @see #getInterpolatedState()
349   * @since 3.2
350   * @exception MaxCountExceededException if the number of functions evaluations is exceeded
351   */
352  public double[] getInterpolatedSecondaryState(final int secondaryStateIndex)
353    throws MaxCountExceededException {
354    return steps.get(index).getInterpolatedSecondaryState(secondaryStateIndex);
355  }
356
357  /** Compare a step interval and a double.
358   * @param time point to locate
359   * @param interval step interval
360   * @return -1 if the double is before the interval, 0 if it is in
361   * the interval, and +1 if it is after the interval, according to
362   * the interval direction
363   */
364  private int locatePoint(final double time, final StepInterpolator interval) {
365    if (forward) {
366      if (time < interval.getPreviousTime()) {
367        return -1;
368      } else if (time > interval.getCurrentTime()) {
369        return +1;
370      } else {
371        return 0;
372      }
373    }
374    if (time > interval.getPreviousTime()) {
375      return -1;
376    } else if (time < interval.getCurrentTime()) {
377      return +1;
378    } else {
379      return 0;
380    }
381  }
382
383}