• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
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.DerivativeException;
26 import org.apache.commons.math.exception.util.LocalizedFormats;
27 import org.apache.commons.math.ode.sampling.StepHandler;
28 import org.apache.commons.math.ode.sampling.StepInterpolator;
29 import org.apache.commons.math.util.FastMath;
30 
31 /**
32  * This class stores all information provided by an ODE integrator
33  * during the integration process and build a continuous model of the
34  * solution from this.
35  *
36  * <p>This class act as a step handler from the integrator point of
37  * view. It is called iteratively during the integration process and
38  * stores a copy of all steps information in a sorted collection for
39  * later use. Once the integration process is over, the user can use
40  * the {@link #setInterpolatedTime setInterpolatedTime} and {@link
41  * #getInterpolatedState getInterpolatedState} to retrieve this
42  * information at any time. It is important to wait for the
43  * integration to be over before attempting to call {@link
44  * #setInterpolatedTime setInterpolatedTime} because some internal
45  * variables are set only once the last step has been handled.</p>
46  *
47  * <p>This is useful for example if the main loop of the user
48  * application should remain independent from the integration process
49  * or if one needs to mimic the behaviour of an analytical model
50  * despite a numerical model is used (i.e. one needs the ability to
51  * get the model value at any time or to navigate through the
52  * data).</p>
53  *
54  * <p>If problem modeling is done with several separate
55  * integration phases for contiguous intervals, the same
56  * ContinuousOutputModel can be used as step handler for all
57  * integration phases as long as they are performed in order and in
58  * the same direction. As an example, one can extrapolate the
59  * trajectory of a satellite with one model (i.e. one set of
60  * differential equations) up to the beginning of a maneuver, use
61  * another more complex model including thrusters modeling and
62  * accurate attitude control during the maneuver, and revert to the
63  * first model after the end of the maneuver. If the same continuous
64  * output model handles the steps of all integration phases, the user
65  * do not need to bother when the maneuver begins or ends, he has all
66  * the data available in a transparent manner.</p>
67  *
68  * <p>An important feature of this class is that it implements the
69  * <code>Serializable</code> interface. This means that the result of
70  * an integration can be serialized and reused later (if stored into a
71  * persistent medium like a filesystem or a database) or elsewhere (if
72  * sent to another application). Only the result of the integration is
73  * stored, there is no reference to the integrated problem by
74  * itself.</p>
75  *
76  * <p>One should be aware that the amount of data stored in a
77  * ContinuousOutputModel instance can be important if the state vector
78  * is large, if the integration interval is long or if the steps are
79  * small (which can result from small tolerance settings in {@link
80  * org.apache.commons.math.ode.nonstiff.AdaptiveStepsizeIntegrator adaptive
81  * step size integrators}).</p>
82  *
83  * @see StepHandler
84  * @see StepInterpolator
85  * @version $Revision: 1073158 $ $Date: 2011-02-21 22:46:52 +0100 (lun. 21 févr. 2011) $
86  * @since 1.2
87  */
88 
89 public class ContinuousOutputModel
90   implements StepHandler, Serializable {
91 
92     /** Serializable version identifier */
93     private static final long serialVersionUID = -1417964919405031606L;
94 
95     /** Initial integration time. */
96     private double initialTime;
97 
98     /** Final integration time. */
99     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    */
ContinuousOutputModel()113   public ContinuousOutputModel() {
114     steps = new ArrayList<StepInterpolator>();
115     reset();
116   }
117 
118   /** Append another model at the end of the instance.
119    * @param model model to add at the end of the instance
120    * @exception DerivativeException if user code called from step interpolator
121    * finalization triggers one
122    * @exception IllegalArgumentException if the model to append is not
123    * compatible with the instance (dimension of the state vector,
124    * propagation direction, hole between the dates)
125    */
append(final ContinuousOutputModel model)126   public void append(final ContinuousOutputModel model)
127     throws DerivativeException {
128 
129     if (model.steps.size() == 0) {
130       return;
131     }
132 
133     if (steps.size() == 0) {
134       initialTime = model.initialTime;
135       forward     = model.forward;
136     } else {
137 
138       if (getInterpolatedState().length != model.getInterpolatedState().length) {
139           throw MathRuntimeException.createIllegalArgumentException(
140                 LocalizedFormats.DIMENSIONS_MISMATCH_SIMPLE,
141                 getInterpolatedState().length, model.getInterpolatedState().length);
142       }
143 
144       if (forward ^ model.forward) {
145           throw MathRuntimeException.createIllegalArgumentException(
146                 LocalizedFormats.PROPAGATION_DIRECTION_MISMATCH);
147       }
148 
149       final StepInterpolator lastInterpolator = steps.get(index);
150       final double current  = lastInterpolator.getCurrentTime();
151       final double previous = lastInterpolator.getPreviousTime();
152       final double step = current - previous;
153       final double gap = model.getInitialTime() - current;
154       if (FastMath.abs(gap) > 1.0e-3 * FastMath.abs(step)) {
155         throw MathRuntimeException.createIllegalArgumentException(
156               LocalizedFormats.HOLE_BETWEEN_MODELS_TIME_RANGES, FastMath.abs(gap));
157       }
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 
170   /** Determines whether this handler needs dense output.
171    * <p>The essence of this class is to provide dense output over all
172    * steps, hence it requires the internal steps to provide themselves
173    * dense output. The method therefore returns always true.</p>
174    * @return always true
175    */
requiresDenseOutput()176   public boolean requiresDenseOutput() {
177     return true;
178   }
179 
180   /** Reset the step handler.
181    * Initialize the internal data as required before the first step is
182    * handled.
183    */
reset()184   public void reset() {
185     initialTime = Double.NaN;
186     finalTime   = Double.NaN;
187     forward     = true;
188     index       = 0;
189     steps.clear();
190    }
191 
192   /** Handle the last accepted step.
193    * A copy of the information provided by the last step is stored in
194    * the instance for later use.
195    * @param interpolator interpolator for the last accepted step.
196    * @param isLast true if the step is the last one
197    * @exception DerivativeException if user code called from step interpolator
198    * finalization triggers one
199    */
handleStep(final StepInterpolator interpolator, final boolean isLast)200   public void handleStep(final StepInterpolator interpolator, final boolean isLast)
201     throws DerivativeException {
202 
203     if (steps.size() == 0) {
204       initialTime = interpolator.getPreviousTime();
205       forward     = interpolator.isForward();
206     }
207 
208     steps.add(interpolator.copy());
209 
210     if (isLast) {
211       finalTime = interpolator.getCurrentTime();
212       index     = steps.size() - 1;
213     }
214 
215   }
216 
217   /**
218    * Get the initial integration time.
219    * @return initial integration time
220    */
getInitialTime()221   public double getInitialTime() {
222     return initialTime;
223   }
224 
225   /**
226    * Get the final integration time.
227    * @return final integration time
228    */
getFinalTime()229   public double getFinalTime() {
230     return finalTime;
231   }
232 
233   /**
234    * Get the time of the interpolated point.
235    * If {@link #setInterpolatedTime} has not been called, it returns
236    * the final integration time.
237    * @return interpolation point time
238    */
getInterpolatedTime()239   public double getInterpolatedTime() {
240     return steps.get(index).getInterpolatedTime();
241   }
242 
243   /** Set the time of the interpolated point.
244    * <p>This method should <strong>not</strong> be called before the
245    * integration is over because some internal variables are set only
246    * once the last step has been handled.</p>
247    * <p>Setting the time outside of the integration interval is now
248    * allowed (it was not allowed up to version 5.9 of Mantissa), but
249    * should be used with care since the accuracy of the interpolator
250    * will probably be very poor far from this interval. This allowance
251    * has been added to simplify implementation of search algorithms
252    * near the interval endpoints.</p>
253    * @param time time of the interpolated point
254    */
setInterpolatedTime(final double time)255   public void setInterpolatedTime(final double time) {
256 
257       // initialize the search with the complete steps table
258       int iMin = 0;
259       final StepInterpolator sMin = steps.get(iMin);
260       double tMin = 0.5 * (sMin.getPreviousTime() + sMin.getCurrentTime());
261 
262       int iMax = steps.size() - 1;
263       final StepInterpolator sMax = steps.get(iMax);
264       double tMax = 0.5 * (sMax.getPreviousTime() + sMax.getCurrentTime());
265 
266       // handle points outside of the integration interval
267       // or in the first and last step
268       if (locatePoint(time, sMin) <= 0) {
269         index = iMin;
270         sMin.setInterpolatedTime(time);
271         return;
272       }
273       if (locatePoint(time, sMax) >= 0) {
274         index = iMax;
275         sMax.setInterpolatedTime(time);
276         return;
277       }
278 
279       // reduction of the table slice size
280       while (iMax - iMin > 5) {
281 
282         // use the last estimated index as the splitting index
283         final StepInterpolator si = steps.get(index);
284         final int location = locatePoint(time, si);
285         if (location < 0) {
286           iMax = index;
287           tMax = 0.5 * (si.getPreviousTime() + si.getCurrentTime());
288         } else if (location > 0) {
289           iMin = index;
290           tMin = 0.5 * (si.getPreviousTime() + si.getCurrentTime());
291         } else {
292           // we have found the target step, no need to continue searching
293           si.setInterpolatedTime(time);
294           return;
295         }
296 
297         // compute a new estimate of the index in the reduced table slice
298         final int iMed = (iMin + iMax) / 2;
299         final StepInterpolator sMed = steps.get(iMed);
300         final double tMed = 0.5 * (sMed.getPreviousTime() + sMed.getCurrentTime());
301 
302         if ((FastMath.abs(tMed - tMin) < 1e-6) || (FastMath.abs(tMax - tMed) < 1e-6)) {
303           // too close to the bounds, we estimate using a simple dichotomy
304           index = iMed;
305         } else {
306           // estimate the index using a reverse quadratic polynom
307           // (reverse means we have i = P(t), thus allowing to simply
308           // compute index = P(time) rather than solving a quadratic equation)
309           final double d12 = tMax - tMed;
310           final double d23 = tMed - tMin;
311           final double d13 = tMax - tMin;
312           final double dt1 = time - tMax;
313           final double dt2 = time - tMed;
314           final double dt3 = time - tMin;
315           final double iLagrange = ((dt2 * dt3 * d23) * iMax -
316                                     (dt1 * dt3 * d13) * iMed +
317                                     (dt1 * dt2 * d12) * iMin) /
318                                    (d12 * d23 * d13);
319           index = (int) FastMath.rint(iLagrange);
320         }
321 
322         // force the next size reduction to be at least one tenth
323         final int low  = FastMath.max(iMin + 1, (9 * iMin + iMax) / 10);
324         final int high = FastMath.min(iMax - 1, (iMin + 9 * iMax) / 10);
325         if (index < low) {
326           index = low;
327         } else if (index > high) {
328           index = high;
329         }
330 
331       }
332 
333       // now the table slice is very small, we perform an iterative search
334       index = iMin;
335       while ((index <= iMax) && (locatePoint(time, steps.get(index)) > 0)) {
336         ++index;
337       }
338 
339       steps.get(index).setInterpolatedTime(time);
340 
341   }
342 
343   /**
344    * Get the state vector of the interpolated point.
345    * @return state vector at time {@link #getInterpolatedTime}
346    * @exception DerivativeException if user code called from step interpolator
347    * finalization triggers one
348    */
getInterpolatedState()349   public double[] getInterpolatedState() throws DerivativeException {
350     return steps.get(index).getInterpolatedState();
351   }
352 
353   /** Compare a step interval and a double.
354    * @param time point to locate
355    * @param interval step interval
356    * @return -1 if the double is before the interval, 0 if it is in
357    * the interval, and +1 if it is after the interval, according to
358    * the interval direction
359    */
locatePoint(final double time, final StepInterpolator interval)360   private int locatePoint(final double time, final StepInterpolator interval) {
361     if (forward) {
362       if (time < interval.getPreviousTime()) {
363         return -1;
364       } else if (time > interval.getCurrentTime()) {
365         return +1;
366       } else {
367         return 0;
368       }
369     }
370     if (time > interval.getPreviousTime()) {
371       return -1;
372     } else if (time < interval.getCurrentTime()) {
373       return +1;
374     } else {
375       return 0;
376     }
377   }
378 
379 }
380