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