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 fvr. 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 */ 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 */ 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 */ 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 */ 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} */ 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 */ 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 */ 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 */ 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 */ 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 */ 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 */ 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 */ 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 */ 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 */ 304 public double getCurrentTime() { 305 return softCurrentTime; 306 } 307 308 /** {@inheritDoc} */ 309 public double getInterpolatedTime() { 310 return interpolatedTime; 311 } 312 313 /** {@inheritDoc} */ 314 public void setInterpolatedTime(final double time) { 315 interpolatedTime = time; 316 dirtyState = true; 317 } 318 319 /** {@inheritDoc} */ 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 */ 334 protected abstract void computeInterpolatedStateAndDerivatives(double theta, 335 double oneMinusThetaH) 336 throws DerivativeException; 337 338 /** {@inheritDoc} */ 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} */ 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 */ 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 */ 423 protected void doFinalize() 424 throws DerivativeException { 425 } 426 427 /** {@inheritDoc} */ 428 public abstract void writeExternal(ObjectOutput out) 429 throws IOException; 430 431 /** {@inheritDoc} */ 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 */ 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 */ 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
