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.math4.legacy.ode;
019
020import java.io.Serializable;
021import java.util.ArrayList;
022import java.util.List;
023
024import org.apache.commons.math4.legacy.exception.DimensionMismatchException;
025import org.apache.commons.math4.legacy.exception.MathIllegalArgumentException;
026import org.apache.commons.math4.legacy.exception.MaxCountExceededException;
027import org.apache.commons.math4.legacy.exception.util.LocalizedFormats;
028import org.apache.commons.math4.legacy.ode.sampling.StepHandler;
029import org.apache.commons.math4.legacy.ode.sampling.StepInterpolator;
030import org.apache.commons.math4.core.jdkmath.JdkMath;
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 file system 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.math4.legacy.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<>();
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.isEmpty()) {
133      return;
134    }
135
136    if (steps.isEmpty()) {
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 (JdkMath.abs(gap) > 1.0e-3 * JdkMath.abs(step)) {
156        throw new MathIllegalArgumentException(LocalizedFormats.HOLE_BETWEEN_MODELS_TIME_RANGES,
157                                               JdkMath.abs(gap));
158      }
159    }
160
161    for (StepInterpolator interpolator : model.steps) {
162      steps.add(interpolator.copy());
163    }
164
165    index = steps.size() - 1;
166    finalTime = (steps.get(index)).getCurrentTime();
167  }
168
169  /** {@inheritDoc} */
170  @Override
171  public void init(double t0, double[] y0, double t) {
172    initialTime = Double.NaN;
173    finalTime   = Double.NaN;
174    forward     = true;
175    index       = 0;
176    steps.clear();
177  }
178
179  /** Handle the last accepted step.
180   * A copy of the information provided by the last step is stored in
181   * the instance for later use.
182   * @param interpolator interpolator for the last accepted step.
183   * @param isLast true if the step is the last one
184   * @exception MaxCountExceededException if the number of functions evaluations is exceeded
185   * during step finalization
186   */
187  @Override
188public void handleStep(final StepInterpolator interpolator, final boolean isLast)
189      throws MaxCountExceededException {
190
191    if (steps.isEmpty()) {
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   * Get the initial integration time.
206   * @return initial integration time
207   */
208  public double getInitialTime() {
209    return initialTime;
210  }
211
212  /**
213   * Get the final integration time.
214   * @return final integration time
215   */
216  public double getFinalTime() {
217    return finalTime;
218  }
219
220  /**
221   * Get the time of the interpolated point.
222   * If {@link #setInterpolatedTime} has not been called, it returns
223   * the final integration time.
224   * @return interpolation point time
225   */
226  public double getInterpolatedTime() {
227    return steps.get(index).getInterpolatedTime();
228  }
229
230  /** Set the time of the interpolated point.
231   * <p>This method should <strong>not</strong> be called before the
232   * integration is over because some internal variables are set only
233   * once the last step has been handled.</p>
234   * <p>Setting the time outside of the integration interval is now
235   * allowed, but should be used with care since the accuracy of the
236   * interpolator will probably be very poor far from this interval.
237   * This allowance has been added to simplify implementation of search
238   * algorithms near the interval endpoints.</p>
239   * <p>Note that each time this method is called, the internal arrays
240   * returned in {@link #getInterpolatedState()}, {@link
241   * #getInterpolatedDerivatives()} and {@link #getInterpolatedSecondaryState(int)}
242   * <em>will</em> be overwritten. So if their content must be preserved
243   * across several calls, user must copy them.</p>
244   * @param time time of the interpolated point
245   * @see #getInterpolatedState()
246   * @see #getInterpolatedDerivatives()
247   * @see #getInterpolatedSecondaryState(int)
248   */
249  public void setInterpolatedTime(final double time) {
250
251      // initialize the search with the complete steps table
252      int iMin = 0;
253      final StepInterpolator sMin = steps.get(iMin);
254      double tMin = 0.5 * (sMin.getPreviousTime() + sMin.getCurrentTime());
255
256      int iMax = steps.size() - 1;
257      final StepInterpolator sMax = steps.get(iMax);
258      double tMax = 0.5 * (sMax.getPreviousTime() + sMax.getCurrentTime());
259
260      // handle points outside of the integration interval
261      // or in the first and last step
262      if (locatePoint(time, sMin) <= 0) {
263        index = iMin;
264        sMin.setInterpolatedTime(time);
265        return;
266      }
267      if (locatePoint(time, sMax) >= 0) {
268        index = iMax;
269        sMax.setInterpolatedTime(time);
270        return;
271      }
272
273      // reduction of the table slice size
274      while (iMax - iMin > 5) {
275
276        // use the last estimated index as the splitting index
277        final StepInterpolator si = steps.get(index);
278        final int location = locatePoint(time, si);
279        if (location < 0) {
280          iMax = index;
281          tMax = 0.5 * (si.getPreviousTime() + si.getCurrentTime());
282        } else if (location > 0) {
283          iMin = index;
284          tMin = 0.5 * (si.getPreviousTime() + si.getCurrentTime());
285        } else {
286          // we have found the target step, no need to continue searching
287          si.setInterpolatedTime(time);
288          return;
289        }
290
291        // compute a new estimate of the index in the reduced table slice
292        final int iMed = (iMin + iMax) / 2;
293        final StepInterpolator sMed = steps.get(iMed);
294        final double tMed = 0.5 * (sMed.getPreviousTime() + sMed.getCurrentTime());
295
296        if (JdkMath.abs(tMed - tMin) < 1e-6 || JdkMath.abs(tMax - tMed) < 1e-6) {
297          // too close to the bounds, we estimate using a simple dichotomy
298          index = iMed;
299        } else {
300          // estimate the index using a reverse quadratic polynom
301          // (reverse means we have i = P(t), thus allowing to simply
302          // compute index = P(time) rather than solving a quadratic equation)
303          final double d12 = tMax - tMed;
304          final double d23 = tMed - tMin;
305          final double d13 = tMax - tMin;
306          final double dt1 = time - tMax;
307          final double dt2 = time - tMed;
308          final double dt3 = time - tMin;
309          final double iLagrange = ((dt2 * dt3 * d23) * iMax -
310                                    (dt1 * dt3 * d13) * iMed +
311                                    (dt1 * dt2 * d12) * iMin) /
312                                   (d12 * d23 * d13);
313          index = (int) JdkMath.rint(iLagrange);
314        }
315
316        // force the next size reduction to be at least one tenth
317        final int low  = JdkMath.max(iMin + 1, (9 * iMin + iMax) / 10);
318        final int high = JdkMath.min(iMax - 1, (iMin + 9 * iMax) / 10);
319        if (index < low) {
320          index = low;
321        } else if (index > high) {
322          index = high;
323        }
324      }
325
326      // now the table slice is very small, we perform an iterative search
327      index = iMin;
328      while (index <= iMax && locatePoint(time, steps.get(index)) > 0) {
329        ++index;
330      }
331
332      steps.get(index).setInterpolatedTime(time);
333  }
334
335  /**
336   * Get the state vector of the interpolated point.
337   * <p>The returned vector is a reference to a reused array, so
338   * it should not be modified and it should be copied if it needs
339   * to be preserved across several calls to the associated
340   * {@link #setInterpolatedTime(double)} method.</p>
341   * @return state vector at time {@link #getInterpolatedTime}
342   * @exception MaxCountExceededException if the number of functions evaluations is exceeded
343   * @see #setInterpolatedTime(double)
344   * @see #getInterpolatedDerivatives()
345   * @see #getInterpolatedSecondaryState(int)
346   * @see #getInterpolatedSecondaryDerivatives(int)
347   */
348  public double[] getInterpolatedState() throws MaxCountExceededException {
349    return steps.get(index).getInterpolatedState();
350  }
351
352  /**
353   * Get the derivatives of the state vector of the interpolated point.
354   * <p>The returned vector is a reference to a reused array, so
355   * it should not be modified and it should be copied if it needs
356   * to be preserved across several calls to the associated
357   * {@link #setInterpolatedTime(double)} method.</p>
358   * @return derivatives of the state vector at time {@link #getInterpolatedTime}
359   * @exception MaxCountExceededException if the number of functions evaluations is exceeded
360   * @see #setInterpolatedTime(double)
361   * @see #getInterpolatedState()
362   * @see #getInterpolatedSecondaryState(int)
363   * @see #getInterpolatedSecondaryDerivatives(int)
364   * @since 3.4
365   */
366  public double[] getInterpolatedDerivatives() throws MaxCountExceededException {
367    return steps.get(index).getInterpolatedDerivatives();
368  }
369
370  /** Get the interpolated secondary state corresponding to the secondary equations.
371   * <p>The returned vector is a reference to a reused array, so
372   * it should not be modified and it should be copied if it needs
373   * to be preserved across several calls to the associated
374   * {@link #setInterpolatedTime(double)} method.</p>
375   * @param secondaryStateIndex index of the secondary set, as returned by {@link
376   * org.apache.commons.math4.legacy.ode.ExpandableStatefulODE#addSecondaryEquations(
377   * org.apache.commons.math4.legacy.ode.SecondaryEquations)
378   * ExpandableStatefulODE.addSecondaryEquations(SecondaryEquations)}
379   * @return interpolated secondary state at the current interpolation date
380   * @see #setInterpolatedTime(double)
381   * @see #getInterpolatedState()
382   * @see #getInterpolatedDerivatives()
383   * @see #getInterpolatedSecondaryDerivatives(int)
384   * @since 3.2
385   * @exception MaxCountExceededException if the number of functions evaluations is exceeded
386   */
387  public double[] getInterpolatedSecondaryState(final int secondaryStateIndex)
388    throws MaxCountExceededException {
389    return steps.get(index).getInterpolatedSecondaryState(secondaryStateIndex);
390  }
391
392  /** Get the interpolated secondary derivatives corresponding to the secondary equations.
393   * <p>The returned vector is a reference to a reused array, so
394   * it should not be modified and it should be copied if it needs
395   * to be preserved across several calls to the associated
396   * {@link #setInterpolatedTime(double)} method.</p>
397   * @param secondaryStateIndex index of the secondary set, as returned by {@link
398   * org.apache.commons.math4.legacy.ode.ExpandableStatefulODE#addSecondaryEquations(
399   * org.apache.commons.math4.legacy.ode.SecondaryEquations)
400   * ExpandableStatefulODE.addSecondaryEquations(SecondaryEquations)}
401   * @return interpolated secondary derivatives at the current interpolation date
402   * @see #setInterpolatedTime(double)
403   * @see #getInterpolatedState()
404   * @see #getInterpolatedDerivatives()
405   * @see #getInterpolatedSecondaryState(int)
406   * @since 3.4
407   * @exception MaxCountExceededException if the number of functions evaluations is exceeded
408   */
409  public double[] getInterpolatedSecondaryDerivatives(final int secondaryStateIndex)
410    throws MaxCountExceededException {
411    return steps.get(index).getInterpolatedSecondaryDerivatives(secondaryStateIndex);
412  }
413
414  /** Compare a step interval and a double.
415   * @param time point to locate
416   * @param interval step interval
417   * @return -1 if the double is before the interval, 0 if it is in
418   * the interval, and +1 if it is after the interval, according to
419   * the interval direction
420   */
421  private int locatePoint(final double time, final StepInterpolator interval) {
422    if (forward) {
423      if (time < interval.getPreviousTime()) {
424        return -1;
425      } else if (time > interval.getCurrentTime()) {
426        return +1;
427      } else {
428        return 0;
429      }
430    }
431    if (time > interval.getPreviousTime()) {
432      return -1;
433    } else if (time < interval.getCurrentTime()) {
434      return +1;
435    } else {
436      return 0;
437    }
438  }
439}