View Javadoc

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