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
018package org.apache.commons.math3.ode.sampling;
019
020import java.io.IOException;
021import java.io.ObjectInput;
022import java.io.ObjectOutput;
023
024import org.apache.commons.math3.exception.MaxCountExceededException;
025import 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
044public 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}