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 }