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