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 */ FirstOrderIntegratorWithJacobians(final FirstOrderIntegrator integrator, final ParameterizedODE ode, final double[] p, final double[] hY, final double[] hP)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 */ FirstOrderIntegratorWithJacobians(final FirstOrderIntegrator integrator, final ODEWithJacobians ode)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 */ addStepHandler(StepHandlerWithJacobians handler)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 */ getStepHandlers()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 */ clearStepHandlers()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 */ addEventHandler(EventHandlerWithJacobians handler, double maxCheckInterval, double convergence, int maxIterationCount)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 */ getEventHandlers()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 */ clearEventHandlers()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 */ integrate(final double t0, final double[] y0, final double[][] dY0dP, final double t, final double[] y, final double[][] dYdY0, final double[][] dYdP)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 */ dispatchCompoundState(final double[] z, final double[] y, final double[][] dydy0, final double[][] dydp)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 */ getCurrentStepStart()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 */ getCurrentSignedStepsize()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 */ setMaxEvaluations(int maxEvaluations)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 */ getMaxEvaluations()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 */ getEvaluations()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 */ checkDimension(final int expected, final Object array)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 */ MappingWrapper()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} */ getDimension()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} */ getMainSetDimension()395 public int getMainSetDimension() { 396 return ode.getDimension(); 397 } 398 399 /** {@inheritDoc} */ computeDerivatives(final double t, final double[] z, final double[] zDot)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 */ FiniteDifferencesWrapper(final ParameterizedODE ode, final double[] p, final double[] hY, final double[] hP)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} */ getDimension()486 public int getDimension() { 487 return ode.getDimension(); 488 } 489 490 /** {@inheritDoc} */ computeDerivatives(double t, double[] y, double[] yDot)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} */ getParametersDimension()498 public int getParametersDimension() { 499 return ode.getParametersDimension(); 500 } 501 502 /** {@inheritDoc} */ computeJacobians(double t, double[] y, double[] yDot, double[][] dFdY, double[][] dFdP)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 */ StepHandlerWrapper(final StepHandlerWithJacobians handler, final int n, final int k)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 */ getHandler()567 public StepHandlerWithJacobians getHandler() { 568 return handler; 569 } 570 571 /** {@inheritDoc} */ handleStep(StepInterpolator interpolator, boolean isLast)572 public void handleStep(StepInterpolator interpolator, boolean isLast) 573 throws DerivativeException { 574 handler.handleStep(new StepInterpolatorWrapper(interpolator, n, k), isLast); 575 } 576 577 /** {@inheritDoc} */ requiresDenseOutput()578 public boolean requiresDenseOutput() { 579 return handler.requiresDenseOutput(); 580 } 581 582 /** {@inheritDoc} */ reset()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") StepInterpolatorWrapper()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 */ StepInterpolatorWrapper(final StepInterpolator interpolator, final int n, final int k)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} */ setInterpolatedTime(double time)638 public void setInterpolatedTime(double time) { 639 interpolator.setInterpolatedTime(time); 640 } 641 642 /** {@inheritDoc} */ isForward()643 public boolean isForward() { 644 return interpolator.isForward(); 645 } 646 647 /** {@inheritDoc} */ getPreviousTime()648 public double getPreviousTime() { 649 return interpolator.getPreviousTime(); 650 } 651 652 /** {@inheritDoc} */ getInterpolatedTime()653 public double getInterpolatedTime() { 654 return interpolator.getInterpolatedTime(); 655 } 656 657 /** {@inheritDoc} */ getInterpolatedY()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} */ getInterpolatedDyDy0()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} */ getInterpolatedDyDp()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} */ getInterpolatedYDot()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} */ getInterpolatedDyDy0Dot()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} */ getInterpolatedDyDpDot()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} */ getCurrentTime()722 public double getCurrentTime() { 723 return interpolator.getCurrentTime(); 724 } 725 726 /** {@inheritDoc} */ copy()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} */ writeExternal(ObjectOutput out)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} */ readExternal(ObjectInput in)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 */ copyArray(final double[] src, final double[] dest)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 */ copyArray(final double[][] src, final double[][] dest)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 */ writeArray(final ObjectOutput out, final double[] array)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 */ writeArray(final ObjectOutput out, final double[][] array)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 */ readArray(final ObjectInput in, final double[] array)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 */ readArray(final ObjectInput in, final double[][] array)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 */ EventHandlerWrapper(final EventHandlerWithJacobians handler, final int n, final int k)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 */ getHandler()872 public EventHandlerWithJacobians getHandler() { 873 return handler; 874 } 875 876 /** {@inheritDoc} */ eventOccurred(double t, double[] z, boolean increasing)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} */ g(double t, double[] z)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} */ resetState(double t, double[] z)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 } 900