Raymond | dee0849 | 2015-04-02 10:43:13 -0700 | [diff] [blame] | 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 | */ |
| 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 | } |