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