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; 19 20 import java.util.ArrayList; 21 import java.util.List; 22 import java.io.Serializable; 23 24 import org.apache.commons.math.MathRuntimeException; 25 import org.apache.commons.math.ode.nonstiff.AdaptiveStepsizeIntegrator; 26 import org.apache.commons.math.ode.sampling.StepHandler; 27 import org.apache.commons.math.ode.sampling.StepInterpolator; 28 29 /** 30 * This class stores all information provided by an ODE integrator 31 * during the integration process and build a continuous model of the 32 * solution from this. 33 * 34 * <p>This class act as a step handler from the integrator point of 35 * view. It is called iteratively during the integration process and 36 * stores a copy of all steps information in a sorted collection for 37 * later use. Once the integration process is over, the user can use 38 * the {@link #setInterpolatedTime setInterpolatedTime} and {@link 39 * #getInterpolatedState getInterpolatedState} to retrieve this 40 * information at any time. It is important to wait for the 41 * integration to be over before attempting to call {@link 42 * #setInterpolatedTime setInterpolatedTime} because some internal 43 * variables are set only once the last step has been handled.</p> 44 * 45 * <p>This is useful for example if the main loop of the user 46 * application should remain independent from the integration process 47 * or if one needs to mimic the behaviour of an analytical model 48 * despite a numerical model is used (i.e. one needs the ability to 49 * get the model value at any time or to navigate through the 50 * data).</p> 51 * 52 * <p>If problem modeling is done with several separate 53 * integration phases for contiguous intervals, the same 54 * ContinuousOutputModel can be used as step handler for all 55 * integration phases as long as they are performed in order and in 56 * the same direction. As an example, one can extrapolate the 57 * trajectory of a satellite with one model (i.e. one set of 58 * differential equations) up to the beginning of a maneuver, use 59 * another more complex model including thrusters modeling and 60 * accurate attitude control during the maneuver, and revert to the 61 * first model after the end of the maneuver. If the same continuous 62 * output model handles the steps of all integration phases, the user 63 * do not need to bother when the maneuver begins or ends, he has all 64 * the data available in a transparent manner.</p> 65 * 66 * <p>An important feature of this class is that it implements the 67 * <code>Serializable</code> interface. This means that the result of 68 * an integration can be serialized and reused later (if stored into a 69 * persistent medium like a filesystem or a database) or elsewhere (if 70 * sent to another application). Only the result of the integration is 71 * stored, there is no reference to the integrated problem by 72 * itself.</p> 73 * 74 * <p>One should be aware that the amount of data stored in a 75 * ContinuousOutputModel instance can be important if the state vector 76 * is large, if the integration interval is long or if the steps are 77 * small (which can result from small tolerance settings in {@link 78 * AdaptiveStepsizeIntegrator adaptive step size integrators}).</p> 79 * 80 * @see StepHandler 81 * @see StepInterpolator 82 * @version $Revision: 782431 $ $Date: 2009-06-07 15:04:37 -0400 (Sun, 07 Jun 2009) $ 83 * @since 1.2 84 */ 85 86 public class ContinuousOutputModel 87 implements StepHandler, Serializable { 88 89 /** Simple constructor. 90 * Build an empty continuous output model. 91 */ 92 public ContinuousOutputModel() { 93 steps = new ArrayList<StepInterpolator>(); 94 reset(); 95 } 96 97 /** Append another model at the end of the instance. 98 * @param model model to add at the end of the instance 99 * @exception DerivativeException if some step interpolators from 100 * the appended model cannot be copied 101 * @exception IllegalArgumentException if the model to append is not 102 * compatible with the instance (dimension of the state vector, 103 * propagation direction, hole between the dates) 104 */ 105 public void append(final ContinuousOutputModel model) 106 throws DerivativeException { 107 108 if (model.steps.size() == 0) { 109 return; 110 } 111 112 if (steps.size() == 0) { 113 initialTime = model.initialTime; 114 forward = model.forward; 115 } else { 116 117 if (getInterpolatedState().length != model.getInterpolatedState().length) { 118 throw MathRuntimeException.createIllegalArgumentException( 119 "dimension mismatch {0} != {1}", 120 getInterpolatedState().length, model.getInterpolatedState().length); 121 } 122 123 if (forward ^ model.forward) { 124 throw MathRuntimeException.createIllegalArgumentException( 125 "propagation direction mismatch"); 126 } 127 128 final StepInterpolator lastInterpolator = steps.get(index); 129 final double current = lastInterpolator.getCurrentTime(); 130 final double previous = lastInterpolator.getPreviousTime(); 131 final double step = current - previous; 132 final double gap = model.getInitialTime() - current; 133 if (Math.abs(gap) > 1.0e-3 * Math.abs(step)) { 134 throw MathRuntimeException.createIllegalArgumentException( 135 "{0} wide hole between models time ranges", Math.abs(gap)); 136 } 137 138 } 139 140 for (StepInterpolator interpolator : model.steps) { 141 steps.add(interpolator.copy()); 142 } 143 144 index = steps.size() - 1; 145 finalTime = (steps.get(index)).getCurrentTime(); 146 147 } 148 149 /** Determines whether this handler needs dense output. 150 * <p>The essence of this class is to provide dense output over all 151 * steps, hence it requires the internal steps to provide themselves 152 * dense output. The method therefore returns always true.</p> 153 * @return always true 154 */ 155 public boolean requiresDenseOutput() { 156 return true; 157 } 158 159 /** Reset the step handler. 160 * Initialize the internal data as required before the first step is 161 * handled. 162 */ 163 public void reset() { 164 initialTime = Double.NaN; 165 finalTime = Double.NaN; 166 forward = true; 167 index = 0; 168 steps.clear(); 169 } 170 171 /** Handle the last accepted step. 172 * A copy of the information provided by the last step is stored in 173 * the instance for later use. 174 * @param interpolator interpolator for the last accepted step. 175 * @param isLast true if the step is the last one 176 * @throws DerivativeException this exception is propagated to the 177 * caller if the underlying user function triggers one 178 */ 179 public void handleStep(final StepInterpolator interpolator, final boolean isLast) 180 throws DerivativeException { 181 182 if (steps.size() == 0) { 183 initialTime = interpolator.getPreviousTime(); 184 forward = interpolator.isForward(); 185 } 186 187 steps.add(interpolator.copy()); 188 189 if (isLast) { 190 finalTime = interpolator.getCurrentTime(); 191 index = steps.size() - 1; 192 } 193 194 } 195 196 /** 197 * Get the initial integration time. 198 * @return initial integration time 199 */ 200 public double getInitialTime() { 201 return initialTime; 202 } 203 204 /** 205 * Get the final integration time. 206 * @return final integration time 207 */ 208 public double getFinalTime() { 209 return finalTime; 210 } 211 212 /** 213 * Get the time of the interpolated point. 214 * If {@link #setInterpolatedTime} has not been called, it returns 215 * the final integration time. 216 * @return interpolation point time 217 */ 218 public double getInterpolatedTime() { 219 return steps.get(index).getInterpolatedTime(); 220 } 221 222 /** Set the time of the interpolated point. 223 * <p>This method should <strong>not</strong> be called before the 224 * integration is over because some internal variables are set only 225 * once the last step has been handled.</p> 226 * <p>Setting the time outside of the integration interval is now 227 * allowed (it was not allowed up to version 5.9 of Mantissa), but 228 * should be used with care since the accuracy of the interpolator 229 * will probably be very poor far from this interval. This allowance 230 * has been added to simplify implementation of search algorithms 231 * near the interval endpoints.</p> 232 * @param time time of the interpolated point 233 */ 234 public void setInterpolatedTime(final double time) { 235 236 // initialize the search with the complete steps table 237 int iMin = 0; 238 final StepInterpolator sMin = steps.get(iMin); 239 double tMin = 0.5 * (sMin.getPreviousTime() + sMin.getCurrentTime()); 240 241 int iMax = steps.size() - 1; 242 final StepInterpolator sMax = steps.get(iMax); 243 double tMax = 0.5 * (sMax.getPreviousTime() + sMax.getCurrentTime()); 244 245 // handle points outside of the integration interval 246 // or in the first and last step 247 if (locatePoint(time, sMin) <= 0) { 248 index = iMin; 249 sMin.setInterpolatedTime(time); 250 return; 251 } 252 if (locatePoint(time, sMax) >= 0) { 253 index = iMax; 254 sMax.setInterpolatedTime(time); 255 return; 256 } 257 258 // reduction of the table slice size 259 while (iMax - iMin > 5) { 260 261 // use the last estimated index as the splitting index 262 final StepInterpolator si = steps.get(index); 263 final int location = locatePoint(time, si); 264 if (location < 0) { 265 iMax = index; 266 tMax = 0.5 * (si.getPreviousTime() + si.getCurrentTime()); 267 } else if (location > 0) { 268 iMin = index; 269 tMin = 0.5 * (si.getPreviousTime() + si.getCurrentTime()); 270 } else { 271 // we have found the target step, no need to continue searching 272 si.setInterpolatedTime(time); 273 return; 274 } 275 276 // compute a new estimate of the index in the reduced table slice 277 final int iMed = (iMin + iMax) / 2; 278 final StepInterpolator sMed = steps.get(iMed); 279 final double tMed = 0.5 * (sMed.getPreviousTime() + sMed.getCurrentTime()); 280 281 if ((Math.abs(tMed - tMin) < 1e-6) || (Math.abs(tMax - tMed) < 1e-6)) { 282 // too close to the bounds, we estimate using a simple dichotomy 283 index = iMed; 284 } else { 285 // estimate the index using a reverse quadratic polynom 286 // (reverse means we have i = P(t), thus allowing to simply 287 // compute index = P(time) rather than solving a quadratic equation) 288 final double d12 = tMax - tMed; 289 final double d23 = tMed - tMin; 290 final double d13 = tMax - tMin; 291 final double dt1 = time - tMax; 292 final double dt2 = time - tMed; 293 final double dt3 = time - tMin; 294 final double iLagrange = ((dt2 * dt3 * d23) * iMax - 295 (dt1 * dt3 * d13) * iMed + 296 (dt1 * dt2 * d12) * iMin) / 297 (d12 * d23 * d13); 298 index = (int) Math.rint(iLagrange); 299 } 300 301 // force the next size reduction to be at least one tenth 302 final int low = Math.max(iMin + 1, (9 * iMin + iMax) / 10); 303 final int high = Math.min(iMax - 1, (iMin + 9 * iMax) / 10); 304 if (index < low) { 305 index = low; 306 } else if (index > high) { 307 index = high; 308 } 309 310 } 311 312 // now the table slice is very small, we perform an iterative search 313 index = iMin; 314 while ((index <= iMax) && (locatePoint(time, steps.get(index)) > 0)) { 315 ++index; 316 } 317 318 steps.get(index).setInterpolatedTime(time); 319 320 } 321 322 /** 323 * Get the state vector of the interpolated point. 324 * @return state vector at time {@link #getInterpolatedTime} 325 * @throws DerivativeException if this call induces an automatic 326 * step finalization that throws one 327 */ 328 public double[] getInterpolatedState() throws DerivativeException { 329 return steps.get(index).getInterpolatedState(); 330 } 331 332 /** Compare a step interval and a double. 333 * @param time point to locate 334 * @param interval step interval 335 * @return -1 if the double is before the interval, 0 if it is in 336 * the interval, and +1 if it is after the interval, according to 337 * the interval direction 338 */ 339 private int locatePoint(final double time, final StepInterpolator interval) { 340 if (forward) { 341 if (time < interval.getPreviousTime()) { 342 return -1; 343 } else if (time > interval.getCurrentTime()) { 344 return +1; 345 } else { 346 return 0; 347 } 348 } 349 if (time > interval.getPreviousTime()) { 350 return -1; 351 } else if (time < interval.getCurrentTime()) { 352 return +1; 353 } else { 354 return 0; 355 } 356 } 357 358 /** Initial integration time. */ 359 private double initialTime; 360 361 /** Final integration time. */ 362 private double finalTime; 363 364 /** Integration direction indicator. */ 365 private boolean forward; 366 367 /** Current interpolator index. */ 368 private int index; 369 370 /** Steps table. */ 371 private List<StepInterpolator> steps; 372 373 /** Serializable version identifier */ 374 private static final long serialVersionUID = -1417964919405031606L; 375 376 }