001/* 002 * Licensed to the Apache Software Foundation (ASF) under one or more 003 * contributor license agreements. See the NOTICE file distributed with 004 * this work for additional information regarding copyright ownership. 005 * The ASF licenses this file to You under the Apache License, Version 2.0 006 * (the "License"); you may not use this file except in compliance with 007 * the License. You may obtain a copy of the License at 008 * 009 * http://www.apache.org/licenses/LICENSE-2.0 010 * 011 * Unless required by applicable law or agreed to in writing, software 012 * distributed under the License is distributed on an "AS IS" BASIS, 013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 014 * See the License for the specific language governing permissions and 015 * limitations under the License. 016 */ 017 018package org.apache.commons.math4.legacy.ode.sampling; 019 020import java.io.IOException; 021import java.io.ObjectInput; 022import java.io.ObjectOutput; 023 024import org.apache.commons.math4.legacy.exception.MaxCountExceededException; 025import org.apache.commons.math4.legacy.ode.EquationsMapper; 026 027/** This abstract class represents an interpolator over the last step 028 * during an ODE integration. 029 * 030 * <p>The various ODE integrators provide objects extending this class 031 * to the step handlers. The handlers can use these objects to 032 * retrieve the state vector at intermediate times between the 033 * previous and the current grid points (dense output).</p> 034 * 035 * @see org.apache.commons.math4.legacy.ode.FirstOrderIntegrator 036 * @see org.apache.commons.math4.legacy.ode.SecondOrderIntegrator 037 * @see StepHandler 038 * 039 * @since 1.2 040 * 041 */ 042 043public abstract class AbstractStepInterpolator 044 implements StepInterpolator { 045 046 /** current time step. */ 047 protected double h; 048 049 /** current state. */ 050 protected double[] currentState; 051 052 /** interpolated time. */ 053 protected double interpolatedTime; 054 055 /** interpolated state. */ 056 protected double[] interpolatedState; 057 058 /** interpolated derivatives. */ 059 protected double[] interpolatedDerivatives; 060 061 /** interpolated primary state. */ 062 protected double[] interpolatedPrimaryState; 063 064 /** interpolated primary derivatives. */ 065 protected double[] interpolatedPrimaryDerivatives; 066 067 /** interpolated secondary state. */ 068 protected double[][] interpolatedSecondaryState; 069 070 /** interpolated secondary derivatives. */ 071 protected double[][] interpolatedSecondaryDerivatives; 072 073 /** global previous time. */ 074 private double globalPreviousTime; 075 076 /** global current time. */ 077 private double globalCurrentTime; 078 079 /** soft previous time. */ 080 private double softPreviousTime; 081 082 /** soft current time. */ 083 private double softCurrentTime; 084 085 /** indicate if the step has been finalized or not. */ 086 private boolean finalized; 087 088 /** integration direction. */ 089 private boolean forward; 090 091 /** indicator for dirty state. */ 092 private boolean dirtyState; 093 094 /** Equations mapper for the primary equations set. */ 095 private EquationsMapper primaryMapper; 096 097 /** Equations mappers for the secondary equations sets. */ 098 private EquationsMapper[] secondaryMappers; 099 100 /** Simple constructor. 101 * This constructor builds an instance that is not usable yet, the 102 * {@link #reinitialize} method should be called before using the 103 * instance in order to initialize the internal arrays. This 104 * constructor is used only in order to delay the initialization in 105 * some cases. As an example, the {@link 106 * org.apache.commons.math4.legacy.ode.nonstiff.EmbeddedRungeKuttaIntegrator} 107 * class uses the prototyping design pattern to create the step 108 * interpolators by cloning an uninitialized model and latter 109 * initializing the copy. 110 */ 111 protected AbstractStepInterpolator() { 112 globalPreviousTime = Double.NaN; 113 globalCurrentTime = Double.NaN; 114 softPreviousTime = Double.NaN; 115 softCurrentTime = Double.NaN; 116 h = Double.NaN; 117 interpolatedTime = Double.NaN; 118 currentState = null; 119 finalized = false; 120 this.forward = true; 121 this.dirtyState = true; 122 primaryMapper = null; 123 secondaryMappers = null; 124 allocateInterpolatedArrays(-1); 125 } 126 127 /** Simple constructor. 128 * @param y reference to the integrator array holding the state at 129 * the end of the step 130 * @param forward integration direction indicator 131 * @param primaryMapper equations mapper for the primary equations set 132 * @param secondaryMappers equations mappers for the secondary equations sets 133 */ 134 protected AbstractStepInterpolator(final double[] y, final boolean forward, 135 final EquationsMapper primaryMapper, 136 final EquationsMapper[] secondaryMappers) { 137 138 globalPreviousTime = Double.NaN; 139 globalCurrentTime = Double.NaN; 140 softPreviousTime = Double.NaN; 141 softCurrentTime = Double.NaN; 142 h = Double.NaN; 143 interpolatedTime = Double.NaN; 144 currentState = y; 145 finalized = false; 146 this.forward = forward; 147 this.dirtyState = true; 148 this.primaryMapper = primaryMapper; 149 this.secondaryMappers = (secondaryMappers == null) ? null : secondaryMappers.clone(); 150 allocateInterpolatedArrays(y.length); 151 } 152 153 /** Copy constructor. 154 155 * <p>The copied interpolator should have been finalized before the 156 * copy, otherwise the copy will not be able to perform correctly 157 * any derivative computation and will throw a {@link 158 * NullPointerException} later. Since we don't want this constructor 159 * to throw the exceptions finalization may involve and since we 160 * don't want this method to modify the state of the copied 161 * interpolator, finalization is <strong>not</strong> done 162 * automatically, it remains under user control.</p> 163 164 * <p>The copy is a deep copy: its arrays are separated from the 165 * original arrays of the instance.</p> 166 167 * @param interpolator interpolator to copy from. 168 169 */ 170 protected AbstractStepInterpolator(final AbstractStepInterpolator interpolator) { 171 172 globalPreviousTime = interpolator.globalPreviousTime; 173 globalCurrentTime = interpolator.globalCurrentTime; 174 softPreviousTime = interpolator.softPreviousTime; 175 softCurrentTime = interpolator.softCurrentTime; 176 h = interpolator.h; 177 interpolatedTime = interpolator.interpolatedTime; 178 179 if (interpolator.currentState == null) { 180 currentState = null; 181 primaryMapper = null; 182 secondaryMappers = null; 183 allocateInterpolatedArrays(-1); 184 } else { 185 currentState = interpolator.currentState.clone(); 186 interpolatedState = interpolator.interpolatedState.clone(); 187 interpolatedDerivatives = interpolator.interpolatedDerivatives.clone(); 188 interpolatedPrimaryState = interpolator.interpolatedPrimaryState.clone(); 189 interpolatedPrimaryDerivatives = interpolator.interpolatedPrimaryDerivatives.clone(); 190 interpolatedSecondaryState = new double[interpolator.interpolatedSecondaryState.length][]; 191 interpolatedSecondaryDerivatives = new double[interpolator.interpolatedSecondaryDerivatives.length][]; 192 for (int i = 0; i < interpolatedSecondaryState.length; ++i) { 193 interpolatedSecondaryState[i] = interpolator.interpolatedSecondaryState[i].clone(); 194 interpolatedSecondaryDerivatives[i] = interpolator.interpolatedSecondaryDerivatives[i].clone(); 195 } 196 } 197 198 finalized = interpolator.finalized; 199 forward = interpolator.forward; 200 dirtyState = interpolator.dirtyState; 201 primaryMapper = interpolator.primaryMapper; 202 secondaryMappers = (interpolator.secondaryMappers == null) ? 203 null : interpolator.secondaryMappers.clone(); 204 } 205 206 /** Allocate the various interpolated states arrays. 207 * @param dimension total dimension (negative if arrays should be set to null) 208 */ 209 private void allocateInterpolatedArrays(final int dimension) { 210 if (dimension < 0) { 211 interpolatedState = null; 212 interpolatedDerivatives = null; 213 interpolatedPrimaryState = null; 214 interpolatedPrimaryDerivatives = null; 215 interpolatedSecondaryState = null; 216 interpolatedSecondaryDerivatives = null; 217 } else { 218 interpolatedState = new double[dimension]; 219 interpolatedDerivatives = new double[dimension]; 220 interpolatedPrimaryState = new double[primaryMapper.getDimension()]; 221 interpolatedPrimaryDerivatives = new double[primaryMapper.getDimension()]; 222 if (secondaryMappers == null) { 223 interpolatedSecondaryState = null; 224 interpolatedSecondaryDerivatives = null; 225 } else { 226 interpolatedSecondaryState = new double[secondaryMappers.length][]; 227 interpolatedSecondaryDerivatives = new double[secondaryMappers.length][]; 228 for (int i = 0; i < secondaryMappers.length; ++i) { 229 interpolatedSecondaryState[i] = new double[secondaryMappers[i].getDimension()]; 230 interpolatedSecondaryDerivatives[i] = new double[secondaryMappers[i].getDimension()]; 231 } 232 } 233 } 234 } 235 236 /** Reinitialize the instance. 237 * @param y reference to the integrator array holding the state at the end of the step 238 * @param isForward integration direction indicator 239 * @param primary equations mapper for the primary equations set 240 * @param secondary equations mappers for the secondary equations sets 241 */ 242 protected void reinitialize(final double[] y, final boolean isForward, 243 final EquationsMapper primary, 244 final EquationsMapper[] secondary) { 245 246 globalPreviousTime = Double.NaN; 247 globalCurrentTime = Double.NaN; 248 softPreviousTime = Double.NaN; 249 softCurrentTime = Double.NaN; 250 h = Double.NaN; 251 interpolatedTime = Double.NaN; 252 currentState = y; 253 finalized = false; 254 this.forward = isForward; 255 this.dirtyState = true; 256 this.primaryMapper = primary; 257 this.secondaryMappers = secondary.clone(); 258 allocateInterpolatedArrays(y.length); 259 } 260 261 /** {@inheritDoc} */ 262 @Override 263 public StepInterpolator copy() throws MaxCountExceededException { 264 265 // finalize the step before performing copy 266 finalizeStep(); 267 268 // create the new independent instance 269 return doCopy(); 270 } 271 272 /** Really copy the finalized instance. 273 * <p>This method is called by {@link #copy()} after the 274 * step has been finalized. It must perform a deep copy 275 * to have an new instance completely independent for the 276 * original instance. 277 * @return a copy of the finalized instance 278 */ 279 protected abstract StepInterpolator doCopy(); 280 281 /** Shift one step forward. 282 * Copy the current time into the previous time, hence preparing the 283 * interpolator for future calls to {@link #storeTime storeTime} 284 */ 285 public void shift() { 286 globalPreviousTime = globalCurrentTime; 287 softPreviousTime = globalPreviousTime; 288 softCurrentTime = globalCurrentTime; 289 } 290 291 /** Store the current step time. 292 * @param t current time 293 */ 294 public void storeTime(final double t) { 295 296 globalCurrentTime = t; 297 softCurrentTime = globalCurrentTime; 298 h = globalCurrentTime - globalPreviousTime; 299 setInterpolatedTime(t); 300 301 // the step is not finalized anymore 302 finalized = false; 303 } 304 305 /** Restrict step range to a limited part of the global step. 306 * <p> 307 * This method can be used to restrict a step and make it appear 308 * as if the original step was smaller. Calling this method 309 * <em>only</em> changes the value returned by {@link #getPreviousTime()}, 310 * it does not change any other property 311 * </p> 312 * @param softPreviousTime start of the restricted step 313 * @since 2.2 314 */ 315 public void setSoftPreviousTime(final double softPreviousTime) { 316 this.softPreviousTime = softPreviousTime; 317 } 318 319 /** Restrict step range to a limited part of the global step. 320 * <p> 321 * This method can be used to restrict a step and make it appear 322 * as if the original step was smaller. Calling this method 323 * <em>only</em> changes the value returned by {@link #getCurrentTime()}, 324 * it does not change any other property 325 * </p> 326 * @param softCurrentTime end of the restricted step 327 * @since 2.2 328 */ 329 public void setSoftCurrentTime(final double softCurrentTime) { 330 this.softCurrentTime = softCurrentTime; 331 } 332 333 /** 334 * Get the previous global grid point time. 335 * @return previous global grid point time 336 */ 337 public double getGlobalPreviousTime() { 338 return globalPreviousTime; 339 } 340 341 /** 342 * Get the current global grid point time. 343 * @return current global grid point time 344 */ 345 public double getGlobalCurrentTime() { 346 return globalCurrentTime; 347 } 348 349 /** 350 * Get the previous soft grid point time. 351 * @return previous soft grid point time 352 * @see #setSoftPreviousTime(double) 353 */ 354 @Override 355public double getPreviousTime() { 356 return softPreviousTime; 357 } 358 359 /** 360 * Get the current soft grid point time. 361 * @return current soft grid point time 362 * @see #setSoftCurrentTime(double) 363 */ 364 @Override 365public double getCurrentTime() { 366 return softCurrentTime; 367 } 368 369 /** {@inheritDoc} */ 370 @Override 371 public double getInterpolatedTime() { 372 return interpolatedTime; 373 } 374 375 /** {@inheritDoc} */ 376 @Override 377 public void setInterpolatedTime(final double time) { 378 interpolatedTime = time; 379 dirtyState = true; 380 } 381 382 /** {@inheritDoc} */ 383 @Override 384 public boolean isForward() { 385 return forward; 386 } 387 388 /** Compute the state and derivatives at the interpolated time. 389 * This is the main processing method that should be implemented by 390 * the derived classes to perform the interpolation. 391 * @param theta normalized interpolation abscissa within the step 392 * (theta is zero at the previous time step and one at the current time step) 393 * @param oneMinusThetaH time gap between the interpolated time and 394 * the current time 395 * @exception MaxCountExceededException if the number of functions evaluations is exceeded 396 */ 397 protected abstract void computeInterpolatedStateAndDerivatives(double theta, 398 double oneMinusThetaH) 399 throws MaxCountExceededException; 400 401 /** Lazy evaluation of complete interpolated state. 402 * @exception MaxCountExceededException if the number of functions evaluations is exceeded 403 */ 404 private void evaluateCompleteInterpolatedState() 405 throws MaxCountExceededException { 406 // lazy evaluation of the state 407 if (dirtyState) { 408 final double oneMinusThetaH = globalCurrentTime - interpolatedTime; 409 final double theta = (h == 0) ? 0 : (h - oneMinusThetaH) / h; 410 computeInterpolatedStateAndDerivatives(theta, oneMinusThetaH); 411 dirtyState = false; 412 } 413 } 414 415 /** {@inheritDoc} */ 416 @Override 417 public double[] getInterpolatedState() throws MaxCountExceededException { 418 evaluateCompleteInterpolatedState(); 419 primaryMapper.extractEquationData(interpolatedState, 420 interpolatedPrimaryState); 421 return interpolatedPrimaryState; 422 } 423 424 /** {@inheritDoc} */ 425 @Override 426 public double[] getInterpolatedDerivatives() throws MaxCountExceededException { 427 evaluateCompleteInterpolatedState(); 428 primaryMapper.extractEquationData(interpolatedDerivatives, 429 interpolatedPrimaryDerivatives); 430 return interpolatedPrimaryDerivatives; 431 } 432 433 /** {@inheritDoc} */ 434 @Override 435 public double[] getInterpolatedSecondaryState(final int index) throws MaxCountExceededException { 436 evaluateCompleteInterpolatedState(); 437 secondaryMappers[index].extractEquationData(interpolatedState, 438 interpolatedSecondaryState[index]); 439 return interpolatedSecondaryState[index]; 440 } 441 442 /** {@inheritDoc} */ 443 @Override 444 public double[] getInterpolatedSecondaryDerivatives(final int index) throws MaxCountExceededException { 445 evaluateCompleteInterpolatedState(); 446 secondaryMappers[index].extractEquationData(interpolatedDerivatives, 447 interpolatedSecondaryDerivatives[index]); 448 return interpolatedSecondaryDerivatives[index]; 449 } 450 451 /** 452 * Finalize the step. 453 454 * <p>Some embedded Runge-Kutta integrators need fewer functions 455 * evaluations than their counterpart step interpolators. These 456 * interpolators should perform the last evaluations they need by 457 * themselves only if they need them. This method triggers these 458 * extra evaluations. It can be called directly by the user step 459 * handler and it is called automatically if {@link 460 * #setInterpolatedTime} is called.</p> 461 462 * <p>Once this method has been called, <strong>no</strong> other 463 * evaluation will be performed on this step. If there is a need to 464 * have some side effects between the step handler and the 465 * differential equations (for example update some data in the 466 * equations once the step has been done), it is advised to call 467 * this method explicitly from the step handler before these side 468 * effects are set up. If the step handler induces no side effect, 469 * then this method can safely be ignored, it will be called 470 * transparently as needed.</p> 471 472 * <p><strong>Warning</strong>: since the step interpolator provided 473 * to the step handler as a parameter of the {@link 474 * StepHandler#handleStep handleStep} is valid only for the duration 475 * of the {@link StepHandler#handleStep handleStep} call, one cannot 476 * simply store a reference and reuse it later. One should first 477 * finalize the instance, then copy this finalized instance into a 478 * new object that can be kept.</p> 479 480 * <p>This method calls the protected <code>doFinalize</code> method 481 * if it has never been called during this step and set a flag 482 * indicating that it has been called once. It is the <code> 483 * doFinalize</code> method which should perform the evaluations. 484 * This wrapping prevents from calling <code>doFinalize</code> several 485 * times and hence evaluating the differential equations too often. 486 * Therefore, subclasses are not allowed not reimplement it, they 487 * should rather reimplement <code>doFinalize</code>.</p> 488 489 * @exception MaxCountExceededException if the number of functions evaluations is exceeded 490 491 */ 492 public final void finalizeStep() throws MaxCountExceededException { 493 if (! finalized) { 494 doFinalize(); 495 finalized = true; 496 } 497 } 498 499 /** 500 * Really finalize the step. 501 * The default implementation of this method does nothing. 502 * @exception MaxCountExceededException if the number of functions evaluations is exceeded 503 */ 504 protected void doFinalize() throws MaxCountExceededException { 505 } 506 507 /** {@inheritDoc} */ 508 @Override 509 public abstract void writeExternal(ObjectOutput out) 510 throws IOException; 511 512 /** {@inheritDoc} */ 513 @Override 514 public abstract void readExternal(ObjectInput in) 515 throws IOException, ClassNotFoundException; 516 517 /** Save the base state of the instance. 518 * This method performs step finalization if it has not been done 519 * before. 520 * @param out stream where to save the state 521 * @exception IOException in case of write error 522 */ 523 protected void writeBaseExternal(final ObjectOutput out) 524 throws IOException { 525 526 if (currentState == null) { 527 out.writeInt(-1); 528 } else { 529 out.writeInt(currentState.length); 530 } 531 out.writeDouble(globalPreviousTime); 532 out.writeDouble(globalCurrentTime); 533 out.writeDouble(softPreviousTime); 534 out.writeDouble(softCurrentTime); 535 out.writeDouble(h); 536 out.writeBoolean(forward); 537 out.writeObject(primaryMapper); 538 out.write(secondaryMappers.length); 539 for (final EquationsMapper mapper : secondaryMappers) { 540 out.writeObject(mapper); 541 } 542 543 if (currentState != null) { 544 for (int i = 0; i < currentState.length; ++i) { 545 out.writeDouble(currentState[i]); 546 } 547 } 548 549 out.writeDouble(interpolatedTime); 550 551 // we do not store the interpolated state, 552 // it will be recomputed as needed after reading 553 554 try { 555 // finalize the step (and don't bother saving the now true flag) 556 finalizeStep(); 557 } catch (MaxCountExceededException mcee) { 558 final IOException ioe = new IOException(mcee.getLocalizedMessage()); 559 ioe.initCause(mcee); 560 throw ioe; 561 } 562 } 563 564 /** Read the base state of the instance. 565 * This method does <strong>neither</strong> set the interpolated 566 * time nor state. It is up to the derived class to reset it 567 * properly calling the {@link #setInterpolatedTime} method later, 568 * once all rest of the object state has been set up properly. 569 * @param in stream where to read the state from 570 * @return interpolated time to be set later by the caller 571 * @exception IOException in case of read error 572 * @exception ClassNotFoundException if an equation mapper class 573 * cannot be found 574 */ 575 protected double readBaseExternal(final ObjectInput in) 576 throws IOException, ClassNotFoundException { 577 578 final int dimension = in.readInt(); 579 globalPreviousTime = in.readDouble(); 580 globalCurrentTime = in.readDouble(); 581 softPreviousTime = in.readDouble(); 582 softCurrentTime = in.readDouble(); 583 h = in.readDouble(); 584 forward = in.readBoolean(); 585 primaryMapper = (EquationsMapper) in.readObject(); 586 secondaryMappers = new EquationsMapper[in.read()]; 587 for (int i = 0; i < secondaryMappers.length; ++i) { 588 secondaryMappers[i] = (EquationsMapper) in.readObject(); 589 } 590 dirtyState = true; 591 592 if (dimension < 0) { 593 currentState = null; 594 } else { 595 currentState = new double[dimension]; 596 for (int i = 0; i < currentState.length; ++i) { 597 currentState[i] = in.readDouble(); 598 } 599 } 600 601 // we do NOT handle the interpolated time and state here 602 interpolatedTime = Double.NaN; 603 allocateInterpolatedArrays(dimension); 604 605 finalized = true; 606 607 return in.readDouble(); 608 } 609}