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