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.math3.ode; 019 020import java.io.Serializable; 021import java.util.ArrayList; 022import java.util.List; 023 024import org.apache.commons.math3.exception.DimensionMismatchException; 025import org.apache.commons.math3.exception.MathIllegalArgumentException; 026import org.apache.commons.math3.exception.MaxCountExceededException; 027import org.apache.commons.math3.exception.util.LocalizedFormats; 028import org.apache.commons.math3.ode.sampling.StepHandler; 029import org.apache.commons.math3.ode.sampling.StepInterpolator; 030import org.apache.commons.math3.util.FastMath; 031 032/** 033 * This class stores all information provided by an ODE integrator 034 * during the integration process and build a continuous model of the 035 * solution from this. 036 * 037 * <p>This class act as a step handler from the integrator point of 038 * view. It is called iteratively during the integration process and 039 * stores a copy of all steps information in a sorted collection for 040 * later use. Once the integration process is over, the user can use 041 * the {@link #setInterpolatedTime setInterpolatedTime} and {@link 042 * #getInterpolatedState getInterpolatedState} to retrieve this 043 * information at any time. It is important to wait for the 044 * integration to be over before attempting to call {@link 045 * #setInterpolatedTime setInterpolatedTime} because some internal 046 * variables are set only once the last step has been handled.</p> 047 * 048 * <p>This is useful for example if the main loop of the user 049 * application should remain independent from the integration process 050 * or if one needs to mimic the behaviour of an analytical model 051 * despite a numerical model is used (i.e. one needs the ability to 052 * get the model value at any time or to navigate through the 053 * data).</p> 054 * 055 * <p>If problem modeling is done with several separate 056 * integration phases for contiguous intervals, the same 057 * ContinuousOutputModel can be used as step handler for all 058 * integration phases as long as they are performed in order and in 059 * the same direction. As an example, one can extrapolate the 060 * trajectory of a satellite with one model (i.e. one set of 061 * differential equations) up to the beginning of a maneuver, use 062 * another more complex model including thrusters modeling and 063 * accurate attitude control during the maneuver, and revert to the 064 * first model after the end of the maneuver. If the same continuous 065 * output model handles the steps of all integration phases, the user 066 * do not need to bother when the maneuver begins or ends, he has all 067 * the data available in a transparent manner.</p> 068 * 069 * <p>An important feature of this class is that it implements the 070 * <code>Serializable</code> interface. This means that the result of 071 * an integration can be serialized and reused later (if stored into a 072 * persistent medium like a filesystem or a database) or elsewhere (if 073 * sent to another application). Only the result of the integration is 074 * stored, there is no reference to the integrated problem by 075 * itself.</p> 076 * 077 * <p>One should be aware that the amount of data stored in a 078 * ContinuousOutputModel instance can be important if the state vector 079 * is large, if the integration interval is long or if the steps are 080 * small (which can result from small tolerance settings in {@link 081 * org.apache.commons.math3.ode.nonstiff.AdaptiveStepsizeIntegrator adaptive 082 * step size integrators}).</p> 083 * 084 * @see StepHandler 085 * @see StepInterpolator 086 * @since 1.2 087 */ 088 089public class ContinuousOutputModel 090 implements StepHandler, Serializable { 091 092 /** Serializable version identifier */ 093 private static final long serialVersionUID = -1417964919405031606L; 094 095 /** Initial integration time. */ 096 private double initialTime; 097 098 /** Final integration time. */ 099 private double finalTime; 100 101 /** Integration direction indicator. */ 102 private boolean forward; 103 104 /** Current interpolator index. */ 105 private int index; 106 107 /** Steps table. */ 108 private List<StepInterpolator> steps; 109 110 /** Simple constructor. 111 * Build an empty continuous output model. 112 */ 113 public ContinuousOutputModel() { 114 steps = new ArrayList<StepInterpolator>(); 115 initialTime = Double.NaN; 116 finalTime = Double.NaN; 117 forward = true; 118 index = 0; 119 } 120 121 /** Append another model at the end of the instance. 122 * @param model model to add at the end of the instance 123 * @exception MathIllegalArgumentException if the model to append is not 124 * compatible with the instance (dimension of the state vector, 125 * propagation direction, hole between the dates) 126 * @exception MaxCountExceededException if the number of functions evaluations is exceeded 127 * during step finalization 128 */ 129 public void append(final ContinuousOutputModel model) 130 throws MathIllegalArgumentException, MaxCountExceededException { 131 132 if (model.steps.size() == 0) { 133 return; 134 } 135 136 if (steps.size() == 0) { 137 initialTime = model.initialTime; 138 forward = model.forward; 139 } else { 140 141 if (getInterpolatedState().length != model.getInterpolatedState().length) { 142 throw new DimensionMismatchException(model.getInterpolatedState().length, 143 getInterpolatedState().length); 144 } 145 146 if (forward ^ model.forward) { 147 throw new MathIllegalArgumentException(LocalizedFormats.PROPAGATION_DIRECTION_MISMATCH); 148 } 149 150 final StepInterpolator lastInterpolator = steps.get(index); 151 final double current = lastInterpolator.getCurrentTime(); 152 final double previous = lastInterpolator.getPreviousTime(); 153 final double step = current - previous; 154 final double gap = model.getInitialTime() - current; 155 if (FastMath.abs(gap) > 1.0e-3 * FastMath.abs(step)) { 156 throw new MathIllegalArgumentException(LocalizedFormats.HOLE_BETWEEN_MODELS_TIME_RANGES, 157 FastMath.abs(gap)); 158 } 159 160 } 161 162 for (StepInterpolator interpolator : model.steps) { 163 steps.add(interpolator.copy()); 164 } 165 166 index = steps.size() - 1; 167 finalTime = (steps.get(index)).getCurrentTime(); 168 169 } 170 171 /** {@inheritDoc} */ 172 public void init(double t0, double[] y0, double t) { 173 initialTime = Double.NaN; 174 finalTime = Double.NaN; 175 forward = true; 176 index = 0; 177 steps.clear(); 178 } 179 180 /** Handle the last accepted step. 181 * A copy of the information provided by the last step is stored in 182 * the instance for later use. 183 * @param interpolator interpolator for the last accepted step. 184 * @param isLast true if the step is the last one 185 * @exception MaxCountExceededException if the number of functions evaluations is exceeded 186 * during step finalization 187 */ 188 public void handleStep(final StepInterpolator interpolator, final boolean isLast) 189 throws MaxCountExceededException { 190 191 if (steps.size() == 0) { 192 initialTime = interpolator.getPreviousTime(); 193 forward = interpolator.isForward(); 194 } 195 196 steps.add(interpolator.copy()); 197 198 if (isLast) { 199 finalTime = interpolator.getCurrentTime(); 200 index = steps.size() - 1; 201 } 202 203 } 204 205 /** 206 * Get the initial integration time. 207 * @return initial integration time 208 */ 209 public double getInitialTime() { 210 return initialTime; 211 } 212 213 /** 214 * Get the final integration time. 215 * @return final integration time 216 */ 217 public double getFinalTime() { 218 return finalTime; 219 } 220 221 /** 222 * Get the time of the interpolated point. 223 * If {@link #setInterpolatedTime} has not been called, it returns 224 * the final integration time. 225 * @return interpolation point time 226 */ 227 public double getInterpolatedTime() { 228 return steps.get(index).getInterpolatedTime(); 229 } 230 231 /** Set the time of the interpolated point. 232 * <p>This method should <strong>not</strong> be called before the 233 * integration is over because some internal variables are set only 234 * once the last step has been handled.</p> 235 * <p>Setting the time outside of the integration interval is now 236 * allowed, but should be used with care since the accuracy of the 237 * interpolator will probably be very poor far from this interval. 238 * This allowance has been added to simplify implementation of search 239 * algorithms near the interval endpoints.</p> 240 * <p>Note that each time this method is called, the internal arrays 241 * returned in {@link #getInterpolatedState()}, {@link 242 * #getInterpolatedDerivatives()} and {@link #getInterpolatedSecondaryState(int)} 243 * <em>will</em> be overwritten. So if their content must be preserved 244 * across several calls, user must copy them.</p> 245 * @param time time of the interpolated point 246 * @see #getInterpolatedState() 247 * @see #getInterpolatedDerivatives() 248 * @see #getInterpolatedSecondaryState(int) 249 */ 250 public void setInterpolatedTime(final double time) { 251 252 // initialize the search with the complete steps table 253 int iMin = 0; 254 final StepInterpolator sMin = steps.get(iMin); 255 double tMin = 0.5 * (sMin.getPreviousTime() + sMin.getCurrentTime()); 256 257 int iMax = steps.size() - 1; 258 final StepInterpolator sMax = steps.get(iMax); 259 double tMax = 0.5 * (sMax.getPreviousTime() + sMax.getCurrentTime()); 260 261 // handle points outside of the integration interval 262 // or in the first and last step 263 if (locatePoint(time, sMin) <= 0) { 264 index = iMin; 265 sMin.setInterpolatedTime(time); 266 return; 267 } 268 if (locatePoint(time, sMax) >= 0) { 269 index = iMax; 270 sMax.setInterpolatedTime(time); 271 return; 272 } 273 274 // reduction of the table slice size 275 while (iMax - iMin > 5) { 276 277 // use the last estimated index as the splitting index 278 final StepInterpolator si = steps.get(index); 279 final int location = locatePoint(time, si); 280 if (location < 0) { 281 iMax = index; 282 tMax = 0.5 * (si.getPreviousTime() + si.getCurrentTime()); 283 } else if (location > 0) { 284 iMin = index; 285 tMin = 0.5 * (si.getPreviousTime() + si.getCurrentTime()); 286 } else { 287 // we have found the target step, no need to continue searching 288 si.setInterpolatedTime(time); 289 return; 290 } 291 292 // compute a new estimate of the index in the reduced table slice 293 final int iMed = (iMin + iMax) / 2; 294 final StepInterpolator sMed = steps.get(iMed); 295 final double tMed = 0.5 * (sMed.getPreviousTime() + sMed.getCurrentTime()); 296 297 if ((FastMath.abs(tMed - tMin) < 1e-6) || (FastMath.abs(tMax - tMed) < 1e-6)) { 298 // too close to the bounds, we estimate using a simple dichotomy 299 index = iMed; 300 } else { 301 // estimate the index using a reverse quadratic polynom 302 // (reverse means we have i = P(t), thus allowing to simply 303 // compute index = P(time) rather than solving a quadratic equation) 304 final double d12 = tMax - tMed; 305 final double d23 = tMed - tMin; 306 final double d13 = tMax - tMin; 307 final double dt1 = time - tMax; 308 final double dt2 = time - tMed; 309 final double dt3 = time - tMin; 310 final double iLagrange = ((dt2 * dt3 * d23) * iMax - 311 (dt1 * dt3 * d13) * iMed + 312 (dt1 * dt2 * d12) * iMin) / 313 (d12 * d23 * d13); 314 index = (int) FastMath.rint(iLagrange); 315 } 316 317 // force the next size reduction to be at least one tenth 318 final int low = FastMath.max(iMin + 1, (9 * iMin + iMax) / 10); 319 final int high = FastMath.min(iMax - 1, (iMin + 9 * iMax) / 10); 320 if (index < low) { 321 index = low; 322 } else if (index > high) { 323 index = high; 324 } 325 326 } 327 328 // now the table slice is very small, we perform an iterative search 329 index = iMin; 330 while ((index <= iMax) && (locatePoint(time, steps.get(index)) > 0)) { 331 ++index; 332 } 333 334 steps.get(index).setInterpolatedTime(time); 335 336 } 337 338 /** 339 * Get the state vector of the interpolated point. 340 * <p>The returned vector is a reference to a reused array, so 341 * it should not be modified and it should be copied if it needs 342 * to be preserved across several calls to the associated 343 * {@link #setInterpolatedTime(double)} method.</p> 344 * @return state vector at time {@link #getInterpolatedTime} 345 * @exception MaxCountExceededException if the number of functions evaluations is exceeded 346 * @see #setInterpolatedTime(double) 347 * @see #getInterpolatedDerivatives() 348 * @see #getInterpolatedSecondaryState(int) 349 * @see #getInterpolatedSecondaryDerivatives(int) 350 */ 351 public double[] getInterpolatedState() throws MaxCountExceededException { 352 return steps.get(index).getInterpolatedState(); 353 } 354 355 /** 356 * Get the derivatives of the state vector of the interpolated point. 357 * <p>The returned vector is a reference to a reused array, so 358 * it should not be modified and it should be copied if it needs 359 * to be preserved across several calls to the associated 360 * {@link #setInterpolatedTime(double)} method.</p> 361 * @return derivatives of the state vector at time {@link #getInterpolatedTime} 362 * @exception MaxCountExceededException if the number of functions evaluations is exceeded 363 * @see #setInterpolatedTime(double) 364 * @see #getInterpolatedState() 365 * @see #getInterpolatedSecondaryState(int) 366 * @see #getInterpolatedSecondaryDerivatives(int) 367 * @since 3.4 368 */ 369 public double[] getInterpolatedDerivatives() throws MaxCountExceededException { 370 return steps.get(index).getInterpolatedDerivatives(); 371 } 372 373 /** Get the interpolated secondary state corresponding to the secondary equations. 374 * <p>The returned vector is a reference to a reused array, so 375 * it should not be modified and it should be copied if it needs 376 * to be preserved across several calls to the associated 377 * {@link #setInterpolatedTime(double)} method.</p> 378 * @param secondaryStateIndex index of the secondary set, as returned by {@link 379 * org.apache.commons.math3.ode.ExpandableStatefulODE#addSecondaryEquations( 380 * org.apache.commons.math3.ode.SecondaryEquations) 381 * ExpandableStatefulODE.addSecondaryEquations(SecondaryEquations)} 382 * @return interpolated secondary state at the current interpolation date 383 * @see #setInterpolatedTime(double) 384 * @see #getInterpolatedState() 385 * @see #getInterpolatedDerivatives() 386 * @see #getInterpolatedSecondaryDerivatives(int) 387 * @since 3.2 388 * @exception MaxCountExceededException if the number of functions evaluations is exceeded 389 */ 390 public double[] getInterpolatedSecondaryState(final int secondaryStateIndex) 391 throws MaxCountExceededException { 392 return steps.get(index).getInterpolatedSecondaryState(secondaryStateIndex); 393 } 394 395 /** Get the interpolated secondary derivatives corresponding to the secondary equations. 396 * <p>The returned vector is a reference to a reused array, so 397 * it should not be modified and it should be copied if it needs 398 * to be preserved across several calls to the associated 399 * {@link #setInterpolatedTime(double)} method.</p> 400 * @param secondaryStateIndex index of the secondary set, as returned by {@link 401 * org.apache.commons.math3.ode.ExpandableStatefulODE#addSecondaryEquations( 402 * org.apache.commons.math3.ode.SecondaryEquations) 403 * ExpandableStatefulODE.addSecondaryEquations(SecondaryEquations)} 404 * @return interpolated secondary derivatives at the current interpolation date 405 * @see #setInterpolatedTime(double) 406 * @see #getInterpolatedState() 407 * @see #getInterpolatedDerivatives() 408 * @see #getInterpolatedSecondaryState(int) 409 * @since 3.4 410 * @exception MaxCountExceededException if the number of functions evaluations is exceeded 411 */ 412 public double[] getInterpolatedSecondaryDerivatives(final int secondaryStateIndex) 413 throws MaxCountExceededException { 414 return steps.get(index).getInterpolatedSecondaryDerivatives(secondaryStateIndex); 415 } 416 417 /** Compare a step interval and a double. 418 * @param time point to locate 419 * @param interval step interval 420 * @return -1 if the double is before the interval, 0 if it is in 421 * the interval, and +1 if it is after the interval, according to 422 * the interval direction 423 */ 424 private int locatePoint(final double time, final StepInterpolator interval) { 425 if (forward) { 426 if (time < interval.getPreviousTime()) { 427 return -1; 428 } else if (time > interval.getCurrentTime()) { 429 return +1; 430 } else { 431 return 0; 432 } 433 } 434 if (time > interval.getPreviousTime()) { 435 return -1; 436 } else if (time < interval.getCurrentTime()) { 437 return +1; 438 } else { 439 return 0; 440 } 441 } 442 443}