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