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.nonstiff.AdaptiveStepsizeIntegrator;
026 import org.apache.commons.math.ode.sampling.StepHandler;
027 import org.apache.commons.math.ode.sampling.StepInterpolator;
028
029 /**
030 * This class stores all information provided by an ODE integrator
031 * during the integration process and build a continuous model of the
032 * solution from this.
033 *
034 * <p>This class act as a step handler from the integrator point of
035 * view. It is called iteratively during the integration process and
036 * stores a copy of all steps information in a sorted collection for
037 * later use. Once the integration process is over, the user can use
038 * the {@link #setInterpolatedTime setInterpolatedTime} and {@link
039 * #getInterpolatedState getInterpolatedState} to retrieve this
040 * information at any time. It is important to wait for the
041 * integration to be over before attempting to call {@link
042 * #setInterpolatedTime setInterpolatedTime} because some internal
043 * variables are set only once the last step has been handled.</p>
044 *
045 * <p>This is useful for example if the main loop of the user
046 * application should remain independent from the integration process
047 * or if one needs to mimic the behaviour of an analytical model
048 * despite a numerical model is used (i.e. one needs the ability to
049 * get the model value at any time or to navigate through the
050 * data).</p>
051 *
052 * <p>If problem modeling is done with several separate
053 * integration phases for contiguous intervals, the same
054 * ContinuousOutputModel can be used as step handler for all
055 * integration phases as long as they are performed in order and in
056 * the same direction. As an example, one can extrapolate the
057 * trajectory of a satellite with one model (i.e. one set of
058 * differential equations) up to the beginning of a maneuver, use
059 * another more complex model including thrusters modeling and
060 * accurate attitude control during the maneuver, and revert to the
061 * first model after the end of the maneuver. If the same continuous
062 * output model handles the steps of all integration phases, the user
063 * do not need to bother when the maneuver begins or ends, he has all
064 * the data available in a transparent manner.</p>
065 *
066 * <p>An important feature of this class is that it implements the
067 * <code>Serializable</code> interface. This means that the result of
068 * an integration can be serialized and reused later (if stored into a
069 * persistent medium like a filesystem or a database) or elsewhere (if
070 * sent to another application). Only the result of the integration is
071 * stored, there is no reference to the integrated problem by
072 * itself.</p>
073 *
074 * <p>One should be aware that the amount of data stored in a
075 * ContinuousOutputModel instance can be important if the state vector
076 * is large, if the integration interval is long or if the steps are
077 * small (which can result from small tolerance settings in {@link
078 * AdaptiveStepsizeIntegrator adaptive step size integrators}).</p>
079 *
080 * @see StepHandler
081 * @see StepInterpolator
082 * @version $Revision: 782431 $ $Date: 2009-06-07 15:04:37 -0400 (Sun, 07 Jun 2009) $
083 * @since 1.2
084 */
085
086 public class ContinuousOutputModel
087 implements StepHandler, Serializable {
088
089 /** Simple constructor.
090 * Build an empty continuous output model.
091 */
092 public ContinuousOutputModel() {
093 steps = new ArrayList<StepInterpolator>();
094 reset();
095 }
096
097 /** Append another model at the end of the instance.
098 * @param model model to add at the end of the instance
099 * @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 }