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.DerivativeException;
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 $Revision: 1073158 $ $Date: 2011-02-21 22:46:52 +0100 (lun. 21 f??vr. 2011) $
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 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 */
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 */
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 */
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 */
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 */
221 public double getInitialTime() {
222 return initialTime;
223 }
224
225 /**
226 * Get the final integration time.
227 * @return final integration time
228 */
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 */
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 */
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 */
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 */
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 }