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