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