View Javadoc

1   /*
2    * Licensed to the Apache Software Foundation (ASF) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * The ASF licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *      http://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  
18  package org.apache.commons.math.ode;
19  
20  import java.util.ArrayList;
21  import java.util.List;
22  import java.io.Serializable;
23  
24  import org.apache.commons.math.MathRuntimeException;
25  import org.apache.commons.math.ode.nonstiff.AdaptiveStepsizeIntegrator;
26  import org.apache.commons.math.ode.sampling.StepHandler;
27  import org.apache.commons.math.ode.sampling.StepInterpolator;
28  
29  /**
30   * This class stores all information provided by an ODE integrator
31   * during the integration process and build a continuous model of the
32   * solution from this.
33   *
34   * <p>This class act as a step handler from the integrator point of
35   * view. It is called iteratively during the integration process and
36   * stores a copy of all steps information in a sorted collection for
37   * later use. Once the integration process is over, the user can use
38   * the {@link #setInterpolatedTime setInterpolatedTime} and {@link
39   * #getInterpolatedState getInterpolatedState} to retrieve this
40   * information at any time. It is important to wait for the
41   * integration to be over before attempting to call {@link
42   * #setInterpolatedTime setInterpolatedTime} because some internal
43   * variables are set only once the last step has been handled.</p>
44   *
45   * <p>This is useful for example if the main loop of the user
46   * application should remain independent from the integration process
47   * or if one needs to mimic the behaviour of an analytical model
48   * despite a numerical model is used (i.e. one needs the ability to
49   * get the model value at any time or to navigate through the
50   * data).</p>
51   *
52   * <p>If problem modeling is done with several separate
53   * integration phases for contiguous intervals, the same
54   * ContinuousOutputModel can be used as step handler for all
55   * integration phases as long as they are performed in order and in
56   * the same direction. As an example, one can extrapolate the
57   * trajectory of a satellite with one model (i.e. one set of
58   * differential equations) up to the beginning of a maneuver, use
59   * another more complex model including thrusters modeling and
60   * accurate attitude control during the maneuver, and revert to the
61   * first model after the end of the maneuver. If the same continuous
62   * output model handles the steps of all integration phases, the user
63   * do not need to bother when the maneuver begins or ends, he has all
64   * the data available in a transparent manner.</p>
65   *
66   * <p>An important feature of this class is that it implements the
67   * <code>Serializable</code> interface. This means that the result of
68   * an integration can be serialized and reused later (if stored into a
69   * persistent medium like a filesystem or a database) or elsewhere (if
70   * sent to another application). Only the result of the integration is
71   * stored, there is no reference to the integrated problem by
72   * itself.</p>
73   *
74   * <p>One should be aware that the amount of data stored in a
75   * ContinuousOutputModel instance can be important if the state vector
76   * is large, if the integration interval is long or if the steps are
77   * small (which can result from small tolerance settings in {@link
78   * AdaptiveStepsizeIntegrator adaptive step size integrators}).</p>
79   *
80   * @see StepHandler
81   * @see StepInterpolator
82   * @version $Revision: 782431 $ $Date: 2009-06-07 15:04:37 -0400 (Sun, 07 Jun 2009) $
83   * @since 1.2
84   */
85  
86  public class ContinuousOutputModel
87    implements StepHandler, Serializable {
88  
89    /** Simple constructor.
90     * Build an empty continuous output model.
91     */
92    public ContinuousOutputModel() {
93      steps = new ArrayList<StepInterpolator>();
94      reset();
95    }
96  
97    /** Append another model at the end of the instance.
98     * @param model model to add at the end of the instance
99     * @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 }