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