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.sampling;
019
020 import java.io.IOException;
021 import java.io.ObjectInput;
022 import java.io.ObjectOutput;
023
024 import org.apache.commons.math.MathRuntimeException;
025 import org.apache.commons.math.ode.DerivativeException;
026 import org.apache.commons.math.ode.FirstOrderIntegrator;
027 import org.apache.commons.math.ode.SecondOrderIntegrator;
028 import org.apache.commons.math.ode.nonstiff.EmbeddedRungeKuttaIntegrator;
029
030 /** This abstract class represents an interpolator over the last step
031 * during an ODE integration.
032 *
033 * <p>The various ODE integrators provide objects extending this class
034 * to the step handlers. The handlers can use these objects to
035 * retrieve the state vector at intermediate times between the
036 * previous and the current grid points (dense output).</p>
037 *
038 * @see FirstOrderIntegrator
039 * @see SecondOrderIntegrator
040 * @see StepHandler
041 *
042 * @version $Revision: 789358 $ $Date: 2009-06-29 11:20:22 -0400 (Mon, 29 Jun 2009) $
043 * @since 1.2
044 *
045 */
046
047 public abstract class AbstractStepInterpolator
048 implements StepInterpolator {
049
050 /** previous time */
051 protected double previousTime;
052
053 /** current time */
054 protected double currentTime;
055
056 /** current time step */
057 protected double h;
058
059 /** current state */
060 protected double[] currentState;
061
062 /** interpolated time */
063 protected double interpolatedTime;
064
065 /** interpolated state */
066 protected double[] interpolatedState;
067
068 /** interpolated derivatives */
069 protected double[] interpolatedDerivatives;
070
071 /** indicate if the step has been finalized or not. */
072 private boolean finalized;
073
074 /** integration direction. */
075 private boolean forward;
076
077 /** indicator for dirty state. */
078 private boolean dirtyState;
079
080
081 /** Simple constructor.
082 * This constructor builds an instance that is not usable yet, the
083 * {@link #reinitialize} method should be called before using the
084 * instance in order to initialize the internal arrays. This
085 * constructor is used only in order to delay the initialization in
086 * some cases. As an example, the {@link
087 * EmbeddedRungeKuttaIntegrator} uses the prototyping design pattern
088 * to create the step interpolators by cloning an uninitialized
089 * model and latter initializing the copy.
090 */
091 protected AbstractStepInterpolator() {
092 previousTime = Double.NaN;
093 currentTime = Double.NaN;
094 h = Double.NaN;
095 interpolatedTime = Double.NaN;
096 currentState = null;
097 interpolatedState = null;
098 interpolatedDerivatives = null;
099 finalized = false;
100 this.forward = true;
101 this.dirtyState = true;
102 }
103
104 /** Simple constructor.
105 * @param y reference to the integrator array holding the state at
106 * the end of the step
107 * @param forward integration direction indicator
108 */
109 protected AbstractStepInterpolator(final double[] y, final boolean forward) {
110
111 previousTime = Double.NaN;
112 currentTime = Double.NaN;
113 h = Double.NaN;
114 interpolatedTime = Double.NaN;
115
116 currentState = y;
117 interpolatedState = new double[y.length];
118 interpolatedDerivatives = new double[y.length];
119
120 finalized = false;
121 this.forward = forward;
122 this.dirtyState = true;
123
124 }
125
126 /** Copy constructor.
127
128 * <p>The copied interpolator should have been finalized before the
129 * copy, otherwise the copy will not be able to perform correctly
130 * any derivative computation and will throw a {@link
131 * NullPointerException} later. Since we don't want this constructor
132 * to throw the exceptions finalization may involve and since we
133 * don't want this method to modify the state of the copied
134 * interpolator, finalization is <strong>not</strong> done
135 * automatically, it remains under user control.</p>
136
137 * <p>The copy is a deep copy: its arrays are separated from the
138 * original arrays of the instance.</p>
139
140 * @param interpolator interpolator to copy from.
141
142 */
143 protected AbstractStepInterpolator(final AbstractStepInterpolator interpolator) {
144
145 previousTime = interpolator.previousTime;
146 currentTime = interpolator.currentTime;
147 h = interpolator.h;
148 interpolatedTime = interpolator.interpolatedTime;
149
150 if (interpolator.currentState != null) {
151 currentState = interpolator.currentState.clone();
152 interpolatedState = interpolator.interpolatedState.clone();
153 interpolatedDerivatives = interpolator.interpolatedDerivatives.clone();
154 } else {
155 currentState = null;
156 interpolatedState = null;
157 interpolatedDerivatives = null;
158 }
159
160 finalized = interpolator.finalized;
161 forward = interpolator.forward;
162 dirtyState = interpolator.dirtyState;
163
164 }
165
166 /** Reinitialize the instance
167 * @param y reference to the integrator array holding the state at
168 * the end of the step
169 * @param forward integration direction indicator
170 */
171 protected void reinitialize(final double[] y, final boolean forward) {
172
173 previousTime = Double.NaN;
174 currentTime = Double.NaN;
175 h = Double.NaN;
176 interpolatedTime = Double.NaN;
177
178 currentState = y;
179 interpolatedState = new double[y.length];
180 interpolatedDerivatives = new double[y.length];
181
182 finalized = false;
183 this.forward = forward;
184 this.dirtyState = true;
185
186 }
187
188 /** {@inheritDoc} */
189 public StepInterpolator copy() throws DerivativeException {
190
191 // finalize the step before performing copy
192 finalizeStep();
193
194 // create the new independent instance
195 return doCopy();
196
197 }
198
199 /** Really copy the finalized instance.
200 * <p>This method is called by {@link #copy()} after the
201 * step has been finalized. It must perform a deep copy
202 * to have an new instance completely independent for the
203 * original instance.
204 * @return a copy of the finalized instance
205 */
206 protected abstract StepInterpolator doCopy();
207
208 /** Shift one step forward.
209 * Copy the current time into the previous time, hence preparing the
210 * interpolator for future calls to {@link #storeTime storeTime}
211 */
212 public void shift() {
213 previousTime = currentTime;
214 }
215
216 /** Store the current step time.
217 * @param t current time
218 */
219 public void storeTime(final double t) {
220
221 currentTime = t;
222 h = currentTime - previousTime;
223 setInterpolatedTime(t);
224
225 // the step is not finalized anymore
226 finalized = false;
227
228 }
229
230 /** {@inheritDoc} */
231 public double getPreviousTime() {
232 return previousTime;
233 }
234
235 /** {@inheritDoc} */
236 public double getCurrentTime() {
237 return currentTime;
238 }
239
240 /** {@inheritDoc} */
241 public double getInterpolatedTime() {
242 return interpolatedTime;
243 }
244
245 /** {@inheritDoc} */
246 public void setInterpolatedTime(final double time) {
247 interpolatedTime = time;
248 dirtyState = true;
249 }
250
251 /** {@inheritDoc} */
252 public boolean isForward() {
253 return forward;
254 }
255
256 /** Compute the state and derivatives at the interpolated time.
257 * This is the main processing method that should be implemented by
258 * the derived classes to perform the interpolation.
259 * @param theta normalized interpolation abscissa within the step
260 * (theta is zero at the previous time step and one at the current time step)
261 * @param oneMinusThetaH time gap between the interpolated time and
262 * the current time
263 * @throws DerivativeException this exception is propagated to the caller if the
264 * underlying user function triggers one
265 */
266 protected abstract void computeInterpolatedStateAndDerivatives(double theta,
267 double oneMinusThetaH)
268 throws DerivativeException;
269
270 /** {@inheritDoc} */
271 public double[] getInterpolatedState() throws DerivativeException {
272
273 // lazy evaluation of the state
274 if (dirtyState) {
275 final double oneMinusThetaH = currentTime - interpolatedTime;
276 final double theta = (h == 0) ? 0 : (h - oneMinusThetaH) / h;
277 computeInterpolatedStateAndDerivatives(theta, oneMinusThetaH);
278 dirtyState = false;
279 }
280
281 return interpolatedState;
282
283 }
284
285 /** {@inheritDoc} */
286 public double[] getInterpolatedDerivatives() throws DerivativeException {
287
288 // lazy evaluation of the state
289 if (dirtyState) {
290 final double oneMinusThetaH = currentTime - interpolatedTime;
291 final double theta = (h == 0) ? 0 : (h - oneMinusThetaH) / h;
292 computeInterpolatedStateAndDerivatives(theta, oneMinusThetaH);
293 dirtyState = false;
294 }
295
296 return interpolatedDerivatives;
297
298 }
299
300 /**
301 * Finalize the step.
302
303 * <p>Some embedded Runge-Kutta integrators need fewer functions
304 * evaluations than their counterpart step interpolators. These
305 * interpolators should perform the last evaluations they need by
306 * themselves only if they need them. This method triggers these
307 * extra evaluations. It can be called directly by the user step
308 * handler and it is called automatically if {@link
309 * #setInterpolatedTime} is called.</p>
310
311 * <p>Once this method has been called, <strong>no</strong> other
312 * evaluation will be performed on this step. If there is a need to
313 * have some side effects between the step handler and the
314 * differential equations (for example update some data in the
315 * equations once the step has been done), it is advised to call
316 * this method explicitly from the step handler before these side
317 * effects are set up. If the step handler induces no side effect,
318 * then this method can safely be ignored, it will be called
319 * transparently as needed.</p>
320
321 * <p><strong>Warning</strong>: since the step interpolator provided
322 * to the step handler as a parameter of the {@link
323 * StepHandler#handleStep handleStep} is valid only for the duration
324 * of the {@link StepHandler#handleStep handleStep} call, one cannot
325 * simply store a reference and reuse it later. One should first
326 * finalize the instance, then copy this finalized instance into a
327 * new object that can be kept.</p>
328
329 * <p>This method calls the protected <code>doFinalize</code> method
330 * if it has never been called during this step and set a flag
331 * indicating that it has been called once. It is the <code>
332 * doFinalize</code> method which should perform the evaluations.
333 * This wrapping prevents from calling <code>doFinalize</code> several
334 * times and hence evaluating the differential equations too often.
335 * Therefore, subclasses are not allowed not reimplement it, they
336 * should rather reimplement <code>doFinalize</code>.</p>
337
338 * @throws DerivativeException this exception is propagated to the
339 * caller if the underlying user function triggers one
340 */
341 public final void finalizeStep()
342 throws DerivativeException {
343 if (! finalized) {
344 doFinalize();
345 finalized = true;
346 }
347 }
348
349 /**
350 * Really finalize the step.
351 * The default implementation of this method does nothing.
352 * @throws DerivativeException this exception is propagated to the
353 * caller if the underlying user function triggers one
354 */
355 protected void doFinalize()
356 throws DerivativeException {
357 }
358
359 /** {@inheritDoc} */
360 public abstract void writeExternal(ObjectOutput out)
361 throws IOException;
362
363 /** {@inheritDoc} */
364 public abstract void readExternal(ObjectInput in)
365 throws IOException, ClassNotFoundException;
366
367 /** Save the base state of the instance.
368 * This method performs step finalization if it has not been done
369 * before.
370 * @param out stream where to save the state
371 * @exception IOException in case of write error
372 */
373 protected void writeBaseExternal(final ObjectOutput out)
374 throws IOException {
375
376 if (currentState == null) {
377 out.writeInt(-1);
378 } else {
379 out.writeInt(currentState.length);
380 }
381 out.writeDouble(previousTime);
382 out.writeDouble(currentTime);
383 out.writeDouble(h);
384 out.writeBoolean(forward);
385
386 if (currentState != null) {
387 for (int i = 0; i < currentState.length; ++i) {
388 out.writeDouble(currentState[i]);
389 }
390 }
391
392 out.writeDouble(interpolatedTime);
393
394 // we do not store the interpolated state,
395 // it will be recomputed as needed after reading
396
397 // finalize the step (and don't bother saving the now true flag)
398 try {
399 finalizeStep();
400 } catch (DerivativeException e) {
401 throw MathRuntimeException.createIOException(e);
402 }
403
404 }
405
406 /** Read the base state of the instance.
407 * This method does <strong>neither</strong> set the interpolated
408 * time nor state. It is up to the derived class to reset it
409 * properly calling the {@link #setInterpolatedTime} method later,
410 * once all rest of the object state has been set up properly.
411 * @param in stream where to read the state from
412 * @return interpolated time be set later by the caller
413 * @exception IOException in case of read error
414 */
415 protected double readBaseExternal(final ObjectInput in)
416 throws IOException {
417
418 final int dimension = in.readInt();
419 previousTime = in.readDouble();
420 currentTime = in.readDouble();
421 h = in.readDouble();
422 forward = in.readBoolean();
423 dirtyState = true;
424
425 if (dimension < 0) {
426 currentState = null;
427 } else {
428 currentState = new double[dimension];
429 for (int i = 0; i < currentState.length; ++i) {
430 currentState[i] = in.readDouble();
431 }
432 }
433
434 // we do NOT handle the interpolated time and state here
435 interpolatedTime = Double.NaN;
436 interpolatedState = (dimension < 0) ? null : new double[dimension];
437 interpolatedDerivatives = (dimension < 0) ? null : new double[dimension];
438
439 finalized = true;
440
441 return in.readDouble();
442
443 }
444
445 }