• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
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.math.ode.sampling;
19 
20 import java.io.IOException;
21 import java.io.ObjectInput;
22 import java.io.ObjectOutput;
23 
24 import org.apache.commons.math.ode.DerivativeException;
25 
26 /** This abstract class represents an interpolator over the last step
27  * during an ODE integration.
28  *
29  * <p>The various ODE integrators provide objects extending this class
30  * to the step handlers. The handlers can use these objects to
31  * retrieve the state vector at intermediate times between the
32  * previous and the current grid points (dense output).</p>
33  *
34  * @see org.apache.commons.math.ode.FirstOrderIntegrator
35  * @see org.apache.commons.math.ode.SecondOrderIntegrator
36  * @see StepHandler
37  *
38  * @version $Revision: 1073158 $ $Date: 2011-02-21 22:46:52 +0100 (lun. 21 févr. 2011) $
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   /** global previous time */
62   private double globalPreviousTime;
63 
64   /** global current time */
65   private double globalCurrentTime;
66 
67   /** soft previous time */
68   private double softPreviousTime;
69 
70   /** soft current time */
71   private double softCurrentTime;
72 
73   /** indicate if the step has been finalized or not. */
74   private boolean finalized;
75 
76   /** integration direction. */
77   private boolean forward;
78 
79   /** indicator for dirty state. */
80   private boolean dirtyState;
81 
82 
83   /** Simple constructor.
84    * This constructor builds an instance that is not usable yet, the
85    * {@link #reinitialize} method should be called before using the
86    * instance in order to initialize the internal arrays. This
87    * constructor is used only in order to delay the initialization in
88    * some cases. As an example, the {@link
89    * org.apache.commons.math.ode.nonstiff.EmbeddedRungeKuttaIntegrator}
90    * class uses the prototyping design pattern to create the step
91    * interpolators by cloning an uninitialized model and latter
92    * initializing the copy.
93    */
AbstractStepInterpolator()94   protected AbstractStepInterpolator() {
95     globalPreviousTime      = Double.NaN;
96     globalCurrentTime       = Double.NaN;
97     softPreviousTime        = Double.NaN;
98     softCurrentTime         = Double.NaN;
99     h                       = Double.NaN;
100     interpolatedTime        = Double.NaN;
101     currentState            = null;
102     interpolatedState       = null;
103     interpolatedDerivatives = null;
104     finalized               = false;
105     this.forward            = true;
106     this.dirtyState         = true;
107   }
108 
109   /** Simple constructor.
110    * @param y reference to the integrator array holding the state at
111    * the end of the step
112    * @param forward integration direction indicator
113    */
AbstractStepInterpolator(final double[] y, final boolean forward)114   protected AbstractStepInterpolator(final double[] y, final boolean forward) {
115 
116     globalPreviousTime = Double.NaN;
117     globalCurrentTime  = Double.NaN;
118     softPreviousTime   = Double.NaN;
119     softCurrentTime    = Double.NaN;
120     h                  = Double.NaN;
121     interpolatedTime   = Double.NaN;
122 
123     currentState            = y;
124     interpolatedState       = new double[y.length];
125     interpolatedDerivatives = new double[y.length];
126 
127     finalized         = false;
128     this.forward      = forward;
129     this.dirtyState   = true;
130 
131   }
132 
133   /** Copy constructor.
134 
135    * <p>The copied interpolator should have been finalized before the
136    * copy, otherwise the copy will not be able to perform correctly
137    * any derivative computation and will throw a {@link
138    * NullPointerException} later. Since we don't want this constructor
139    * to throw the exceptions finalization may involve and since we
140    * don't want this method to modify the state of the copied
141    * interpolator, finalization is <strong>not</strong> done
142    * automatically, it remains under user control.</p>
143    *
144    * <p>The copy is a deep copy: its arrays are separated from the
145    * original arrays of the instance.</p>
146    *
147    * @param interpolator interpolator to copy from.
148    *
149    */
AbstractStepInterpolator(final AbstractStepInterpolator interpolator)150   protected AbstractStepInterpolator(final AbstractStepInterpolator interpolator) {
151 
152     globalPreviousTime = interpolator.globalPreviousTime;
153     globalCurrentTime  = interpolator.globalCurrentTime;
154     softPreviousTime   = interpolator.softPreviousTime;
155     softCurrentTime    = interpolator.softCurrentTime;
156     h                  = interpolator.h;
157     interpolatedTime   = interpolator.interpolatedTime;
158 
159     if (interpolator.currentState != null) {
160       currentState            = interpolator.currentState.clone();
161       interpolatedState       = interpolator.interpolatedState.clone();
162       interpolatedDerivatives = interpolator.interpolatedDerivatives.clone();
163     } else {
164       currentState            = null;
165       interpolatedState       = null;
166       interpolatedDerivatives = null;
167     }
168 
169     finalized  = interpolator.finalized;
170     forward    = interpolator.forward;
171     dirtyState = interpolator.dirtyState;
172 
173   }
174 
175   /** Reinitialize the instance
176    * @param y reference to the integrator array holding the state at
177    * the end of the step
178    * @param isForward integration direction indicator
179    */
reinitialize(final double[] y, final boolean isForward)180   protected void reinitialize(final double[] y, final boolean isForward) {
181 
182     globalPreviousTime = Double.NaN;
183     globalCurrentTime  = Double.NaN;
184     softPreviousTime   = Double.NaN;
185     softCurrentTime    = Double.NaN;
186     h                  = Double.NaN;
187     interpolatedTime   = Double.NaN;
188 
189     currentState            = y;
190     interpolatedState       = new double[y.length];
191     interpolatedDerivatives = new double[y.length];
192 
193     finalized         = false;
194     this.forward      = isForward;
195     this.dirtyState   = true;
196 
197   }
198 
199   /** {@inheritDoc} */
copy()200    public StepInterpolator copy() throws DerivativeException {
201 
202      // finalize the step before performing copy
203      finalizeStep();
204 
205      // create the new independent instance
206      return doCopy();
207 
208    }
209 
210    /** Really copy the finalized instance.
211     * <p>This method is called by {@link #copy()} after the
212     * step has been finalized. It must perform a deep copy
213     * to have an new instance completely independent for the
214     * original instance.
215     * @return a copy of the finalized instance
216     */
doCopy()217    protected abstract StepInterpolator doCopy();
218 
219   /** Shift one step forward.
220    * Copy the current time into the previous time, hence preparing the
221    * interpolator for future calls to {@link #storeTime storeTime}
222    */
shift()223   public void shift() {
224     globalPreviousTime = globalCurrentTime;
225     softPreviousTime   = globalPreviousTime;
226     softCurrentTime    = globalCurrentTime;
227   }
228 
229   /** Store the current step time.
230    * @param t current time
231    */
storeTime(final double t)232   public void storeTime(final double t) {
233 
234     globalCurrentTime = t;
235     softCurrentTime   = globalCurrentTime;
236     h                 = globalCurrentTime - globalPreviousTime;
237     setInterpolatedTime(t);
238 
239     // the step is not finalized anymore
240     finalized  = false;
241 
242   }
243 
244   /** Restrict step range to a limited part of the global step.
245    * <p>
246    * This method can be used to restrict a step and make it appear
247    * as if the original step was smaller. Calling this method
248    * <em>only</em> changes the value returned by {@link #getPreviousTime()},
249    * it does not change any other property
250    * </p>
251    * @param softPreviousTime start of the restricted step
252    * @since 2.2
253    */
setSoftPreviousTime(final double softPreviousTime)254   public void setSoftPreviousTime(final double softPreviousTime) {
255       this.softPreviousTime = softPreviousTime;
256   }
257 
258   /** Restrict step range to a limited part of the global step.
259    * <p>
260    * This method can be used to restrict a step and make it appear
261    * as if the original step was smaller. Calling this method
262    * <em>only</em> changes the value returned by {@link #getCurrentTime()},
263    * it does not change any other property
264    * </p>
265    * @param softCurrentTime end of the restricted step
266    * @since 2.2
267    */
setSoftCurrentTime(final double softCurrentTime)268   public void setSoftCurrentTime(final double softCurrentTime) {
269       this.softCurrentTime  = softCurrentTime;
270   }
271 
272   /**
273    * Get the previous global grid point time.
274    * @return previous global grid point time
275    * @since 2.2
276    */
getGlobalPreviousTime()277   public double getGlobalPreviousTime() {
278     return globalPreviousTime;
279   }
280 
281   /**
282    * Get the current global grid point time.
283    * @return current global grid point time
284    * @since 2.2
285    */
getGlobalCurrentTime()286   public double getGlobalCurrentTime() {
287     return globalCurrentTime;
288   }
289 
290   /**
291    * Get the previous soft grid point time.
292    * @return previous soft grid point time
293    * @see #setSoftPreviousTime(double)
294    */
getPreviousTime()295   public double getPreviousTime() {
296     return softPreviousTime;
297   }
298 
299   /**
300    * Get the current soft grid point time.
301    * @return current soft grid point time
302    * @see #setSoftCurrentTime(double)
303    */
getCurrentTime()304   public double getCurrentTime() {
305     return softCurrentTime;
306   }
307 
308   /** {@inheritDoc} */
getInterpolatedTime()309   public double getInterpolatedTime() {
310     return interpolatedTime;
311   }
312 
313   /** {@inheritDoc} */
setInterpolatedTime(final double time)314   public void setInterpolatedTime(final double time) {
315       interpolatedTime = time;
316       dirtyState       = true;
317   }
318 
319   /** {@inheritDoc} */
isForward()320   public boolean isForward() {
321     return forward;
322   }
323 
324   /** Compute the state and derivatives at the interpolated time.
325    * This is the main processing method that should be implemented by
326    * the derived classes to perform the interpolation.
327    * @param theta normalized interpolation abscissa within the step
328    * (theta is zero at the previous time step and one at the current time step)
329    * @param oneMinusThetaH time gap between the interpolated time and
330    * the current time
331    * @throws DerivativeException this exception is propagated to the caller if the
332    * underlying user function triggers one
333    */
computeInterpolatedStateAndDerivatives(double theta, double oneMinusThetaH)334   protected abstract void computeInterpolatedStateAndDerivatives(double theta,
335                                                                  double oneMinusThetaH)
336     throws DerivativeException;
337 
338   /** {@inheritDoc} */
getInterpolatedState()339   public double[] getInterpolatedState() throws DerivativeException {
340 
341       // lazy evaluation of the state
342       if (dirtyState) {
343           final double oneMinusThetaH = globalCurrentTime - interpolatedTime;
344           final double theta = (h == 0) ? 0 : (h - oneMinusThetaH) / h;
345           computeInterpolatedStateAndDerivatives(theta, oneMinusThetaH);
346           dirtyState = false;
347       }
348 
349       return interpolatedState;
350 
351   }
352 
353   /** {@inheritDoc} */
getInterpolatedDerivatives()354   public double[] getInterpolatedDerivatives() throws DerivativeException {
355 
356       // lazy evaluation of the state
357       if (dirtyState) {
358           final double oneMinusThetaH = globalCurrentTime - interpolatedTime;
359           final double theta = (h == 0) ? 0 : (h - oneMinusThetaH) / h;
360           computeInterpolatedStateAndDerivatives(theta, oneMinusThetaH);
361           dirtyState = false;
362       }
363 
364       return interpolatedDerivatives;
365 
366   }
367 
368   /**
369    * Finalize the step.
370    *
371    * <p>Some embedded Runge-Kutta integrators need fewer functions
372    * evaluations than their counterpart step interpolators. These
373    * interpolators should perform the last evaluations they need by
374    * themselves only if they need them. This method triggers these
375    * extra evaluations. It can be called directly by the user step
376    * handler and it is called automatically if {@link
377    * #setInterpolatedTime} is called.</p>
378    *
379    * <p>Once this method has been called, <strong>no</strong> other
380    * evaluation will be performed on this step. If there is a need to
381    * have some side effects between the step handler and the
382    * differential equations (for example update some data in the
383    * equations once the step has been done), it is advised to call
384    * this method explicitly from the step handler before these side
385    * effects are set up. If the step handler induces no side effect,
386    * then this method can safely be ignored, it will be called
387    * transparently as needed.</p>
388    *
389    * <p><strong>Warning</strong>: since the step interpolator provided
390    * to the step handler as a parameter of the {@link
391    * StepHandler#handleStep handleStep} is valid only for the duration
392    * of the {@link StepHandler#handleStep handleStep} call, one cannot
393    * simply store a reference and reuse it later. One should first
394    * finalize the instance, then copy this finalized instance into a
395    * new object that can be kept.</p>
396    *
397    * <p>This method calls the protected <code>doFinalize</code> method
398    * if it has never been called during this step and set a flag
399    * indicating that it has been called once. It is the <code>
400    * doFinalize</code> method which should perform the evaluations.
401    * This wrapping prevents from calling <code>doFinalize</code> several
402    * times and hence evaluating the differential equations too often.
403    * Therefore, subclasses are not allowed not reimplement it, they
404    * should rather reimplement <code>doFinalize</code>.</p>
405    *
406    * @throws DerivativeException this exception is propagated to the
407    * caller if the underlying user function triggers one
408    */
finalizeStep()409   public final void finalizeStep()
410     throws DerivativeException {
411     if (! finalized) {
412       doFinalize();
413       finalized = true;
414     }
415   }
416 
417   /**
418    * Really finalize the step.
419    * The default implementation of this method does nothing.
420    * @throws DerivativeException this exception is propagated to the
421    * caller if the underlying user function triggers one
422    */
doFinalize()423   protected void doFinalize()
424     throws DerivativeException {
425   }
426 
427   /** {@inheritDoc} */
writeExternal(ObjectOutput out)428   public abstract void writeExternal(ObjectOutput out)
429     throws IOException;
430 
431   /** {@inheritDoc} */
readExternal(ObjectInput in)432   public abstract void readExternal(ObjectInput in)
433     throws IOException, ClassNotFoundException;
434 
435   /** Save the base state of the instance.
436    * This method performs step finalization if it has not been done
437    * before.
438    * @param out stream where to save the state
439    * @exception IOException in case of write error
440    */
writeBaseExternal(final ObjectOutput out)441   protected void writeBaseExternal(final ObjectOutput out)
442     throws IOException {
443 
444     if (currentState == null) {
445         out.writeInt(-1);
446     } else {
447         out.writeInt(currentState.length);
448     }
449     out.writeDouble(globalPreviousTime);
450     out.writeDouble(globalCurrentTime);
451     out.writeDouble(softPreviousTime);
452     out.writeDouble(softCurrentTime);
453     out.writeDouble(h);
454     out.writeBoolean(forward);
455 
456     if (currentState != null) {
457         for (int i = 0; i < currentState.length; ++i) {
458             out.writeDouble(currentState[i]);
459         }
460     }
461 
462     out.writeDouble(interpolatedTime);
463 
464     // we do not store the interpolated state,
465     // it will be recomputed as needed after reading
466 
467     // finalize the step (and don't bother saving the now true flag)
468     try {
469       finalizeStep();
470     } catch (DerivativeException e) {
471         IOException ioe = new IOException(e.getLocalizedMessage());
472         ioe.initCause(e);
473         throw ioe;
474     }
475 
476   }
477 
478   /** Read the base state of the instance.
479    * This method does <strong>neither</strong> set the interpolated
480    * time nor state. It is up to the derived class to reset it
481    * properly calling the {@link #setInterpolatedTime} method later,
482    * once all rest of the object state has been set up properly.
483    * @param in stream where to read the state from
484    * @return interpolated time be set later by the caller
485    * @exception IOException in case of read error
486    */
readBaseExternal(final ObjectInput in)487   protected double readBaseExternal(final ObjectInput in)
488     throws IOException {
489 
490     final int dimension = in.readInt();
491     globalPreviousTime  = in.readDouble();
492     globalCurrentTime   = in.readDouble();
493     softPreviousTime    = in.readDouble();
494     softCurrentTime     = in.readDouble();
495     h                   = in.readDouble();
496     forward             = in.readBoolean();
497     dirtyState          = true;
498 
499     if (dimension < 0) {
500         currentState = null;
501     } else {
502         currentState  = new double[dimension];
503         for (int i = 0; i < currentState.length; ++i) {
504             currentState[i] = in.readDouble();
505         }
506     }
507 
508     // we do NOT handle the interpolated time and state here
509     interpolatedTime        = Double.NaN;
510     interpolatedState       = (dimension < 0) ? null : new double[dimension];
511     interpolatedDerivatives = (dimension < 0) ? null : new double[dimension];
512 
513     finalized = true;
514 
515     return in.readDouble();
516 
517   }
518 
519 }
520