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