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