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.ode.EquationsMapper;
025
026 /** This abstract class represents an interpolator over the last step
027 * during an ODE integration.
028 *
029 * <p>The various ODE integrators provide objects extending this class
030 * to the step handlers. The handlers can use these objects to
031 * retrieve the state vector at intermediate times between the
032 * previous and the current grid points (dense output).</p>
033 *
034 * @see org.apache.commons.math.ode.FirstOrderIntegrator
035 * @see org.apache.commons.math.ode.SecondOrderIntegrator
036 * @see StepHandler
037 *
038 * @version $Id: AbstractStepInterpolator.java 1178235 2011-10-02 19:43:17Z luc $
039 * @since 1.2
040 *
041 */
042
043 public abstract class AbstractStepInterpolator
044 implements StepInterpolator {
045
046 /** current time step */
047 protected double h;
048
049 /** current state */
050 protected double[] currentState;
051
052 /** interpolated time */
053 protected double interpolatedTime;
054
055 /** interpolated state */
056 protected double[] interpolatedState;
057
058 /** interpolated derivatives */
059 protected double[] interpolatedDerivatives;
060
061 /** interpolated primary state */
062 protected double[] interpolatedPrimaryState;
063
064 /** interpolated primary derivatives */
065 protected double[] interpolatedPrimaryDerivatives;
066
067 /** interpolated secondary state */
068 protected double[][] interpolatedSecondaryState;
069
070 /** interpolated secondary derivatives */
071 protected double[][] interpolatedSecondaryDerivatives;
072
073 /** global previous time */
074 private double globalPreviousTime;
075
076 /** global current time */
077 private double globalCurrentTime;
078
079 /** soft previous time */
080 private double softPreviousTime;
081
082 /** soft current time */
083 private double softCurrentTime;
084
085 /** indicate if the step has been finalized or not. */
086 private boolean finalized;
087
088 /** integration direction. */
089 private boolean forward;
090
091 /** indicator for dirty state. */
092 private boolean dirtyState;
093
094 /** Equations mapper for the primary equations set. */
095 private EquationsMapper primaryMapper;
096
097 /** Equations mappers for the secondary equations sets. */
098 private EquationsMapper[] secondaryMappers;
099
100 /** Simple constructor.
101 * This constructor builds an instance that is not usable yet, the
102 * {@link #reinitialize} method should be called before using the
103 * instance in order to initialize the internal arrays. This
104 * constructor is used only in order to delay the initialization in
105 * some cases. As an example, the {@link
106 * org.apache.commons.math.ode.nonstiff.EmbeddedRungeKuttaIntegrator}
107 * class uses the prototyping design pattern to create the step
108 * interpolators by cloning an uninitialized model and latter
109 * initializing the copy.
110 */
111 protected AbstractStepInterpolator() {
112 globalPreviousTime = Double.NaN;
113 globalCurrentTime = Double.NaN;
114 softPreviousTime = Double.NaN;
115 softCurrentTime = Double.NaN;
116 h = Double.NaN;
117 interpolatedTime = Double.NaN;
118 currentState = null;
119 finalized = false;
120 this.forward = true;
121 this.dirtyState = true;
122 primaryMapper = null;
123 secondaryMappers = null;
124 allocateInterpolatedArrays(-1);
125 }
126
127 /** Simple constructor.
128 * @param y reference to the integrator array holding the state at
129 * the end of the step
130 * @param forward integration direction indicator
131 * @param primaryMapper equations mapper for the primary equations set
132 * @param secondaryMappers equations mappers for the secondary equations sets
133 */
134 protected AbstractStepInterpolator(final double[] y, final boolean forward,
135 final EquationsMapper primaryMapper,
136 final EquationsMapper[] secondaryMappers) {
137
138 globalPreviousTime = Double.NaN;
139 globalCurrentTime = Double.NaN;
140 softPreviousTime = Double.NaN;
141 softCurrentTime = Double.NaN;
142 h = Double.NaN;
143 interpolatedTime = Double.NaN;
144 currentState = y;
145 finalized = false;
146 this.forward = forward;
147 this.dirtyState = true;
148 this.primaryMapper = primaryMapper;
149 this.secondaryMappers = (secondaryMappers == null) ? null : secondaryMappers.clone();
150 allocateInterpolatedArrays(y.length);
151
152 }
153
154 /** Copy constructor.
155
156 * <p>The copied interpolator should have been finalized before the
157 * copy, otherwise the copy will not be able to perform correctly
158 * any derivative computation and will throw a {@link
159 * NullPointerException} later. Since we don't want this constructor
160 * to throw the exceptions finalization may involve and since we
161 * don't want this method to modify the state of the copied
162 * interpolator, finalization is <strong>not</strong> done
163 * automatically, it remains under user control.</p>
164
165 * <p>The copy is a deep copy: its arrays are separated from the
166 * original arrays of the instance.</p>
167
168 * @param interpolator interpolator to copy from.
169
170 */
171 protected AbstractStepInterpolator(final AbstractStepInterpolator interpolator) {
172
173 globalPreviousTime = interpolator.globalPreviousTime;
174 globalCurrentTime = interpolator.globalCurrentTime;
175 softPreviousTime = interpolator.softPreviousTime;
176 softCurrentTime = interpolator.softCurrentTime;
177 h = interpolator.h;
178 interpolatedTime = interpolator.interpolatedTime;
179
180 if (interpolator.currentState == null) {
181 currentState = null;
182 primaryMapper = null;
183 secondaryMappers = null;
184 allocateInterpolatedArrays(-1);
185 } else {
186 currentState = interpolator.currentState.clone();
187 interpolatedState = interpolator.interpolatedState.clone();
188 interpolatedDerivatives = interpolator.interpolatedDerivatives.clone();
189 interpolatedPrimaryState = interpolator.interpolatedPrimaryState.clone();
190 interpolatedPrimaryDerivatives = interpolator.interpolatedPrimaryDerivatives.clone();
191 interpolatedSecondaryState = new double[interpolator.interpolatedSecondaryState.length][];
192 interpolatedSecondaryDerivatives = new double[interpolator.interpolatedSecondaryDerivatives.length][];
193 for (int i = 0; i < interpolatedSecondaryState.length; ++i) {
194 interpolatedSecondaryState[i] = interpolator.interpolatedSecondaryState[i].clone();
195 interpolatedSecondaryDerivatives[i] = interpolator.interpolatedSecondaryDerivatives[i].clone();
196 }
197 }
198
199 finalized = interpolator.finalized;
200 forward = interpolator.forward;
201 dirtyState = interpolator.dirtyState;
202 primaryMapper = interpolator.primaryMapper;
203 secondaryMappers = (interpolator.secondaryMappers == null) ?
204 null : interpolator.secondaryMappers.clone();
205
206 }
207
208 /** Allocate the various interpolated states arrays.
209 * @param dimension total dimension (negative if arrays should be set to null)
210 */
211 private void allocateInterpolatedArrays(final int dimension) {
212 if (dimension < 0) {
213 interpolatedState = null;
214 interpolatedDerivatives = null;
215 interpolatedPrimaryState = null;
216 interpolatedPrimaryDerivatives = null;
217 interpolatedSecondaryState = null;
218 interpolatedSecondaryDerivatives = null;
219 } else {
220 interpolatedState = new double[dimension];
221 interpolatedDerivatives = new double[dimension];
222 interpolatedPrimaryState = new double[primaryMapper.getDimension()];
223 interpolatedPrimaryDerivatives = new double[primaryMapper.getDimension()];
224 if (secondaryMappers == null) {
225 interpolatedSecondaryState = null;
226 interpolatedSecondaryDerivatives = null;
227 } else {
228 interpolatedSecondaryState = new double[secondaryMappers.length][];
229 interpolatedSecondaryDerivatives = new double[secondaryMappers.length][];
230 for (int i = 0; i < secondaryMappers.length; ++i) {
231 interpolatedSecondaryState[i] = new double[secondaryMappers[i].getDimension()];
232 interpolatedSecondaryDerivatives[i] = new double[secondaryMappers[i].getDimension()];
233 }
234 }
235 }
236 }
237
238 /** Reinitialize the instance
239 * @param y reference to the integrator array holding the state at the end of the step
240 * @param isForward integration direction indicator
241 * @param primary equations mapper for the primary equations set
242 * @param secondary equations mappers for the secondary equations sets
243 */
244 protected void reinitialize(final double[] y, final boolean isForward,
245 final EquationsMapper primary,
246 final EquationsMapper[] secondary) {
247
248 globalPreviousTime = Double.NaN;
249 globalCurrentTime = Double.NaN;
250 softPreviousTime = Double.NaN;
251 softCurrentTime = Double.NaN;
252 h = Double.NaN;
253 interpolatedTime = Double.NaN;
254 currentState = y;
255 finalized = false;
256 this.forward = isForward;
257 this.dirtyState = true;
258 this.primaryMapper = primary;
259 this.secondaryMappers = secondary.clone();
260 allocateInterpolatedArrays(y.length);
261
262 }
263
264 /** {@inheritDoc} */
265 public StepInterpolator copy() {
266
267 // finalize the step before performing copy
268 finalizeStep();
269
270 // create the new independent instance
271 return doCopy();
272
273 }
274
275 /** Really copy the finalized instance.
276 * <p>This method is called by {@link #copy()} after the
277 * step has been finalized. It must perform a deep copy
278 * to have an new instance completely independent for the
279 * original instance.
280 * @return a copy of the finalized instance
281 */
282 protected abstract StepInterpolator doCopy();
283
284 /** Shift one step forward.
285 * Copy the current time into the previous time, hence preparing the
286 * interpolator for future calls to {@link #storeTime storeTime}
287 */
288 public void shift() {
289 globalPreviousTime = globalCurrentTime;
290 softPreviousTime = globalPreviousTime;
291 softCurrentTime = globalCurrentTime;
292 }
293
294 /** Store the current step time.
295 * @param t current time
296 */
297 public void storeTime(final double t) {
298
299 globalCurrentTime = t;
300 softCurrentTime = globalCurrentTime;
301 h = globalCurrentTime - globalPreviousTime;
302 setInterpolatedTime(t);
303
304 // the step is not finalized anymore
305 finalized = false;
306
307 }
308
309 /** Restrict step range to a limited part of the global step.
310 * <p>
311 * This method can be used to restrict a step and make it appear
312 * as if the original step was smaller. Calling this method
313 * <em>only</em> changes the value returned by {@link #getPreviousTime()},
314 * it does not change any other property
315 * </p>
316 * @param softPreviousTime start of the restricted step
317 * @since 2.2
318 */
319 public void setSoftPreviousTime(final double softPreviousTime) {
320 this.softPreviousTime = softPreviousTime;
321 }
322
323 /** Restrict step range to a limited part of the global step.
324 * <p>
325 * This method can be used to restrict a step and make it appear
326 * as if the original step was smaller. Calling this method
327 * <em>only</em> changes the value returned by {@link #getCurrentTime()},
328 * it does not change any other property
329 * </p>
330 * @param softCurrentTime end of the restricted step
331 * @since 2.2
332 */
333 public void setSoftCurrentTime(final double softCurrentTime) {
334 this.softCurrentTime = softCurrentTime;
335 }
336
337 /**
338 * Get the previous global grid point time.
339 * @return previous global grid point time
340 */
341 public double getGlobalPreviousTime() {
342 return globalPreviousTime;
343 }
344
345 /**
346 * Get the current global grid point time.
347 * @return current global grid point time
348 */
349 public double getGlobalCurrentTime() {
350 return globalCurrentTime;
351 }
352
353 /**
354 * Get the previous soft grid point time.
355 * @return previous soft grid point time
356 * @see #setSoftPreviousTime(double)
357 */
358 public double getPreviousTime() {
359 return softPreviousTime;
360 }
361
362 /**
363 * Get the current soft grid point time.
364 * @return current soft grid point time
365 * @see #setSoftCurrentTime(double)
366 */
367 public double getCurrentTime() {
368 return softCurrentTime;
369 }
370
371 /** {@inheritDoc} */
372 public double getInterpolatedTime() {
373 return interpolatedTime;
374 }
375
376 /** {@inheritDoc} */
377 public void setInterpolatedTime(final double time) {
378 interpolatedTime = time;
379 dirtyState = true;
380 }
381
382 /** {@inheritDoc} */
383 public boolean isForward() {
384 return forward;
385 }
386
387 /** Compute the state and derivatives at the interpolated time.
388 * This is the main processing method that should be implemented by
389 * the derived classes to perform the interpolation.
390 * @param theta normalized interpolation abscissa within the step
391 * (theta is zero at the previous time step and one at the current time step)
392 * @param oneMinusThetaH time gap between the interpolated time and
393 * the current time
394 */
395 protected abstract void computeInterpolatedStateAndDerivatives(double theta,
396 double oneMinusThetaH);
397
398 /** Lazy evaluation of complete interpolated state.
399 */
400 private void evaluateCompleteInterpolatedState() {
401 // lazy evaluation of the state
402 if (dirtyState) {
403 final double oneMinusThetaH = globalCurrentTime - interpolatedTime;
404 final double theta = (h == 0) ? 0 : (h - oneMinusThetaH) / h;
405 computeInterpolatedStateAndDerivatives(theta, oneMinusThetaH);
406 dirtyState = false;
407 }
408 }
409
410 /** {@inheritDoc} */
411 public double[] getInterpolatedState() {
412 evaluateCompleteInterpolatedState();
413 primaryMapper.extractEquationData(interpolatedState,
414 interpolatedPrimaryState);
415 return interpolatedPrimaryState;
416 }
417
418 /** {@inheritDoc} */
419 public double[] getInterpolatedDerivatives() {
420 evaluateCompleteInterpolatedState();
421 primaryMapper.extractEquationData(interpolatedDerivatives,
422 interpolatedPrimaryDerivatives);
423 return interpolatedPrimaryDerivatives;
424 }
425
426 /** {@inheritDoc} */
427 public double[] getInterpolatedSecondaryState(final int index) {
428 evaluateCompleteInterpolatedState();
429 secondaryMappers[index].extractEquationData(interpolatedState,
430 interpolatedSecondaryState[index]);
431 return interpolatedSecondaryState[index];
432 }
433
434 /** {@inheritDoc} */
435 public double[] getInterpolatedSecondaryDerivatives(final int index) {
436 evaluateCompleteInterpolatedState();
437 secondaryMappers[index].extractEquationData(interpolatedDerivatives,
438 interpolatedSecondaryDerivatives[index]);
439 return interpolatedSecondaryDerivatives[index];
440 }
441
442 /**
443 * Finalize the step.
444
445 * <p>Some embedded Runge-Kutta integrators need fewer functions
446 * evaluations than their counterpart step interpolators. These
447 * interpolators should perform the last evaluations they need by
448 * themselves only if they need them. This method triggers these
449 * extra evaluations. It can be called directly by the user step
450 * handler and it is called automatically if {@link
451 * #setInterpolatedTime} is called.</p>
452
453 * <p>Once this method has been called, <strong>no</strong> other
454 * evaluation will be performed on this step. If there is a need to
455 * have some side effects between the step handler and the
456 * differential equations (for example update some data in the
457 * equations once the step has been done), it is advised to call
458 * this method explicitly from the step handler before these side
459 * effects are set up. If the step handler induces no side effect,
460 * then this method can safely be ignored, it will be called
461 * transparently as needed.</p>
462
463 * <p><strong>Warning</strong>: since the step interpolator provided
464 * to the step handler as a parameter of the {@link
465 * StepHandler#handleStep handleStep} is valid only for the duration
466 * of the {@link StepHandler#handleStep handleStep} call, one cannot
467 * simply store a reference and reuse it later. One should first
468 * finalize the instance, then copy this finalized instance into a
469 * new object that can be kept.</p>
470
471 * <p>This method calls the protected <code>doFinalize</code> method
472 * if it has never been called during this step and set a flag
473 * indicating that it has been called once. It is the <code>
474 * doFinalize</code> method which should perform the evaluations.
475 * This wrapping prevents from calling <code>doFinalize</code> several
476 * times and hence evaluating the differential equations too often.
477 * Therefore, subclasses are not allowed not reimplement it, they
478 * should rather reimplement <code>doFinalize</code>.</p>
479
480 */
481 public final void finalizeStep() {
482 if (! finalized) {
483 doFinalize();
484 finalized = true;
485 }
486 }
487
488 /**
489 * Really finalize the step.
490 * The default implementation of this method does nothing.
491 */
492 protected void doFinalize() {
493 }
494
495 /** {@inheritDoc} */
496 public abstract void writeExternal(ObjectOutput out)
497 throws IOException;
498
499 /** {@inheritDoc} */
500 public abstract void readExternal(ObjectInput in)
501 throws IOException, ClassNotFoundException;
502
503 /** Save the base state of the instance.
504 * This method performs step finalization if it has not been done
505 * before.
506 * @param out stream where to save the state
507 * @exception IOException in case of write error
508 */
509 protected void writeBaseExternal(final ObjectOutput out)
510 throws IOException {
511
512 if (currentState == null) {
513 out.writeInt(-1);
514 } else {
515 out.writeInt(currentState.length);
516 }
517 out.writeDouble(globalPreviousTime);
518 out.writeDouble(globalCurrentTime);
519 out.writeDouble(softPreviousTime);
520 out.writeDouble(softCurrentTime);
521 out.writeDouble(h);
522 out.writeBoolean(forward);
523 out.writeObject(primaryMapper);
524 out.write(secondaryMappers.length);
525 for (final EquationsMapper mapper : secondaryMappers) {
526 out.writeObject(mapper);
527 }
528
529 if (currentState != null) {
530 for (int i = 0; i < currentState.length; ++i) {
531 out.writeDouble(currentState[i]);
532 }
533 }
534
535 out.writeDouble(interpolatedTime);
536
537 // we do not store the interpolated state,
538 // it will be recomputed as needed after reading
539
540 // finalize the step (and don't bother saving the now true flag)
541 finalizeStep();
542
543 }
544
545 /** Read the base state of the instance.
546 * This method does <strong>neither</strong> set the interpolated
547 * time nor state. It is up to the derived class to reset it
548 * properly calling the {@link #setInterpolatedTime} method later,
549 * once all rest of the object state has been set up properly.
550 * @param in stream where to read the state from
551 * @return interpolated time to be set later by the caller
552 * @exception IOException in case of read error
553 * @exception ClassNotFoundException if an equation mapper class
554 * cannot be found
555 */
556 protected double readBaseExternal(final ObjectInput in)
557 throws IOException, ClassNotFoundException {
558
559 final int dimension = in.readInt();
560 globalPreviousTime = in.readDouble();
561 globalCurrentTime = in.readDouble();
562 softPreviousTime = in.readDouble();
563 softCurrentTime = in.readDouble();
564 h = in.readDouble();
565 forward = in.readBoolean();
566 primaryMapper = (EquationsMapper) in.readObject();
567 secondaryMappers = new EquationsMapper[in.read()];
568 for (int i = 0; i < secondaryMappers.length; ++i) {
569 secondaryMappers[i] = (EquationsMapper) in.readObject();
570 }
571 dirtyState = true;
572
573 if (dimension < 0) {
574 currentState = null;
575 } else {
576 currentState = new double[dimension];
577 for (int i = 0; i < currentState.length; ++i) {
578 currentState[i] = in.readDouble();
579 }
580 }
581
582 // we do NOT handle the interpolated time and state here
583 interpolatedTime = Double.NaN;
584 allocateInterpolatedArrays(dimension);
585
586 finalized = true;
587
588 return in.readDouble();
589
590 }
591
592 }