Raymond | dee0849 | 2015-04-02 10:43:13 -0700 | [diff] [blame] | 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.jacobians; |
| 19 | |
| 20 | import java.io.IOException; |
| 21 | import java.io.ObjectInput; |
| 22 | import java.io.ObjectOutput; |
| 23 | import java.lang.reflect.Array; |
| 24 | import java.util.ArrayList; |
| 25 | import java.util.Collection; |
| 26 | |
| 27 | import org.apache.commons.math.MathRuntimeException; |
| 28 | import org.apache.commons.math.MaxEvaluationsExceededException; |
| 29 | import org.apache.commons.math.ode.DerivativeException; |
| 30 | import org.apache.commons.math.exception.util.LocalizedFormats; |
| 31 | import org.apache.commons.math.ode.ExtendedFirstOrderDifferentialEquations; |
| 32 | import org.apache.commons.math.ode.FirstOrderIntegrator; |
| 33 | import org.apache.commons.math.ode.IntegratorException; |
| 34 | import org.apache.commons.math.ode.events.EventException; |
| 35 | import org.apache.commons.math.ode.events.EventHandler; |
| 36 | import org.apache.commons.math.ode.sampling.StepHandler; |
| 37 | import org.apache.commons.math.ode.sampling.StepInterpolator; |
| 38 | |
| 39 | /** This class enhances a first order integrator for differential equations to |
| 40 | * compute also partial derivatives of the solution with respect to initial state |
| 41 | * and parameters. |
| 42 | * <p>In order to compute both the state and its derivatives, the ODE problem |
| 43 | * is extended with jacobians of the raw ODE and the variational equations are |
| 44 | * added to form a new compound problem of higher dimension. If the original ODE |
| 45 | * problem has dimension n and there are p parameters, the compound problem will |
| 46 | * have dimension n × (1 + n + p).</p> |
| 47 | * @see ParameterizedODE |
| 48 | * @see ODEWithJacobians |
| 49 | * @version $Revision: 1073158 $ $Date: 2011-02-21 22:46:52 +0100 (lun. 21 févr. 2011) $ |
| 50 | * @since 2.1 |
| 51 | * @deprecated as of 2.2 the complete package is deprecated, it will be replaced |
| 52 | * in 3.0 by a completely rewritten implementation |
| 53 | */ |
| 54 | @Deprecated |
| 55 | public class FirstOrderIntegratorWithJacobians { |
| 56 | |
| 57 | /** Underlying integrator for compound problem. */ |
| 58 | private final FirstOrderIntegrator integrator; |
| 59 | |
| 60 | /** Raw equations to integrate. */ |
| 61 | private final ODEWithJacobians ode; |
| 62 | |
| 63 | /** Maximal number of evaluations allowed. */ |
| 64 | private int maxEvaluations; |
| 65 | |
| 66 | /** Number of evaluations already performed. */ |
| 67 | private int evaluations; |
| 68 | |
| 69 | /** Build an enhanced integrator using internal differentiation to compute jacobians. |
| 70 | * @param integrator underlying integrator to solve the compound problem |
| 71 | * @param ode original problem (f in the equation y' = f(t, y)) |
| 72 | * @param p parameters array (may be null if {@link |
| 73 | * ParameterizedODE#getParametersDimension() |
| 74 | * getParametersDimension()} from original problem is zero) |
| 75 | * @param hY step sizes to use for computing the jacobian df/dy, must have the |
| 76 | * same dimension as the original problem |
| 77 | * @param hP step sizes to use for computing the jacobian df/dp, must have the |
| 78 | * same dimension as the original problem parameters dimension |
| 79 | * @see #FirstOrderIntegratorWithJacobians(FirstOrderIntegrator, |
| 80 | * ODEWithJacobians) |
| 81 | */ |
| 82 | public FirstOrderIntegratorWithJacobians(final FirstOrderIntegrator integrator, |
| 83 | final ParameterizedODE ode, |
| 84 | final double[] p, final double[] hY, final double[] hP) { |
| 85 | checkDimension(ode.getDimension(), hY); |
| 86 | checkDimension(ode.getParametersDimension(), p); |
| 87 | checkDimension(ode.getParametersDimension(), hP); |
| 88 | this.integrator = integrator; |
| 89 | this.ode = new FiniteDifferencesWrapper(ode, p, hY, hP); |
| 90 | setMaxEvaluations(-1); |
| 91 | } |
| 92 | |
| 93 | /** Build an enhanced integrator using ODE builtin jacobian computation features. |
| 94 | * @param integrator underlying integrator to solve the compound problem |
| 95 | * @param ode original problem, which can compute the jacobians by itself |
| 96 | * @see #FirstOrderIntegratorWithJacobians(FirstOrderIntegrator, |
| 97 | * ParameterizedODE, double[], double[], double[]) |
| 98 | */ |
| 99 | public FirstOrderIntegratorWithJacobians(final FirstOrderIntegrator integrator, |
| 100 | final ODEWithJacobians ode) { |
| 101 | this.integrator = integrator; |
| 102 | this.ode = ode; |
| 103 | setMaxEvaluations(-1); |
| 104 | } |
| 105 | |
| 106 | /** Add a step handler to this integrator. |
| 107 | * <p>The handler will be called by the integrator for each accepted |
| 108 | * step.</p> |
| 109 | * @param handler handler for the accepted steps |
| 110 | * @see #getStepHandlers() |
| 111 | * @see #clearStepHandlers() |
| 112 | */ |
| 113 | public void addStepHandler(StepHandlerWithJacobians handler) { |
| 114 | final int n = ode.getDimension(); |
| 115 | final int k = ode.getParametersDimension(); |
| 116 | integrator.addStepHandler(new StepHandlerWrapper(handler, n, k)); |
| 117 | } |
| 118 | |
| 119 | /** Get all the step handlers that have been added to the integrator. |
| 120 | * @return an unmodifiable collection of the added events handlers |
| 121 | * @see #addStepHandler(StepHandlerWithJacobians) |
| 122 | * @see #clearStepHandlers() |
| 123 | */ |
| 124 | public Collection<StepHandlerWithJacobians> getStepHandlers() { |
| 125 | final Collection<StepHandlerWithJacobians> handlers = |
| 126 | new ArrayList<StepHandlerWithJacobians>(); |
| 127 | for (final StepHandler handler : integrator.getStepHandlers()) { |
| 128 | if (handler instanceof StepHandlerWrapper) { |
| 129 | handlers.add(((StepHandlerWrapper) handler).getHandler()); |
| 130 | } |
| 131 | } |
| 132 | return handlers; |
| 133 | } |
| 134 | |
| 135 | /** Remove all the step handlers that have been added to the integrator. |
| 136 | * @see #addStepHandler(StepHandlerWithJacobians) |
| 137 | * @see #getStepHandlers() |
| 138 | */ |
| 139 | public void clearStepHandlers() { |
| 140 | integrator.clearStepHandlers(); |
| 141 | } |
| 142 | |
| 143 | /** Add an event handler to the integrator. |
| 144 | * @param handler event handler |
| 145 | * @param maxCheckInterval maximal time interval between switching |
| 146 | * function checks (this interval prevents missing sign changes in |
| 147 | * case the integration steps becomes very large) |
| 148 | * @param convergence convergence threshold in the event time search |
| 149 | * @param maxIterationCount upper limit of the iteration count in |
| 150 | * the event time search |
| 151 | * @see #getEventHandlers() |
| 152 | * @see #clearEventHandlers() |
| 153 | */ |
| 154 | public void addEventHandler(EventHandlerWithJacobians handler, |
| 155 | double maxCheckInterval, |
| 156 | double convergence, |
| 157 | int maxIterationCount) { |
| 158 | final int n = ode.getDimension(); |
| 159 | final int k = ode.getParametersDimension(); |
| 160 | integrator.addEventHandler(new EventHandlerWrapper(handler, n, k), |
| 161 | maxCheckInterval, convergence, maxIterationCount); |
| 162 | } |
| 163 | |
| 164 | /** Get all the event handlers that have been added to the integrator. |
| 165 | * @return an unmodifiable collection of the added events handlers |
| 166 | * @see #addEventHandler(EventHandlerWithJacobians, double, double, int) |
| 167 | * @see #clearEventHandlers() |
| 168 | */ |
| 169 | public Collection<EventHandlerWithJacobians> getEventHandlers() { |
| 170 | final Collection<EventHandlerWithJacobians> handlers = |
| 171 | new ArrayList<EventHandlerWithJacobians>(); |
| 172 | for (final EventHandler handler : integrator.getEventHandlers()) { |
| 173 | if (handler instanceof EventHandlerWrapper) { |
| 174 | handlers.add(((EventHandlerWrapper) handler).getHandler()); |
| 175 | } |
| 176 | } |
| 177 | return handlers; |
| 178 | } |
| 179 | |
| 180 | /** Remove all the event handlers that have been added to the integrator. |
| 181 | * @see #addEventHandler(EventHandlerWithJacobians, double, double, int) |
| 182 | * @see #getEventHandlers() |
| 183 | */ |
| 184 | public void clearEventHandlers() { |
| 185 | integrator.clearEventHandlers(); |
| 186 | } |
| 187 | |
| 188 | /** Integrate the differential equations and the variational equations up to the given time. |
| 189 | * <p>This method solves an Initial Value Problem (IVP) and also computes the derivatives |
| 190 | * of the solution with respect to initial state and parameters. This can be used as |
| 191 | * a basis to solve Boundary Value Problems (BVP).</p> |
| 192 | * <p>Since this method stores some internal state variables made |
| 193 | * available in its public interface during integration ({@link |
| 194 | * #getCurrentSignedStepsize()}), it is <em>not</em> thread-safe.</p> |
| 195 | * @param t0 initial time |
| 196 | * @param y0 initial value of the state vector at t0 |
| 197 | * @param dY0dP initial value of the state vector derivative with respect to the |
| 198 | * parameters at t0 |
| 199 | * @param t target time for the integration |
| 200 | * (can be set to a value smaller than <code>t0</code> for backward integration) |
| 201 | * @param y placeholder where to put the state vector at each successful |
| 202 | * step (and hence at the end of integration), can be the same object as y0 |
| 203 | * @param dYdY0 placeholder where to put the state vector derivative with respect |
| 204 | * to the initial state (dy[i]/dy0[j] is in element array dYdY0[i][j]) at each successful |
| 205 | * step (and hence at the end of integration) |
| 206 | * @param dYdP placeholder where to put the state vector derivative with respect |
| 207 | * to the parameters (dy[i]/dp[j] is in element array dYdP[i][j]) at each successful |
| 208 | * step (and hence at the end of integration) |
| 209 | * @return stop time, will be the same as target time if integration reached its |
| 210 | * target, but may be different if some event handler stops it at some point. |
| 211 | * @throws IntegratorException if the integrator cannot perform integration |
| 212 | * @throws DerivativeException this exception is propagated to the caller if |
| 213 | * the underlying user function triggers one |
| 214 | */ |
| 215 | public double integrate(final double t0, final double[] y0, final double[][] dY0dP, |
| 216 | final double t, final double[] y, |
| 217 | final double[][] dYdY0, final double[][] dYdP) |
| 218 | throws DerivativeException, IntegratorException { |
| 219 | |
| 220 | final int n = ode.getDimension(); |
| 221 | final int k = ode.getParametersDimension(); |
| 222 | checkDimension(n, y0); |
| 223 | checkDimension(n, y); |
| 224 | checkDimension(n, dYdY0); |
| 225 | checkDimension(n, dYdY0[0]); |
| 226 | if (k != 0) { |
| 227 | checkDimension(n, dY0dP); |
| 228 | checkDimension(k, dY0dP[0]); |
| 229 | checkDimension(n, dYdP); |
| 230 | checkDimension(k, dYdP[0]); |
| 231 | } |
| 232 | |
| 233 | // set up initial state, including partial derivatives |
| 234 | // the compound state z contains the raw state y and its derivatives |
| 235 | // with respect to initial state y0 and to parameters p |
| 236 | // y[i] is stored in z[i] |
| 237 | // dy[i]/dy0[j] is stored in z[n + i * n + j] |
| 238 | // dy[i]/dp[j] is stored in z[n * (n + 1) + i * k + j] |
| 239 | final double[] z = new double[n * (1 + n + k)]; |
| 240 | System.arraycopy(y0, 0, z, 0, n); |
| 241 | for (int i = 0; i < n; ++i) { |
| 242 | |
| 243 | // set diagonal element of dy/dy0 to 1.0 at t = t0 |
| 244 | z[i * (1 + n) + n] = 1.0; |
| 245 | |
| 246 | // set initial derivatives with respect to parameters |
| 247 | System.arraycopy(dY0dP[i], 0, z, n * (n + 1) + i * k, k); |
| 248 | |
| 249 | } |
| 250 | |
| 251 | // integrate the compound state variational equations |
| 252 | evaluations = 0; |
| 253 | final double stopTime = integrator.integrate(new MappingWrapper(), t0, z, t, z); |
| 254 | |
| 255 | // dispatch the final compound state into the state and partial derivatives arrays |
| 256 | dispatchCompoundState(z, y, dYdY0, dYdP); |
| 257 | |
| 258 | return stopTime; |
| 259 | |
| 260 | } |
| 261 | |
| 262 | /** Dispatch a compound state array into state and jacobians arrays. |
| 263 | * @param z compound state |
| 264 | * @param y raw state array to fill |
| 265 | * @param dydy0 jacobian array to fill |
| 266 | * @param dydp jacobian array to fill |
| 267 | */ |
| 268 | private static void dispatchCompoundState(final double[] z, final double[] y, |
| 269 | final double[][] dydy0, final double[][] dydp) { |
| 270 | |
| 271 | final int n = y.length; |
| 272 | final int k = dydp[0].length; |
| 273 | |
| 274 | // state |
| 275 | System.arraycopy(z, 0, y, 0, n); |
| 276 | |
| 277 | // jacobian with respect to initial state |
| 278 | for (int i = 0; i < n; ++i) { |
| 279 | System.arraycopy(z, n * (i + 1), dydy0[i], 0, n); |
| 280 | } |
| 281 | |
| 282 | // jacobian with respect to parameters |
| 283 | for (int i = 0; i < n; ++i) { |
| 284 | System.arraycopy(z, n * (n + 1) + i * k, dydp[i], 0, k); |
| 285 | } |
| 286 | |
| 287 | } |
| 288 | |
| 289 | /** Get the current value of the step start time t<sub>i</sub>. |
| 290 | * <p>This method can be called during integration (typically by |
| 291 | * the object implementing the {@link org.apache.commons.math.ode.FirstOrderDifferentialEquations |
| 292 | * differential equations} problem) if the value of the current step that |
| 293 | * is attempted is needed.</p> |
| 294 | * <p>The result is undefined if the method is called outside of |
| 295 | * calls to <code>integrate</code>.</p> |
| 296 | * @return current value of the step start time t<sub>i</sub> |
| 297 | */ |
| 298 | public double getCurrentStepStart() { |
| 299 | return integrator.getCurrentStepStart(); |
| 300 | } |
| 301 | |
| 302 | /** Get the current signed value of the integration stepsize. |
| 303 | * <p>This method can be called during integration (typically by |
| 304 | * the object implementing the {@link org.apache.commons.math.ode.FirstOrderDifferentialEquations |
| 305 | * differential equations} problem) if the signed value of the current stepsize |
| 306 | * that is tried is needed.</p> |
| 307 | * <p>The result is undefined if the method is called outside of |
| 308 | * calls to <code>integrate</code>.</p> |
| 309 | * @return current signed value of the stepsize |
| 310 | */ |
| 311 | public double getCurrentSignedStepsize() { |
| 312 | return integrator.getCurrentSignedStepsize(); |
| 313 | } |
| 314 | |
| 315 | /** Set the maximal number of differential equations function evaluations. |
| 316 | * <p>The purpose of this method is to avoid infinite loops which can occur |
| 317 | * for example when stringent error constraints are set or when lots of |
| 318 | * discrete events are triggered, thus leading to many rejected steps.</p> |
| 319 | * @param maxEvaluations maximal number of function evaluations (negative |
| 320 | * values are silently converted to maximal integer value, thus representing |
| 321 | * almost unlimited evaluations) |
| 322 | */ |
| 323 | public void setMaxEvaluations(int maxEvaluations) { |
| 324 | this.maxEvaluations = (maxEvaluations < 0) ? Integer.MAX_VALUE : maxEvaluations; |
| 325 | } |
| 326 | |
| 327 | /** Get the maximal number of functions evaluations. |
| 328 | * @return maximal number of functions evaluations |
| 329 | */ |
| 330 | public int getMaxEvaluations() { |
| 331 | return maxEvaluations; |
| 332 | } |
| 333 | |
| 334 | /** Get the number of evaluations of the differential equations function. |
| 335 | * <p> |
| 336 | * The number of evaluations corresponds to the last call to the |
| 337 | * <code>integrate</code> method. It is 0 if the method has not been called yet. |
| 338 | * </p> |
| 339 | * @return number of evaluations of the differential equations function |
| 340 | */ |
| 341 | public int getEvaluations() { |
| 342 | return evaluations; |
| 343 | } |
| 344 | |
| 345 | /** Check array dimensions. |
| 346 | * @param expected expected dimension |
| 347 | * @param array (may be null if expected is 0) |
| 348 | * @throws IllegalArgumentException if the array dimension does not match the expected one |
| 349 | */ |
| 350 | private void checkDimension(final int expected, final Object array) |
| 351 | throws IllegalArgumentException { |
| 352 | int arrayDimension = (array == null) ? 0 : Array.getLength(array); |
| 353 | if (arrayDimension != expected) { |
| 354 | throw MathRuntimeException.createIllegalArgumentException( |
| 355 | LocalizedFormats.DIMENSIONS_MISMATCH_SIMPLE, arrayDimension, expected); |
| 356 | } |
| 357 | } |
| 358 | |
| 359 | /** Wrapper class used to map state and jacobians into compound state. */ |
| 360 | private class MappingWrapper implements ExtendedFirstOrderDifferentialEquations { |
| 361 | |
| 362 | /** Current state. */ |
| 363 | private final double[] y; |
| 364 | |
| 365 | /** Time derivative of the current state. */ |
| 366 | private final double[] yDot; |
| 367 | |
| 368 | /** Derivatives of yDot with respect to state. */ |
| 369 | private final double[][] dFdY; |
| 370 | |
| 371 | /** Derivatives of yDot with respect to parameters. */ |
| 372 | private final double[][] dFdP; |
| 373 | |
| 374 | /** Simple constructor. |
| 375 | */ |
| 376 | public MappingWrapper() { |
| 377 | |
| 378 | final int n = ode.getDimension(); |
| 379 | final int k = ode.getParametersDimension(); |
| 380 | y = new double[n]; |
| 381 | yDot = new double[n]; |
| 382 | dFdY = new double[n][n]; |
| 383 | dFdP = new double[n][k]; |
| 384 | |
| 385 | } |
| 386 | |
| 387 | /** {@inheritDoc} */ |
| 388 | public int getDimension() { |
| 389 | final int n = y.length; |
| 390 | final int k = dFdP[0].length; |
| 391 | return n * (1 + n + k); |
| 392 | } |
| 393 | |
| 394 | /** {@inheritDoc} */ |
| 395 | public int getMainSetDimension() { |
| 396 | return ode.getDimension(); |
| 397 | } |
| 398 | |
| 399 | /** {@inheritDoc} */ |
| 400 | public void computeDerivatives(final double t, final double[] z, final double[] zDot) |
| 401 | throws DerivativeException { |
| 402 | |
| 403 | final int n = y.length; |
| 404 | final int k = dFdP[0].length; |
| 405 | |
| 406 | // compute raw ODE and its jacobians: dy/dt, d[dy/dt]/dy0 and d[dy/dt]/dp |
| 407 | System.arraycopy(z, 0, y, 0, n); |
| 408 | if (++evaluations > maxEvaluations) { |
| 409 | throw new DerivativeException(new MaxEvaluationsExceededException(maxEvaluations)); |
| 410 | } |
| 411 | ode.computeDerivatives(t, y, yDot); |
| 412 | ode.computeJacobians(t, y, yDot, dFdY, dFdP); |
| 413 | |
| 414 | // state part of the compound equations |
| 415 | System.arraycopy(yDot, 0, zDot, 0, n); |
| 416 | |
| 417 | // variational equations: from d[dy/dt]/dy0 to d[dy/dy0]/dt |
| 418 | for (int i = 0; i < n; ++i) { |
| 419 | final double[] dFdYi = dFdY[i]; |
| 420 | for (int j = 0; j < n; ++j) { |
| 421 | double s = 0; |
| 422 | final int startIndex = n + j; |
| 423 | int zIndex = startIndex; |
| 424 | for (int l = 0; l < n; ++l) { |
| 425 | s += dFdYi[l] * z[zIndex]; |
| 426 | zIndex += n; |
| 427 | } |
| 428 | zDot[startIndex + i * n] = s; |
| 429 | } |
| 430 | } |
| 431 | |
| 432 | // variational equations: from d[dy/dt]/dy0 and d[dy/dt]/dp to d[dy/dp]/dt |
| 433 | for (int i = 0; i < n; ++i) { |
| 434 | final double[] dFdYi = dFdY[i]; |
| 435 | final double[] dFdPi = dFdP[i]; |
| 436 | for (int j = 0; j < k; ++j) { |
| 437 | double s = dFdPi[j]; |
| 438 | final int startIndex = n * (n + 1) + j; |
| 439 | int zIndex = startIndex; |
| 440 | for (int l = 0; l < n; ++l) { |
| 441 | s += dFdYi[l] * z[zIndex]; |
| 442 | zIndex += k; |
| 443 | } |
| 444 | zDot[startIndex + i * k] = s; |
| 445 | } |
| 446 | } |
| 447 | |
| 448 | } |
| 449 | |
| 450 | } |
| 451 | |
| 452 | /** Wrapper class to compute jacobians by finite differences for ODE which do not compute them themselves. */ |
| 453 | private class FiniteDifferencesWrapper implements ODEWithJacobians { |
| 454 | |
| 455 | /** Raw ODE without jacobians computation. */ |
| 456 | private final ParameterizedODE ode; |
| 457 | |
| 458 | /** Parameters array (may be null if parameters dimension from original problem is zero) */ |
| 459 | private final double[] p; |
| 460 | |
| 461 | /** Step sizes to use for computing the jacobian df/dy. */ |
| 462 | private final double[] hY; |
| 463 | |
| 464 | /** Step sizes to use for computing the jacobian df/dp. */ |
| 465 | private final double[] hP; |
| 466 | |
| 467 | /** Temporary array for state derivatives used to compute jacobians. */ |
| 468 | private final double[] tmpDot; |
| 469 | |
| 470 | /** Simple constructor. |
| 471 | * @param ode original ODE problem, without jacobians computations |
| 472 | * @param p parameters array (may be null if parameters dimension from original problem is zero) |
| 473 | * @param hY step sizes to use for computing the jacobian df/dy |
| 474 | * @param hP step sizes to use for computing the jacobian df/dp |
| 475 | */ |
| 476 | public FiniteDifferencesWrapper(final ParameterizedODE ode, |
| 477 | final double[] p, final double[] hY, final double[] hP) { |
| 478 | this.ode = ode; |
| 479 | this.p = p.clone(); |
| 480 | this.hY = hY.clone(); |
| 481 | this.hP = hP.clone(); |
| 482 | tmpDot = new double[ode.getDimension()]; |
| 483 | } |
| 484 | |
| 485 | /** {@inheritDoc} */ |
| 486 | public int getDimension() { |
| 487 | return ode.getDimension(); |
| 488 | } |
| 489 | |
| 490 | /** {@inheritDoc} */ |
| 491 | public void computeDerivatives(double t, double[] y, double[] yDot) throws DerivativeException { |
| 492 | // this call to computeDerivatives has already been counted, |
| 493 | // we must not increment the counter again |
| 494 | ode.computeDerivatives(t, y, yDot); |
| 495 | } |
| 496 | |
| 497 | /** {@inheritDoc} */ |
| 498 | public int getParametersDimension() { |
| 499 | return ode.getParametersDimension(); |
| 500 | } |
| 501 | |
| 502 | /** {@inheritDoc} */ |
| 503 | public void computeJacobians(double t, double[] y, double[] yDot, |
| 504 | double[][] dFdY, double[][] dFdP) |
| 505 | throws DerivativeException { |
| 506 | |
| 507 | final int n = hY.length; |
| 508 | final int k = hP.length; |
| 509 | |
| 510 | evaluations += n + k; |
| 511 | if (evaluations > maxEvaluations) { |
| 512 | throw new DerivativeException(new MaxEvaluationsExceededException(maxEvaluations)); |
| 513 | } |
| 514 | |
| 515 | // compute df/dy where f is the ODE and y is the state array |
| 516 | for (int j = 0; j < n; ++j) { |
| 517 | final double savedYj = y[j]; |
| 518 | y[j] += hY[j]; |
| 519 | ode.computeDerivatives(t, y, tmpDot); |
| 520 | for (int i = 0; i < n; ++i) { |
| 521 | dFdY[i][j] = (tmpDot[i] - yDot[i]) / hY[j]; |
| 522 | } |
| 523 | y[j] = savedYj; |
| 524 | } |
| 525 | |
| 526 | // compute df/dp where f is the ODE and p is the parameters array |
| 527 | for (int j = 0; j < k; ++j) { |
| 528 | ode.setParameter(j, p[j] + hP[j]); |
| 529 | ode.computeDerivatives(t, y, tmpDot); |
| 530 | for (int i = 0; i < n; ++i) { |
| 531 | dFdP[i][j] = (tmpDot[i] - yDot[i]) / hP[j]; |
| 532 | } |
| 533 | ode.setParameter(j, p[j]); |
| 534 | } |
| 535 | |
| 536 | } |
| 537 | |
| 538 | } |
| 539 | |
| 540 | /** Wrapper for step handlers. */ |
| 541 | private static class StepHandlerWrapper implements StepHandler { |
| 542 | |
| 543 | /** Underlying step handler with jacobians. */ |
| 544 | private final StepHandlerWithJacobians handler; |
| 545 | |
| 546 | /** Dimension of the original ODE. */ |
| 547 | private final int n; |
| 548 | |
| 549 | /** Number of parameters. */ |
| 550 | private final int k; |
| 551 | |
| 552 | /** Simple constructor. |
| 553 | * @param handler underlying step handler with jacobians |
| 554 | * @param n dimension of the original ODE |
| 555 | * @param k number of parameters |
| 556 | */ |
| 557 | public StepHandlerWrapper(final StepHandlerWithJacobians handler, |
| 558 | final int n, final int k) { |
| 559 | this.handler = handler; |
| 560 | this.n = n; |
| 561 | this.k = k; |
| 562 | } |
| 563 | |
| 564 | /** Get the underlying step handler with jacobians. |
| 565 | * @return underlying step handler with jacobians |
| 566 | */ |
| 567 | public StepHandlerWithJacobians getHandler() { |
| 568 | return handler; |
| 569 | } |
| 570 | |
| 571 | /** {@inheritDoc} */ |
| 572 | public void handleStep(StepInterpolator interpolator, boolean isLast) |
| 573 | throws DerivativeException { |
| 574 | handler.handleStep(new StepInterpolatorWrapper(interpolator, n, k), isLast); |
| 575 | } |
| 576 | |
| 577 | /** {@inheritDoc} */ |
| 578 | public boolean requiresDenseOutput() { |
| 579 | return handler.requiresDenseOutput(); |
| 580 | } |
| 581 | |
| 582 | /** {@inheritDoc} */ |
| 583 | public void reset() { |
| 584 | handler.reset(); |
| 585 | } |
| 586 | |
| 587 | } |
| 588 | |
| 589 | /** Wrapper for step interpolators. */ |
| 590 | private static class StepInterpolatorWrapper |
| 591 | implements StepInterpolatorWithJacobians { |
| 592 | |
| 593 | /** Wrapped interpolator. */ |
| 594 | private StepInterpolator interpolator; |
| 595 | |
| 596 | /** State array. */ |
| 597 | private double[] y; |
| 598 | |
| 599 | /** Jacobian with respect to initial state dy/dy0. */ |
| 600 | private double[][] dydy0; |
| 601 | |
| 602 | /** Jacobian with respect to parameters dy/dp. */ |
| 603 | private double[][] dydp; |
| 604 | |
| 605 | /** Time derivative of the state array. */ |
| 606 | private double[] yDot; |
| 607 | |
| 608 | /** Time derivative of the sacobian with respect to initial state dy/dy0. */ |
| 609 | private double[][] dydy0Dot; |
| 610 | |
| 611 | /** Time derivative of the jacobian with respect to parameters dy/dp. */ |
| 612 | private double[][] dydpDot; |
| 613 | |
| 614 | /** Simple constructor. |
| 615 | * <p>This constructor is used only for externalization. It does nothing.</p> |
| 616 | */ |
| 617 | @SuppressWarnings("unused") |
| 618 | public StepInterpolatorWrapper() { |
| 619 | } |
| 620 | |
| 621 | /** Simple constructor. |
| 622 | * @param interpolator wrapped interpolator |
| 623 | * @param n dimension of the original ODE |
| 624 | * @param k number of parameters |
| 625 | */ |
| 626 | public StepInterpolatorWrapper(final StepInterpolator interpolator, |
| 627 | final int n, final int k) { |
| 628 | this.interpolator = interpolator; |
| 629 | y = new double[n]; |
| 630 | dydy0 = new double[n][n]; |
| 631 | dydp = new double[n][k]; |
| 632 | yDot = new double[n]; |
| 633 | dydy0Dot = new double[n][n]; |
| 634 | dydpDot = new double[n][k]; |
| 635 | } |
| 636 | |
| 637 | /** {@inheritDoc} */ |
| 638 | public void setInterpolatedTime(double time) { |
| 639 | interpolator.setInterpolatedTime(time); |
| 640 | } |
| 641 | |
| 642 | /** {@inheritDoc} */ |
| 643 | public boolean isForward() { |
| 644 | return interpolator.isForward(); |
| 645 | } |
| 646 | |
| 647 | /** {@inheritDoc} */ |
| 648 | public double getPreviousTime() { |
| 649 | return interpolator.getPreviousTime(); |
| 650 | } |
| 651 | |
| 652 | /** {@inheritDoc} */ |
| 653 | public double getInterpolatedTime() { |
| 654 | return interpolator.getInterpolatedTime(); |
| 655 | } |
| 656 | |
| 657 | /** {@inheritDoc} */ |
| 658 | public double[] getInterpolatedY() throws DerivativeException { |
| 659 | double[] extendedState = interpolator.getInterpolatedState(); |
| 660 | System.arraycopy(extendedState, 0, y, 0, y.length); |
| 661 | return y; |
| 662 | } |
| 663 | |
| 664 | /** {@inheritDoc} */ |
| 665 | public double[][] getInterpolatedDyDy0() throws DerivativeException { |
| 666 | double[] extendedState = interpolator.getInterpolatedState(); |
| 667 | final int n = y.length; |
| 668 | int start = n; |
| 669 | for (int i = 0; i < n; ++i) { |
| 670 | System.arraycopy(extendedState, start, dydy0[i], 0, n); |
| 671 | start += n; |
| 672 | } |
| 673 | return dydy0; |
| 674 | } |
| 675 | |
| 676 | /** {@inheritDoc} */ |
| 677 | public double[][] getInterpolatedDyDp() throws DerivativeException { |
| 678 | double[] extendedState = interpolator.getInterpolatedState(); |
| 679 | final int n = y.length; |
| 680 | final int k = dydp[0].length; |
| 681 | int start = n * (n + 1); |
| 682 | for (int i = 0; i < n; ++i) { |
| 683 | System.arraycopy(extendedState, start, dydp[i], 0, k); |
| 684 | start += k; |
| 685 | } |
| 686 | return dydp; |
| 687 | } |
| 688 | |
| 689 | /** {@inheritDoc} */ |
| 690 | public double[] getInterpolatedYDot() throws DerivativeException { |
| 691 | double[] extendedDerivatives = interpolator.getInterpolatedDerivatives(); |
| 692 | System.arraycopy(extendedDerivatives, 0, yDot, 0, yDot.length); |
| 693 | return yDot; |
| 694 | } |
| 695 | |
| 696 | /** {@inheritDoc} */ |
| 697 | public double[][] getInterpolatedDyDy0Dot() throws DerivativeException { |
| 698 | double[] extendedDerivatives = interpolator.getInterpolatedDerivatives(); |
| 699 | final int n = y.length; |
| 700 | int start = n; |
| 701 | for (int i = 0; i < n; ++i) { |
| 702 | System.arraycopy(extendedDerivatives, start, dydy0Dot[i], 0, n); |
| 703 | start += n; |
| 704 | } |
| 705 | return dydy0Dot; |
| 706 | } |
| 707 | |
| 708 | /** {@inheritDoc} */ |
| 709 | public double[][] getInterpolatedDyDpDot() throws DerivativeException { |
| 710 | double[] extendedDerivatives = interpolator.getInterpolatedDerivatives(); |
| 711 | final int n = y.length; |
| 712 | final int k = dydpDot[0].length; |
| 713 | int start = n * (n + 1); |
| 714 | for (int i = 0; i < n; ++i) { |
| 715 | System.arraycopy(extendedDerivatives, start, dydpDot[i], 0, k); |
| 716 | start += k; |
| 717 | } |
| 718 | return dydpDot; |
| 719 | } |
| 720 | |
| 721 | /** {@inheritDoc} */ |
| 722 | public double getCurrentTime() { |
| 723 | return interpolator.getCurrentTime(); |
| 724 | } |
| 725 | |
| 726 | /** {@inheritDoc} */ |
| 727 | public StepInterpolatorWithJacobians copy() throws DerivativeException { |
| 728 | final int n = y.length; |
| 729 | final int k = dydp[0].length; |
| 730 | StepInterpolatorWrapper copied = |
| 731 | new StepInterpolatorWrapper(interpolator.copy(), n, k); |
| 732 | copyArray(y, copied.y); |
| 733 | copyArray(dydy0, copied.dydy0); |
| 734 | copyArray(dydp, copied.dydp); |
| 735 | copyArray(yDot, copied.yDot); |
| 736 | copyArray(dydy0Dot, copied.dydy0Dot); |
| 737 | copyArray(dydpDot, copied.dydpDot); |
| 738 | return copied; |
| 739 | } |
| 740 | |
| 741 | /** {@inheritDoc} */ |
| 742 | public void writeExternal(ObjectOutput out) throws IOException { |
| 743 | out.writeObject(interpolator); |
| 744 | out.writeInt(y.length); |
| 745 | out.writeInt(dydp[0].length); |
| 746 | writeArray(out, y); |
| 747 | writeArray(out, dydy0); |
| 748 | writeArray(out, dydp); |
| 749 | writeArray(out, yDot); |
| 750 | writeArray(out, dydy0Dot); |
| 751 | writeArray(out, dydpDot); |
| 752 | } |
| 753 | |
| 754 | /** {@inheritDoc} */ |
| 755 | public void readExternal(ObjectInput in) throws IOException, ClassNotFoundException { |
| 756 | interpolator = (StepInterpolator) in.readObject(); |
| 757 | final int n = in.readInt(); |
| 758 | final int k = in.readInt(); |
| 759 | y = new double[n]; |
| 760 | dydy0 = new double[n][n]; |
| 761 | dydp = new double[n][k]; |
| 762 | yDot = new double[n]; |
| 763 | dydy0Dot = new double[n][n]; |
| 764 | dydpDot = new double[n][k]; |
| 765 | readArray(in, y); |
| 766 | readArray(in, dydy0); |
| 767 | readArray(in, dydp); |
| 768 | readArray(in, yDot); |
| 769 | readArray(in, dydy0Dot); |
| 770 | readArray(in, dydpDot); |
| 771 | } |
| 772 | |
| 773 | /** Copy an array. |
| 774 | * @param src source array |
| 775 | * @param dest destination array |
| 776 | */ |
| 777 | private static void copyArray(final double[] src, final double[] dest) { |
| 778 | System.arraycopy(src, 0, dest, 0, src.length); |
| 779 | } |
| 780 | |
| 781 | /** Copy an array. |
| 782 | * @param src source array |
| 783 | * @param dest destination array |
| 784 | */ |
| 785 | private static void copyArray(final double[][] src, final double[][] dest) { |
| 786 | for (int i = 0; i < src.length; ++i) { |
| 787 | copyArray(src[i], dest[i]); |
| 788 | } |
| 789 | } |
| 790 | |
| 791 | /** Write an array. |
| 792 | * @param out output stream |
| 793 | * @param array array to write |
| 794 | * @exception IOException if array cannot be read |
| 795 | */ |
| 796 | private static void writeArray(final ObjectOutput out, final double[] array) |
| 797 | throws IOException { |
| 798 | for (int i = 0; i < array.length; ++i) { |
| 799 | out.writeDouble(array[i]); |
| 800 | } |
| 801 | } |
| 802 | |
| 803 | /** Write an array. |
| 804 | * @param out output stream |
| 805 | * @param array array to write |
| 806 | * @exception IOException if array cannot be read |
| 807 | */ |
| 808 | private static void writeArray(final ObjectOutput out, final double[][] array) |
| 809 | throws IOException { |
| 810 | for (int i = 0; i < array.length; ++i) { |
| 811 | writeArray(out, array[i]); |
| 812 | } |
| 813 | } |
| 814 | |
| 815 | /** Read an array. |
| 816 | * @param in input stream |
| 817 | * @param array array to read |
| 818 | * @exception IOException if array cannot be read |
| 819 | */ |
| 820 | private static void readArray(final ObjectInput in, final double[] array) |
| 821 | throws IOException { |
| 822 | for (int i = 0; i < array.length; ++i) { |
| 823 | array[i] = in.readDouble(); |
| 824 | } |
| 825 | } |
| 826 | |
| 827 | /** Read an array. |
| 828 | * @param in input stream |
| 829 | * @param array array to read |
| 830 | * @exception IOException if array cannot be read |
| 831 | */ |
| 832 | private static void readArray(final ObjectInput in, final double[][] array) |
| 833 | throws IOException { |
| 834 | for (int i = 0; i < array.length; ++i) { |
| 835 | readArray(in, array[i]); |
| 836 | } |
| 837 | } |
| 838 | |
| 839 | } |
| 840 | |
| 841 | /** Wrapper for event handlers. */ |
| 842 | private static class EventHandlerWrapper implements EventHandler { |
| 843 | |
| 844 | /** Underlying event handler with jacobians. */ |
| 845 | private final EventHandlerWithJacobians handler; |
| 846 | |
| 847 | /** State array. */ |
| 848 | private double[] y; |
| 849 | |
| 850 | /** Jacobian with respect to initial state dy/dy0. */ |
| 851 | private double[][] dydy0; |
| 852 | |
| 853 | /** Jacobian with respect to parameters dy/dp. */ |
| 854 | private double[][] dydp; |
| 855 | |
| 856 | /** Simple constructor. |
| 857 | * @param handler underlying event handler with jacobians |
| 858 | * @param n dimension of the original ODE |
| 859 | * @param k number of parameters |
| 860 | */ |
| 861 | public EventHandlerWrapper(final EventHandlerWithJacobians handler, |
| 862 | final int n, final int k) { |
| 863 | this.handler = handler; |
| 864 | y = new double[n]; |
| 865 | dydy0 = new double[n][n]; |
| 866 | dydp = new double[n][k]; |
| 867 | } |
| 868 | |
| 869 | /** Get the underlying event handler with jacobians. |
| 870 | * @return underlying event handler with jacobians |
| 871 | */ |
| 872 | public EventHandlerWithJacobians getHandler() { |
| 873 | return handler; |
| 874 | } |
| 875 | |
| 876 | /** {@inheritDoc} */ |
| 877 | public int eventOccurred(double t, double[] z, boolean increasing) |
| 878 | throws EventException { |
| 879 | dispatchCompoundState(z, y, dydy0, dydp); |
| 880 | return handler.eventOccurred(t, y, dydy0, dydp, increasing); |
| 881 | } |
| 882 | |
| 883 | /** {@inheritDoc} */ |
| 884 | public double g(double t, double[] z) |
| 885 | throws EventException { |
| 886 | dispatchCompoundState(z, y, dydy0, dydp); |
| 887 | return handler.g(t, y, dydy0, dydp); |
| 888 | } |
| 889 | |
| 890 | /** {@inheritDoc} */ |
| 891 | public void resetState(double t, double[] z) |
| 892 | throws EventException { |
| 893 | dispatchCompoundState(z, y, dydy0, dydp); |
| 894 | handler.resetState(t, y, dydy0, dydp); |
| 895 | } |
| 896 | |
| 897 | } |
| 898 | |
| 899 | } |