1 /*
2 * Licensed to the Apache Software Foundation (ASF) under one or more
3 * contributor license agreements. See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * The ASF licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License. You may obtain a copy of the License at
8 *
9 * http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17
18 package org.apache.commons.math4.legacy.ode.sampling;
19
20 import java.io.IOException;
21 import java.io.ObjectInput;
22 import java.io.ObjectOutput;
23
24 import org.apache.commons.math4.legacy.exception.MaxCountExceededException;
25 import org.apache.commons.math4.legacy.ode.EquationsMapper;
26
27 /** This abstract class represents an interpolator over the last step
28 * during an ODE integration.
29 *
30 * <p>The various ODE integrators provide objects extending this class
31 * to the step handlers. The handlers can use these objects to
32 * retrieve the state vector at intermediate times between the
33 * previous and the current grid points (dense output).</p>
34 *
35 * @see org.apache.commons.math4.legacy.ode.FirstOrderIntegrator
36 * @see org.apache.commons.math4.legacy.ode.SecondOrderIntegrator
37 * @see StepHandler
38 *
39 * @since 1.2
40 *
41 */
42
43 public abstract class AbstractStepInterpolator
44 implements StepInterpolator {
45
46 /** current time step. */
47 protected double h;
48
49 /** current state. */
50 protected double[] currentState;
51
52 /** interpolated time. */
53 protected double interpolatedTime;
54
55 /** interpolated state. */
56 protected double[] interpolatedState;
57
58 /** interpolated derivatives. */
59 protected double[] interpolatedDerivatives;
60
61 /** interpolated primary state. */
62 protected double[] interpolatedPrimaryState;
63
64 /** interpolated primary derivatives. */
65 protected double[] interpolatedPrimaryDerivatives;
66
67 /** interpolated secondary state. */
68 protected double[][] interpolatedSecondaryState;
69
70 /** interpolated secondary derivatives. */
71 protected double[][] interpolatedSecondaryDerivatives;
72
73 /** global previous time. */
74 private double globalPreviousTime;
75
76 /** global current time. */
77 private double globalCurrentTime;
78
79 /** soft previous time. */
80 private double softPreviousTime;
81
82 /** soft current time. */
83 private double softCurrentTime;
84
85 /** indicate if the step has been finalized or not. */
86 private boolean finalized;
87
88 /** integration direction. */
89 private boolean forward;
90
91 /** indicator for dirty state. */
92 private boolean dirtyState;
93
94 /** Equations mapper for the primary equations set. */
95 private EquationsMapper primaryMapper;
96
97 /** Equations mappers for the secondary equations sets. */
98 private EquationsMapper[] secondaryMappers;
99
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.math4.legacy.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 /** Copy constructor.
154
155 * <p>The copied interpolator should have been finalized before the
156 * copy, otherwise the copy will not be able to perform correctly
157 * any derivative computation and will throw a {@link
158 * NullPointerException} later. Since we don't want this constructor
159 * to throw the exceptions finalization may involve and since we
160 * don't want this method to modify the state of the copied
161 * interpolator, finalization is <strong>not</strong> done
162 * automatically, it remains under user control.</p>
163
164 * <p>The copy is a deep copy: its arrays are separated from the
165 * original arrays of the instance.</p>
166
167 * @param interpolator interpolator to copy from.
168
169 */
170 protected AbstractStepInterpolator(final AbstractStepInterpolator interpolator) {
171
172 globalPreviousTime = interpolator.globalPreviousTime;
173 globalCurrentTime = interpolator.globalCurrentTime;
174 softPreviousTime = interpolator.softPreviousTime;
175 softCurrentTime = interpolator.softCurrentTime;
176 h = interpolator.h;
177 interpolatedTime = interpolator.interpolatedTime;
178
179 if (interpolator.currentState == null) {
180 currentState = null;
181 primaryMapper = null;
182 secondaryMappers = null;
183 allocateInterpolatedArrays(-1);
184 } else {
185 currentState = interpolator.currentState.clone();
186 interpolatedState = interpolator.interpolatedState.clone();
187 interpolatedDerivatives = interpolator.interpolatedDerivatives.clone();
188 interpolatedPrimaryState = interpolator.interpolatedPrimaryState.clone();
189 interpolatedPrimaryDerivatives = interpolator.interpolatedPrimaryDerivatives.clone();
190 interpolatedSecondaryState = new double[interpolator.interpolatedSecondaryState.length][];
191 interpolatedSecondaryDerivatives = new double[interpolator.interpolatedSecondaryDerivatives.length][];
192 for (int i = 0; i < interpolatedSecondaryState.length; ++i) {
193 interpolatedSecondaryState[i] = interpolator.interpolatedSecondaryState[i].clone();
194 interpolatedSecondaryDerivatives[i] = interpolator.interpolatedSecondaryDerivatives[i].clone();
195 }
196 }
197
198 finalized = interpolator.finalized;
199 forward = interpolator.forward;
200 dirtyState = interpolator.dirtyState;
201 primaryMapper = interpolator.primaryMapper;
202 secondaryMappers = (interpolator.secondaryMappers == null) ?
203 null : interpolator.secondaryMappers.clone();
204 }
205
206 /** Allocate the various interpolated states arrays.
207 * @param dimension total dimension (negative if arrays should be set to null)
208 */
209 private void allocateInterpolatedArrays(final int dimension) {
210 if (dimension < 0) {
211 interpolatedState = null;
212 interpolatedDerivatives = null;
213 interpolatedPrimaryState = null;
214 interpolatedPrimaryDerivatives = null;
215 interpolatedSecondaryState = null;
216 interpolatedSecondaryDerivatives = null;
217 } else {
218 interpolatedState = new double[dimension];
219 interpolatedDerivatives = new double[dimension];
220 interpolatedPrimaryState = new double[primaryMapper.getDimension()];
221 interpolatedPrimaryDerivatives = new double[primaryMapper.getDimension()];
222 if (secondaryMappers == null) {
223 interpolatedSecondaryState = null;
224 interpolatedSecondaryDerivatives = null;
225 } else {
226 interpolatedSecondaryState = new double[secondaryMappers.length][];
227 interpolatedSecondaryDerivatives = new double[secondaryMappers.length][];
228 for (int i = 0; i < secondaryMappers.length; ++i) {
229 interpolatedSecondaryState[i] = new double[secondaryMappers[i].getDimension()];
230 interpolatedSecondaryDerivatives[i] = new double[secondaryMappers[i].getDimension()];
231 }
232 }
233 }
234 }
235
236 /** Reinitialize the instance.
237 * @param y reference to the integrator array holding the state at the end of the step
238 * @param isForward integration direction indicator
239 * @param primary equations mapper for the primary equations set
240 * @param secondary equations mappers for the secondary equations sets
241 */
242 protected void reinitialize(final double[] y, final boolean isForward,
243 final EquationsMapper primary,
244 final EquationsMapper[] secondary) {
245
246 globalPreviousTime = Double.NaN;
247 globalCurrentTime = Double.NaN;
248 softPreviousTime = Double.NaN;
249 softCurrentTime = Double.NaN;
250 h = Double.NaN;
251 interpolatedTime = Double.NaN;
252 currentState = y;
253 finalized = false;
254 this.forward = isForward;
255 this.dirtyState = true;
256 this.primaryMapper = primary;
257 this.secondaryMappers = secondary.clone();
258 allocateInterpolatedArrays(y.length);
259 }
260
261 /** {@inheritDoc} */
262 @Override
263 public StepInterpolator copy() throws MaxCountExceededException {
264
265 // finalize the step before performing copy
266 finalizeStep();
267
268 // create the new independent instance
269 return doCopy();
270 }
271
272 /** Really copy the finalized instance.
273 * <p>This method is called by {@link #copy()} after the
274 * step has been finalized. It must perform a deep copy
275 * to have an new instance completely independent for the
276 * original instance.
277 * @return a copy of the finalized instance
278 */
279 protected abstract StepInterpolator doCopy();
280
281 /** Shift one step forward.
282 * Copy the current time into the previous time, hence preparing the
283 * interpolator for future calls to {@link #storeTime storeTime}
284 */
285 public void shift() {
286 globalPreviousTime = globalCurrentTime;
287 softPreviousTime = globalPreviousTime;
288 softCurrentTime = globalCurrentTime;
289 }
290
291 /** Store the current step time.
292 * @param t current time
293 */
294 public void storeTime(final double t) {
295
296 globalCurrentTime = t;
297 softCurrentTime = globalCurrentTime;
298 h = globalCurrentTime - globalPreviousTime;
299 setInterpolatedTime(t);
300
301 // the step is not finalized anymore
302 finalized = false;
303 }
304
305 /** Restrict step range to a limited part of the global step.
306 * <p>
307 * This method can be used to restrict a step and make it appear
308 * as if the original step was smaller. Calling this method
309 * <em>only</em> changes the value returned by {@link #getPreviousTime()},
310 * it does not change any other property
311 * </p>
312 * @param softPreviousTime start of the restricted step
313 * @since 2.2
314 */
315 public void setSoftPreviousTime(final double softPreviousTime) {
316 this.softPreviousTime = softPreviousTime;
317 }
318
319 /** Restrict step range to a limited part of the global step.
320 * <p>
321 * This method can be used to restrict a step and make it appear
322 * as if the original step was smaller. Calling this method
323 * <em>only</em> changes the value returned by {@link #getCurrentTime()},
324 * it does not change any other property
325 * </p>
326 * @param softCurrentTime end of the restricted step
327 * @since 2.2
328 */
329 public void setSoftCurrentTime(final double softCurrentTime) {
330 this.softCurrentTime = softCurrentTime;
331 }
332
333 /**
334 * Get the previous global grid point time.
335 * @return previous global grid point time
336 */
337 public double getGlobalPreviousTime() {
338 return globalPreviousTime;
339 }
340
341 /**
342 * Get the current global grid point time.
343 * @return current global grid point time
344 */
345 public double getGlobalCurrentTime() {
346 return globalCurrentTime;
347 }
348
349 /**
350 * Get the previous soft grid point time.
351 * @return previous soft grid point time
352 * @see #setSoftPreviousTime(double)
353 */
354 @Override
355 public double getPreviousTime() {
356 return softPreviousTime;
357 }
358
359 /**
360 * Get the current soft grid point time.
361 * @return current soft grid point time
362 * @see #setSoftCurrentTime(double)
363 */
364 @Override
365 public double getCurrentTime() {
366 return softCurrentTime;
367 }
368
369 /** {@inheritDoc} */
370 @Override
371 public double getInterpolatedTime() {
372 return interpolatedTime;
373 }
374
375 /** {@inheritDoc} */
376 @Override
377 public void setInterpolatedTime(final double time) {
378 interpolatedTime = time;
379 dirtyState = true;
380 }
381
382 /** {@inheritDoc} */
383 @Override
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 @Override
417 public double[] getInterpolatedState() throws MaxCountExceededException {
418 evaluateCompleteInterpolatedState();
419 primaryMapper.extractEquationData(interpolatedState,
420 interpolatedPrimaryState);
421 return interpolatedPrimaryState;
422 }
423
424 /** {@inheritDoc} */
425 @Override
426 public double[] getInterpolatedDerivatives() throws MaxCountExceededException {
427 evaluateCompleteInterpolatedState();
428 primaryMapper.extractEquationData(interpolatedDerivatives,
429 interpolatedPrimaryDerivatives);
430 return interpolatedPrimaryDerivatives;
431 }
432
433 /** {@inheritDoc} */
434 @Override
435 public double[] getInterpolatedSecondaryState(final int index) throws MaxCountExceededException {
436 evaluateCompleteInterpolatedState();
437 secondaryMappers[index].extractEquationData(interpolatedState,
438 interpolatedSecondaryState[index]);
439 return interpolatedSecondaryState[index];
440 }
441
442 /** {@inheritDoc} */
443 @Override
444 public double[] getInterpolatedSecondaryDerivatives(final int index) throws MaxCountExceededException {
445 evaluateCompleteInterpolatedState();
446 secondaryMappers[index].extractEquationData(interpolatedDerivatives,
447 interpolatedSecondaryDerivatives[index]);
448 return interpolatedSecondaryDerivatives[index];
449 }
450
451 /**
452 * Finalize the step.
453
454 * <p>Some embedded Runge-Kutta integrators need fewer functions
455 * evaluations than their counterpart step interpolators. These
456 * interpolators should perform the last evaluations they need by
457 * themselves only if they need them. This method triggers these
458 * extra evaluations. It can be called directly by the user step
459 * handler and it is called automatically if {@link
460 * #setInterpolatedTime} is called.</p>
461
462 * <p>Once this method has been called, <strong>no</strong> other
463 * evaluation will be performed on this step. If there is a need to
464 * have some side effects between the step handler and the
465 * differential equations (for example update some data in the
466 * equations once the step has been done), it is advised to call
467 * this method explicitly from the step handler before these side
468 * effects are set up. If the step handler induces no side effect,
469 * then this method can safely be ignored, it will be called
470 * transparently as needed.</p>
471
472 * <p><strong>Warning</strong>: since the step interpolator provided
473 * to the step handler as a parameter of the {@link
474 * StepHandler#handleStep handleStep} is valid only for the duration
475 * of the {@link StepHandler#handleStep handleStep} call, one cannot
476 * simply store a reference and reuse it later. One should first
477 * finalize the instance, then copy this finalized instance into a
478 * new object that can be kept.</p>
479
480 * <p>This method calls the protected <code>doFinalize</code> method
481 * if it has never been called during this step and set a flag
482 * indicating that it has been called once. It is the <code>
483 * doFinalize</code> method which should perform the evaluations.
484 * This wrapping prevents from calling <code>doFinalize</code> several
485 * times and hence evaluating the differential equations too often.
486 * Therefore, subclasses are not allowed not reimplement it, they
487 * should rather reimplement <code>doFinalize</code>.</p>
488
489 * @exception MaxCountExceededException if the number of functions evaluations is exceeded
490
491 */
492 public final void finalizeStep() throws MaxCountExceededException {
493 if (! finalized) {
494 doFinalize();
495 finalized = true;
496 }
497 }
498
499 /**
500 * Really finalize the step.
501 * The default implementation of this method does nothing.
502 * @exception MaxCountExceededException if the number of functions evaluations is exceeded
503 */
504 protected void doFinalize() throws MaxCountExceededException {
505 }
506
507 /** {@inheritDoc} */
508 @Override
509 public abstract void writeExternal(ObjectOutput out)
510 throws IOException;
511
512 /** {@inheritDoc} */
513 @Override
514 public abstract void readExternal(ObjectInput in)
515 throws IOException, ClassNotFoundException;
516
517 /** Save the base state of the instance.
518 * This method performs step finalization if it has not been done
519 * before.
520 * @param out stream where to save the state
521 * @exception IOException in case of write error
522 */
523 protected void writeBaseExternal(final ObjectOutput out)
524 throws IOException {
525
526 if (currentState == null) {
527 out.writeInt(-1);
528 } else {
529 out.writeInt(currentState.length);
530 }
531 out.writeDouble(globalPreviousTime);
532 out.writeDouble(globalCurrentTime);
533 out.writeDouble(softPreviousTime);
534 out.writeDouble(softCurrentTime);
535 out.writeDouble(h);
536 out.writeBoolean(forward);
537 out.writeObject(primaryMapper);
538 out.write(secondaryMappers.length);
539 for (final EquationsMapper mapper : secondaryMappers) {
540 out.writeObject(mapper);
541 }
542
543 if (currentState != null) {
544 for (int i = 0; i < currentState.length; ++i) {
545 out.writeDouble(currentState[i]);
546 }
547 }
548
549 out.writeDouble(interpolatedTime);
550
551 // we do not store the interpolated state,
552 // it will be recomputed as needed after reading
553
554 try {
555 // finalize the step (and don't bother saving the now true flag)
556 finalizeStep();
557 } catch (MaxCountExceededException mcee) {
558 final IOException ioe = new IOException(mcee.getLocalizedMessage());
559 ioe.initCause(mcee);
560 throw ioe;
561 }
562 }
563
564 /** Read the base state of the instance.
565 * This method does <strong>neither</strong> set the interpolated
566 * time nor state. It is up to the derived class to reset it
567 * properly calling the {@link #setInterpolatedTime} method later,
568 * once all rest of the object state has been set up properly.
569 * @param in stream where to read the state from
570 * @return interpolated time to be set later by the caller
571 * @exception IOException in case of read error
572 * @exception ClassNotFoundException if an equation mapper class
573 * cannot be found
574 */
575 protected double readBaseExternal(final ObjectInput in)
576 throws IOException, ClassNotFoundException {
577
578 final int dimension = in.readInt();
579 globalPreviousTime = in.readDouble();
580 globalCurrentTime = in.readDouble();
581 softPreviousTime = in.readDouble();
582 softCurrentTime = in.readDouble();
583 h = in.readDouble();
584 forward = in.readBoolean();
585 primaryMapper = (EquationsMapper) in.readObject();
586 secondaryMappers = new EquationsMapper[in.read()];
587 for (int i = 0; i < secondaryMappers.length; ++i) {
588 secondaryMappers[i] = (EquationsMapper) in.readObject();
589 }
590 dirtyState = true;
591
592 if (dimension < 0) {
593 currentState = null;
594 } else {
595 currentState = new double[dimension];
596 for (int i = 0; i < currentState.length; ++i) {
597 currentState[i] = in.readDouble();
598 }
599 }
600
601 // we do NOT handle the interpolated time and state here
602 interpolatedTime = Double.NaN;
603 allocateInterpolatedArrays(dimension);
604
605 finalized = true;
606
607 return in.readDouble();
608 }
609 }