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 * @since 1.2
087 */
088
089public class ContinuousOutputModel
090  implements StepHandler, Serializable {
091
092    /** Serializable version identifier */
093    private static final long serialVersionUID = -1417964919405031606L;
094
095    /** Initial integration time. */
096    private double initialTime;
097
098    /** Final integration time. */
099    private double finalTime;
100
101    /** Integration direction indicator. */
102    private boolean forward;
103
104    /** Current interpolator index. */
105    private int index;
106
107    /** Steps table. */
108    private List<StepInterpolator> steps;
109
110  /** Simple constructor.
111   * Build an empty continuous output model.
112   */
113  public ContinuousOutputModel() {
114    steps = new ArrayList<StepInterpolator>();
115    initialTime = Double.NaN;
116    finalTime   = Double.NaN;
117    forward     = true;
118    index       = 0;
119  }
120
121  /** Append another model at the end of the instance.
122   * @param model model to add at the end of the instance
123   * @exception MathIllegalArgumentException if the model to append is not
124   * compatible with the instance (dimension of the state vector,
125   * propagation direction, hole between the dates)
126   * @exception MaxCountExceededException if the number of functions evaluations is exceeded
127   * during step finalization
128   */
129  public void append(final ContinuousOutputModel model)
130    throws MathIllegalArgumentException, MaxCountExceededException {
131
132    if (model.steps.size() == 0) {
133      return;
134    }
135
136    if (steps.size() == 0) {
137      initialTime = model.initialTime;
138      forward     = model.forward;
139    } else {
140
141      if (getInterpolatedState().length != model.getInterpolatedState().length) {
142          throw new DimensionMismatchException(model.getInterpolatedState().length,
143                                               getInterpolatedState().length);
144      }
145
146      if (forward ^ model.forward) {
147          throw new MathIllegalArgumentException(LocalizedFormats.PROPAGATION_DIRECTION_MISMATCH);
148      }
149
150      final StepInterpolator lastInterpolator = steps.get(index);
151      final double current  = lastInterpolator.getCurrentTime();
152      final double previous = lastInterpolator.getPreviousTime();
153      final double step = current - previous;
154      final double gap = model.getInitialTime() - current;
155      if (FastMath.abs(gap) > 1.0e-3 * FastMath.abs(step)) {
156        throw new MathIllegalArgumentException(LocalizedFormats.HOLE_BETWEEN_MODELS_TIME_RANGES,
157                                               FastMath.abs(gap));
158      }
159
160    }
161
162    for (StepInterpolator interpolator : model.steps) {
163      steps.add(interpolator.copy());
164    }
165
166    index = steps.size() - 1;
167    finalTime = (steps.get(index)).getCurrentTime();
168
169  }
170
171  /** {@inheritDoc} */
172  public void init(double t0, double[] y0, double t) {
173    initialTime = Double.NaN;
174    finalTime   = Double.NaN;
175    forward     = true;
176    index       = 0;
177    steps.clear();
178  }
179
180  /** Handle the last accepted step.
181   * A copy of the information provided by the last step is stored in
182   * the instance for later use.
183   * @param interpolator interpolator for the last accepted step.
184   * @param isLast true if the step is the last one
185   * @exception MaxCountExceededException if the number of functions evaluations is exceeded
186   * during step finalization
187   */
188  public void handleStep(final StepInterpolator interpolator, final boolean isLast)
189      throws MaxCountExceededException {
190
191    if (steps.size() == 0) {
192      initialTime = interpolator.getPreviousTime();
193      forward     = interpolator.isForward();
194    }
195
196    steps.add(interpolator.copy());
197
198    if (isLast) {
199      finalTime = interpolator.getCurrentTime();
200      index     = steps.size() - 1;
201    }
202
203  }
204
205  /**
206   * Get the initial integration time.
207   * @return initial integration time
208   */
209  public double getInitialTime() {
210    return initialTime;
211  }
212
213  /**
214   * Get the final integration time.
215   * @return final integration time
216   */
217  public double getFinalTime() {
218    return finalTime;
219  }
220
221  /**
222   * Get the time of the interpolated point.
223   * If {@link #setInterpolatedTime} has not been called, it returns
224   * the final integration time.
225   * @return interpolation point time
226   */
227  public double getInterpolatedTime() {
228    return steps.get(index).getInterpolatedTime();
229  }
230
231  /** Set the time of the interpolated point.
232   * <p>This method should <strong>not</strong> be called before the
233   * integration is over because some internal variables are set only
234   * once the last step has been handled.</p>
235   * <p>Setting the time outside of the integration interval is now
236   * allowed (it was not allowed up to version 5.9 of Mantissa), but
237   * should be used with care since the accuracy of the interpolator
238   * will probably be very poor far from this interval. This allowance
239   * has been added to simplify implementation of search algorithms
240   * near the interval endpoints.</p>
241   * @param time time of the interpolated point
242   */
243  public void setInterpolatedTime(final double time) {
244
245      // initialize the search with the complete steps table
246      int iMin = 0;
247      final StepInterpolator sMin = steps.get(iMin);
248      double tMin = 0.5 * (sMin.getPreviousTime() + sMin.getCurrentTime());
249
250      int iMax = steps.size() - 1;
251      final StepInterpolator sMax = steps.get(iMax);
252      double tMax = 0.5 * (sMax.getPreviousTime() + sMax.getCurrentTime());
253
254      // handle points outside of the integration interval
255      // or in the first and last step
256      if (locatePoint(time, sMin) <= 0) {
257        index = iMin;
258        sMin.setInterpolatedTime(time);
259        return;
260      }
261      if (locatePoint(time, sMax) >= 0) {
262        index = iMax;
263        sMax.setInterpolatedTime(time);
264        return;
265      }
266
267      // reduction of the table slice size
268      while (iMax - iMin > 5) {
269
270        // use the last estimated index as the splitting index
271        final StepInterpolator si = steps.get(index);
272        final int location = locatePoint(time, si);
273        if (location < 0) {
274          iMax = index;
275          tMax = 0.5 * (si.getPreviousTime() + si.getCurrentTime());
276        } else if (location > 0) {
277          iMin = index;
278          tMin = 0.5 * (si.getPreviousTime() + si.getCurrentTime());
279        } else {
280          // we have found the target step, no need to continue searching
281          si.setInterpolatedTime(time);
282          return;
283        }
284
285        // compute a new estimate of the index in the reduced table slice
286        final int iMed = (iMin + iMax) / 2;
287        final StepInterpolator sMed = steps.get(iMed);
288        final double tMed = 0.5 * (sMed.getPreviousTime() + sMed.getCurrentTime());
289
290        if ((FastMath.abs(tMed - tMin) < 1e-6) || (FastMath.abs(tMax - tMed) < 1e-6)) {
291          // too close to the bounds, we estimate using a simple dichotomy
292          index = iMed;
293        } else {
294          // estimate the index using a reverse quadratic polynom
295          // (reverse means we have i = P(t), thus allowing to simply
296          // compute index = P(time) rather than solving a quadratic equation)
297          final double d12 = tMax - tMed;
298          final double d23 = tMed - tMin;
299          final double d13 = tMax - tMin;
300          final double dt1 = time - tMax;
301          final double dt2 = time - tMed;
302          final double dt3 = time - tMin;
303          final double iLagrange = ((dt2 * dt3 * d23) * iMax -
304                                    (dt1 * dt3 * d13) * iMed +
305                                    (dt1 * dt2 * d12) * iMin) /
306                                   (d12 * d23 * d13);
307          index = (int) FastMath.rint(iLagrange);
308        }
309
310        // force the next size reduction to be at least one tenth
311        final int low  = FastMath.max(iMin + 1, (9 * iMin + iMax) / 10);
312        final int high = FastMath.min(iMax - 1, (iMin + 9 * iMax) / 10);
313        if (index < low) {
314          index = low;
315        } else if (index > high) {
316          index = high;
317        }
318
319      }
320
321      // now the table slice is very small, we perform an iterative search
322      index = iMin;
323      while ((index <= iMax) && (locatePoint(time, steps.get(index)) > 0)) {
324        ++index;
325      }
326
327      steps.get(index).setInterpolatedTime(time);
328
329  }
330
331  /**
332   * Get the state vector of the interpolated point.
333   * @return state vector at time {@link #getInterpolatedTime}
334   * @exception MaxCountExceededException if the number of functions evaluations is exceeded
335   * @see #getInterpolatedSecondaryState(int)
336   */
337  public double[] getInterpolatedState() throws MaxCountExceededException {
338    return steps.get(index).getInterpolatedState();
339  }
340
341  /** Get the interpolated secondary state corresponding to the secondary equations.
342   * @param secondaryStateIndex index of the secondary set, as returned by {@link
343   * org.apache.commons.math3.ode.ExpandableStatefulODE#addSecondaryEquations(
344   * org.apache.commons.math3.ode.SecondaryEquations)
345   * ExpandableStatefulODE.addSecondaryEquations(SecondaryEquations)}
346   * @return interpolated secondary state at the current interpolation date
347   * @see #getInterpolatedState()
348   * @since 3.2
349   * @exception MaxCountExceededException if the number of functions evaluations is exceeded
350   */
351  public double[] getInterpolatedSecondaryState(final int secondaryStateIndex)
352    throws MaxCountExceededException {
353    return steps.get(index).getInterpolatedSecondaryState(secondaryStateIndex);
354  }
355
356  /** Compare a step interval and a double.
357   * @param time point to locate
358   * @param interval step interval
359   * @return -1 if the double is before the interval, 0 if it is in
360   * the interval, and +1 if it is after the interval, according to
361   * the interval direction
362   */
363  private int locatePoint(final double time, final StepInterpolator interval) {
364    if (forward) {
365      if (time < interval.getPreviousTime()) {
366        return -1;
367      } else if (time > interval.getCurrentTime()) {
368        return +1;
369      } else {
370        return 0;
371      }
372    }
373    if (time > interval.getPreviousTime()) {
374      return -1;
375    } else if (time < interval.getCurrentTime()) {
376      return +1;
377    } else {
378      return 0;
379    }
380  }
381
382}