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 * @since 1.2
040 *
041 */
042
043public 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.math3.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() throws MaxCountExceededException {
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   * @exception MaxCountExceededException if the number of functions evaluations is exceeded
395   */
396  protected abstract void computeInterpolatedStateAndDerivatives(double theta,
397                                                                 double oneMinusThetaH)
398      throws MaxCountExceededException;
399
400  /** Lazy evaluation of complete interpolated state.
401   * @exception MaxCountExceededException if the number of functions evaluations is exceeded
402   */
403  private void evaluateCompleteInterpolatedState()
404      throws MaxCountExceededException {
405      // lazy evaluation of the state
406      if (dirtyState) {
407          final double oneMinusThetaH = globalCurrentTime - interpolatedTime;
408          final double theta = (h == 0) ? 0 : (h - oneMinusThetaH) / h;
409          computeInterpolatedStateAndDerivatives(theta, oneMinusThetaH);
410          dirtyState = false;
411      }
412  }
413
414  /** {@inheritDoc} */
415  public double[] getInterpolatedState() throws MaxCountExceededException {
416      evaluateCompleteInterpolatedState();
417      primaryMapper.extractEquationData(interpolatedState,
418                                        interpolatedPrimaryState);
419      return interpolatedPrimaryState;
420  }
421
422  /** {@inheritDoc} */
423  public double[] getInterpolatedDerivatives() throws MaxCountExceededException {
424      evaluateCompleteInterpolatedState();
425      primaryMapper.extractEquationData(interpolatedDerivatives,
426                                        interpolatedPrimaryDerivatives);
427      return interpolatedPrimaryDerivatives;
428  }
429
430  /** {@inheritDoc} */
431  public double[] getInterpolatedSecondaryState(final int index) throws MaxCountExceededException {
432      evaluateCompleteInterpolatedState();
433      secondaryMappers[index].extractEquationData(interpolatedState,
434                                                  interpolatedSecondaryState[index]);
435      return interpolatedSecondaryState[index];
436  }
437
438  /** {@inheritDoc} */
439  public double[] getInterpolatedSecondaryDerivatives(final int index) throws MaxCountExceededException {
440      evaluateCompleteInterpolatedState();
441      secondaryMappers[index].extractEquationData(interpolatedDerivatives,
442                                                  interpolatedSecondaryDerivatives[index]);
443      return interpolatedSecondaryDerivatives[index];
444  }
445
446  /**
447   * Finalize the step.
448
449   * <p>Some embedded Runge-Kutta integrators need fewer functions
450   * evaluations than their counterpart step interpolators. These
451   * interpolators should perform the last evaluations they need by
452   * themselves only if they need them. This method triggers these
453   * extra evaluations. It can be called directly by the user step
454   * handler and it is called automatically if {@link
455   * #setInterpolatedTime} is called.</p>
456
457   * <p>Once this method has been called, <strong>no</strong> other
458   * evaluation will be performed on this step. If there is a need to
459   * have some side effects between the step handler and the
460   * differential equations (for example update some data in the
461   * equations once the step has been done), it is advised to call
462   * this method explicitly from the step handler before these side
463   * effects are set up. If the step handler induces no side effect,
464   * then this method can safely be ignored, it will be called
465   * transparently as needed.</p>
466
467   * <p><strong>Warning</strong>: since the step interpolator provided
468   * to the step handler as a parameter of the {@link
469   * StepHandler#handleStep handleStep} is valid only for the duration
470   * of the {@link StepHandler#handleStep handleStep} call, one cannot
471   * simply store a reference and reuse it later. One should first
472   * finalize the instance, then copy this finalized instance into a
473   * new object that can be kept.</p>
474
475   * <p>This method calls the protected <code>doFinalize</code> method
476   * if it has never been called during this step and set a flag
477   * indicating that it has been called once. It is the <code>
478   * doFinalize</code> method which should perform the evaluations.
479   * This wrapping prevents from calling <code>doFinalize</code> several
480   * times and hence evaluating the differential equations too often.
481   * Therefore, subclasses are not allowed not reimplement it, they
482   * should rather reimplement <code>doFinalize</code>.</p>
483
484   * @exception MaxCountExceededException if the number of functions evaluations is exceeded
485
486   */
487  public final void finalizeStep() throws MaxCountExceededException {
488    if (! finalized) {
489      doFinalize();
490      finalized = true;
491    }
492  }
493
494  /**
495   * Really finalize the step.
496   * The default implementation of this method does nothing.
497   * @exception MaxCountExceededException if the number of functions evaluations is exceeded
498   */
499  protected void doFinalize() throws MaxCountExceededException {
500  }
501
502  /** {@inheritDoc} */
503  public abstract void writeExternal(ObjectOutput out)
504    throws IOException;
505
506  /** {@inheritDoc} */
507  public abstract void readExternal(ObjectInput in)
508    throws IOException, ClassNotFoundException;
509
510  /** Save the base state of the instance.
511   * This method performs step finalization if it has not been done
512   * before.
513   * @param out stream where to save the state
514   * @exception IOException in case of write error
515   */
516  protected void writeBaseExternal(final ObjectOutput out)
517    throws IOException {
518
519    if (currentState == null) {
520        out.writeInt(-1);
521    } else {
522        out.writeInt(currentState.length);
523    }
524    out.writeDouble(globalPreviousTime);
525    out.writeDouble(globalCurrentTime);
526    out.writeDouble(softPreviousTime);
527    out.writeDouble(softCurrentTime);
528    out.writeDouble(h);
529    out.writeBoolean(forward);
530    out.writeObject(primaryMapper);
531    out.write(secondaryMappers.length);
532    for (final EquationsMapper  mapper : secondaryMappers) {
533        out.writeObject(mapper);
534    }
535
536    if (currentState != null) {
537        for (int i = 0; i < currentState.length; ++i) {
538            out.writeDouble(currentState[i]);
539        }
540    }
541
542    out.writeDouble(interpolatedTime);
543
544    // we do not store the interpolated state,
545    // it will be recomputed as needed after reading
546
547    try {
548        // finalize the step (and don't bother saving the now true flag)
549        finalizeStep();
550    } catch (MaxCountExceededException mcee) {
551        final IOException ioe = new IOException(mcee.getLocalizedMessage());
552        ioe.initCause(mcee);
553        throw ioe;
554    }
555
556  }
557
558  /** Read the base state of the instance.
559   * This method does <strong>neither</strong> set the interpolated
560   * time nor state. It is up to the derived class to reset it
561   * properly calling the {@link #setInterpolatedTime} method later,
562   * once all rest of the object state has been set up properly.
563   * @param in stream where to read the state from
564   * @return interpolated time to be set later by the caller
565   * @exception IOException in case of read error
566   * @exception ClassNotFoundException if an equation mapper class
567   * cannot be found
568   */
569  protected double readBaseExternal(final ObjectInput in)
570    throws IOException, ClassNotFoundException {
571
572    final int dimension = in.readInt();
573    globalPreviousTime  = in.readDouble();
574    globalCurrentTime   = in.readDouble();
575    softPreviousTime    = in.readDouble();
576    softCurrentTime     = in.readDouble();
577    h                   = in.readDouble();
578    forward             = in.readBoolean();
579    primaryMapper       = (EquationsMapper) in.readObject();
580    secondaryMappers    = new EquationsMapper[in.read()];
581    for (int i = 0; i < secondaryMappers.length; ++i) {
582        secondaryMappers[i] = (EquationsMapper) in.readObject();
583    }
584    dirtyState          = true;
585
586    if (dimension < 0) {
587        currentState = null;
588    } else {
589        currentState  = new double[dimension];
590        for (int i = 0; i < currentState.length; ++i) {
591            currentState[i] = in.readDouble();
592        }
593    }
594
595    // we do NOT handle the interpolated time and state here
596    interpolatedTime = Double.NaN;
597    allocateInterpolatedArrays(dimension);
598
599    finalized = true;
600
601    return in.readDouble();
602
603  }
604
605}