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