EventState.java

  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.  */

  18. package org.apache.commons.math4.legacy.ode.events;

  20. import org.apache.commons.math4.legacy.analysis.UnivariateFunction;
  21. import org.apache.commons.math4.legacy.analysis.solvers.AllowedSolution;
  22. import org.apache.commons.math4.legacy.analysis.solvers.BracketedUnivariateSolver;
  23. import org.apache.commons.math4.legacy.analysis.solvers.PegasusSolver;
  24. import org.apache.commons.math4.legacy.analysis.solvers.UnivariateSolver;
  25. import org.apache.commons.math4.legacy.analysis.solvers.UnivariateSolverUtils;
  26. import org.apache.commons.math4.legacy.exception.MaxCountExceededException;
  27. import org.apache.commons.math4.legacy.exception.NoBracketingException;
  28. import org.apache.commons.math4.legacy.ode.EquationsMapper;
  29. import org.apache.commons.math4.legacy.ode.ExpandableStatefulODE;
  30. import org.apache.commons.math4.legacy.ode.sampling.StepInterpolator;
  31. import org.apache.commons.math4.core.jdkmath.JdkMath;

  33. /** This class handles the state for one {@link EventHandler
  34.  * event handler} during integration steps.
  35.  *
  36.  * <p>Each time the integrator proposes a step, the event handler
  37.  * switching function should be checked. This class handles the state
  38.  * of one handler during one integration step, with references to the
  39.  * state at the end of the preceding step. This information is used to
  40.  * decide if the handler should trigger an event or not during the
  41.  * proposed step.</p>
  42.  *
  43.  * @since 1.2
  44.  */
  45. public class EventState {

  47.     /** Event handler. */
  48.     private final EventHandler handler;

  50.     /** Maximal time interval between events handler checks. */
  51.     private final double maxCheckInterval;

  53.     /** Convergence threshold for event localization. */
  54.     private final double convergence;

  56.     /** Upper limit in the iteration count for event localization. */
  57.     private final int maxIterationCount;

  59.     /** Equation being integrated. */
  60.     private ExpandableStatefulODE expandable;

  62.     /** Time at the beginning of the step. */
  63.     private double t0;

  65.     /** Value of the events handler at the beginning of the step. */
  66.     private double g0;

  68.     /** Simulated sign of g0 (we cheat when crossing events). */
  69.     private boolean g0Positive;

  71.     /** Indicator of event expected during the step. */
  72.     private boolean pendingEvent;

  74.     /** Occurrence time of the pending event. */
  75.     private double pendingEventTime;

  77.     /** Occurrence time of the previous event. */
  78.     private double previousEventTime;

  80.     /** Integration direction. */
  81.     private boolean forward;

  83.     /** Variation direction around pending event.
  84.      *  (this is considered with respect to the integration direction)
  85.      */
  86.     private boolean increasing;

  88.     /** Next action indicator. */
  89.     private EventHandler.Action nextAction;

  91.     /** Root-finding algorithm to use to detect state events. */
  92.     private final UnivariateSolver solver;

  94.     /** Simple constructor.
  95.      * @param handler event handler
  96.      * @param maxCheckInterval maximal time interval between switching
  97.      * function checks (this interval prevents missing sign changes in
  98.      * case the integration steps becomes very large)
  99.      * @param convergence convergence threshold in the event time search
  100.      * @param maxIterationCount upper limit of the iteration count in
  101.      * the event time search
  102.      * @param solver Root-finding algorithm to use to detect state events
  103.      */
  104.     public EventState(final EventHandler handler, final double maxCheckInterval,
  105.                       final double convergence, final int maxIterationCount,
  106.                       final UnivariateSolver solver) {
  107.         this.handler           = handler;
  108.         this.maxCheckInterval  = maxCheckInterval;
  109.         this.convergence       = JdkMath.abs(convergence);
  110.         this.maxIterationCount = maxIterationCount;
  111.         this.solver            = solver;

  113.         // some dummy values ...
  114.         expandable        = null;
  115.         t0                = Double.NaN;
  116.         g0                = Double.NaN;
  117.         g0Positive        = true;
  118.         pendingEvent      = false;
  119.         pendingEventTime  = Double.NaN;
  120.         previousEventTime = Double.NaN;
  121.         increasing        = true;
  122.         nextAction        = EventHandler.Action.CONTINUE;
  123.     }

  125.     /** Get the underlying event handler.
  126.      * @return underlying event handler
  127.      */
  128.     public EventHandler getEventHandler() {
  129.         return handler;
  130.     }

  132.     /** Set the equation.
  133.      * @param expandable equation being integrated
  134.      */
  135.     public void setExpandable(final ExpandableStatefulODE expandable) {
  136.         this.expandable = expandable;
  137.     }

  139.     /** Get the maximal time interval between events handler checks.
  140.      * @return maximal time interval between events handler checks
  141.      */
  142.     public double getMaxCheckInterval() {
  143.         return maxCheckInterval;
  144.     }

  146.     /** Get the convergence threshold for event localization.
  147.      * @return convergence threshold for event localization
  148.      */
  149.     public double getConvergence() {
  150.         return convergence;
  151.     }

  153.     /** Get the upper limit in the iteration count for event localization.
  154.      * @return upper limit in the iteration count for event localization
  155.      */
  156.     public int getMaxIterationCount() {
  157.         return maxIterationCount;
  158.     }

  160.     /** Reinitialize the beginning of the step.
  161.      * @param interpolator valid for the current step
  162.      * @exception MaxCountExceededException if the interpolator throws one because
  163.      * the number of functions evaluations is exceeded
  164.      */
  165.     public void reinitializeBegin(final StepInterpolator interpolator)
  166.         throws MaxCountExceededException {

  168.         t0 = interpolator.getPreviousTime();
  169.         interpolator.setInterpolatedTime(t0);
  170.         g0 = handler.g(t0, getCompleteState(interpolator));
  171.         if (g0 == 0) {
  172.             // excerpt from MATH-421 issue:
  173.             // If an ODE solver is setup with an EventHandler that return STOP
  174.             // when the even is triggered, the integrator stops (which is exactly
  175.             // the expected behavior). If however the user wants to restart the
  176.             // solver from the final state reached at the event with the same
  177.             // configuration (expecting the event to be triggered again at a
  178.             // later time), then the integrator may fail to start. It can get stuck
  179.             // at the previous event. The use case for the bug MATH-421 is fairly
  180.             // general, so events occurring exactly at start in the first step should
  181.             // be ignored.

  183.             // extremely rare case: there is a zero EXACTLY at interval start
  184.             // we will use the sign slightly after step beginning to force ignoring this zero
  185.             final double epsilon = JdkMath.max(solver.getAbsoluteAccuracy(),
  186.                                                 JdkMath.abs(solver.getRelativeAccuracy() * t0));
  187.             final double tStart = t0 + 0.5 * epsilon;
  188.             interpolator.setInterpolatedTime(tStart);
  189.             g0 = handler.g(tStart, getCompleteState(interpolator));
  190.         }
  191.         g0Positive = g0 >= 0;
  192.     }

  194.     /** Get the complete state (primary and secondary).
  195.      * @param interpolator interpolator to use
  196.      * @return complete state
  197.      */
  198.     private double[] getCompleteState(final StepInterpolator interpolator) {

  200.         final double[] complete = new double[expandable.getTotalDimension()];

  202.         expandable.getPrimaryMapper().insertEquationData(interpolator.getInterpolatedState(),
  203.                                                          complete);
  204.         int index = 0;
  205.         for (EquationsMapper secondary : expandable.getSecondaryMappers()) {
  206.             secondary.insertEquationData(interpolator.getInterpolatedSecondaryState(index++),
  207.                                          complete);
  208.         }

  210.         return complete;
  211.     }

  213.     /** Evaluate the impact of the proposed step on the event handler.
  214.      * @param interpolator step interpolator for the proposed step
  215.      * @return true if the event handler triggers an event before
  216.      * the end of the proposed step
  217.      * @exception MaxCountExceededException if the interpolator throws one because
  218.      * the number of functions evaluations is exceeded
  219.      * @exception NoBracketingException if the event cannot be bracketed
  220.      */
  221.     public boolean evaluateStep(final StepInterpolator interpolator)
  222.         throws MaxCountExceededException, NoBracketingException {

  224.         try {
  225.             forward = interpolator.isForward();
  226.             final double t1 = interpolator.getCurrentTime();
  227.             final double dt = t1 - t0;
  228.             if (JdkMath.abs(dt) < convergence) {
  229.                 // we cannot do anything on such a small step, don't trigger any events
  230.                 return false;
  231.             }
  232.             final int    n = JdkMath.max(1, (int) JdkMath.ceil(JdkMath.abs(dt) / maxCheckInterval));
  233.             final double h = dt / n;

  235.             final UnivariateFunction f = new UnivariateFunction() {
  236.                 /** {@inheritDoc} */
  237.                 @Override
  238.                 public double value(final double t) throws LocalMaxCountExceededException {
  239.                     try {
  240.                         interpolator.setInterpolatedTime(t);
  241.                         return handler.g(t, getCompleteState(interpolator));
  242.                     } catch (MaxCountExceededException mcee) {
  243.                         throw new LocalMaxCountExceededException(mcee);
  244.                     }
  245.                 }
  246.             };

  248.             double ta = t0;
  249.             double ga = g0;
  250.             for (int i = 0; i < n; ++i) {

  252.                 // evaluate handler value at the end of the substep
  253.                 final double tb = (i == n - 1) ? t1 : t0 + (i + 1) * h;
  254.                 interpolator.setInterpolatedTime(tb);
  255.                 final double gb = handler.g(tb, getCompleteState(interpolator));

  257.                 // check events occurrence
  258.                 if (g0Positive ^ (gb >= 0)) {
  259.                     // there is a sign change: an event is expected during this step

  261.                     // variation direction, with respect to the integration direction
  262.                     increasing = gb >= ga;

  264.                     // find the event time making sure we select a solution just at or past the exact root
  265.                     final double root;
  266.                     if (solver instanceof BracketedUnivariateSolver<?>) {
  267.                         @SuppressWarnings("unchecked")
  268.                         BracketedUnivariateSolver<UnivariateFunction> bracketing =
  269.                                 (BracketedUnivariateSolver<UnivariateFunction>) solver;
  270.                         root = forward ?
  271.                                bracketing.solve(maxIterationCount, f, ta, tb, AllowedSolution.RIGHT_SIDE) :
  272.                                bracketing.solve(maxIterationCount, f, tb, ta, AllowedSolution.LEFT_SIDE);
  273.                     } else {
  274.                         final double baseRoot = forward ?
  275.                                                 solver.solve(maxIterationCount, f, ta, tb) :
  276.                                                 solver.solve(maxIterationCount, f, tb, ta);
  277.                         final int remainingEval = maxIterationCount - solver.getEvaluations();
  278.                         BracketedUnivariateSolver<UnivariateFunction> bracketing =
  279.                                 new PegasusSolver(solver.getRelativeAccuracy(), solver.getAbsoluteAccuracy());
  280.                         root = forward ?
  281.                                UnivariateSolverUtils.forceSide(remainingEval, f, bracketing,
  282.                                                                    baseRoot, ta, tb, AllowedSolution.RIGHT_SIDE) :
  283.                                UnivariateSolverUtils.forceSide(remainingEval, f, bracketing,
  284.                                                                    baseRoot, tb, ta, AllowedSolution.LEFT_SIDE);
  285.                     }

  287.                     if (!Double.isNaN(previousEventTime) &&
  288.                         JdkMath.abs(root - ta) <= convergence &&
  289.                         JdkMath.abs(root - previousEventTime) <= convergence) {
  290.                         // we have either found nothing or found (again ?) a past event,
  291.                         // retry the substep excluding this value, and taking care to have the
  292.                         // required sign in case the g function is noisy around its zero and
  293.                         // crosses the axis several times
  294.                         do {
  295.                             ta = forward ? ta + convergence : ta - convergence;
  296.                             ga = f.value(ta);
  297.                         } while ((g0Positive ^ (ga >= 0)) && (forward ^ (ta >= tb)));

  299.                         if (forward ^ (ta >= tb)) {
  300.                             // we were able to skip this spurious root
  301.                             --i;
  302.                         } else {
  303.                             // we can't avoid this root before the end of the step,
  304.                             // we have to handle it despite it is close to the former one
  305.                             // maybe we have two very close roots
  306.                             pendingEventTime = root;
  307.                             pendingEvent = true;
  308.                             return true;
  309.                         }
  310.                     } else if (Double.isNaN(previousEventTime) ||
  311.                                JdkMath.abs(previousEventTime - root) > convergence) {
  312.                         pendingEventTime = root;
  313.                         pendingEvent = true;
  314.                         return true;
  315.                     } else {
  316.                         // no sign change: there is no event for now
  317.                         ta = tb;
  318.                         ga = gb;
  319.                     }
  320.                 } else {
  321.                     // no sign change: there is no event for now
  322.                     ta = tb;
  323.                     ga = gb;
  324.                 }
  325.             }

  327.             // no event during the whole step
  328.             pendingEvent     = false;
  329.             pendingEventTime = Double.NaN;
  330.             return false;
  331.         } catch (LocalMaxCountExceededException lmcee) {
  332.             throw lmcee.getException();
  333.         }
  334.     }

  336.     /** Get the occurrence time of the event triggered in the current step.
  337.      * @return occurrence time of the event triggered in the current
  338.      * step or infinity if no events are triggered
  339.      */
  340.     public double getEventTime() {
  341.         return pendingEvent ?
  342.                pendingEventTime :
  343.                (forward ? Double.POSITIVE_INFINITY : Double.NEGATIVE_INFINITY);
  344.     }

  346.     /** Acknowledge the fact the step has been accepted by the integrator.
  347.      * @param t value of the independent <i>time</i> variable at the
  348.      * end of the step
  349.      * @param y array containing the current value of the state vector
  350.      * at the end of the step
  351.      */
  352.     public void stepAccepted(final double t, final double[] y) {

  354.         t0 = t;
  355.         g0 = handler.g(t, y);

  357.         if (pendingEvent && JdkMath.abs(pendingEventTime - t) <= convergence) {
  358.             // force the sign to its value "just after the event"
  359.             previousEventTime = t;
  360.             g0Positive        = increasing;
  361.             nextAction        = handler.eventOccurred(t, y, increasing == forward);
  362.         } else {
  363.             g0Positive = g0 >= 0;
  364.             nextAction = EventHandler.Action.CONTINUE;
  365.         }
  366.     }

  368.     /** Check if the integration should be stopped at the end of the
  369.      * current step.
  370.      * @return true if the integration should be stopped
  371.      */
  372.     public boolean stop() {
  373.         return nextAction == EventHandler.Action.STOP;
  374.     }

  376.     /** Let the event handler reset the state if it wants.
  377.      * @param t value of the independent <i>time</i> variable at the
  378.      * beginning of the next step
  379.      * @param y array were to put the desired state vector at the beginning
  380.      * of the next step
  381.      * @return true if the integrator should reset the derivatives too
  382.      */
  383.     public boolean reset(final double t, final double[] y) {

  385.         if (!(pendingEvent && JdkMath.abs(pendingEventTime - t) <= convergence)) {
  386.             return false;
  387.         }

  389.         if (nextAction == EventHandler.Action.RESET_STATE) {
  390.             handler.resetState(t, y);
  391.         }
  392.         pendingEvent      = false;
  393.         pendingEventTime  = Double.NaN;

  395.         return nextAction == EventHandler.Action.RESET_STATE ||
  396.                nextAction == EventHandler.Action.RESET_DERIVATIVES;
  397.     }

  399.     /** Local wrapper to propagate exceptions. */
  400.     private static final class LocalMaxCountExceededException extends RuntimeException {

  402.         /** Serializable UID. */
  403.         private static final long serialVersionUID = 20120901L;

  405.         /** Wrapped exception. */
  406.         private final MaxCountExceededException wrapped;

  408.         /** Simple constructor.
  409.          * @param exception exception to wrap
  410.          */
  411.         LocalMaxCountExceededException(final MaxCountExceededException exception) {
  412.             wrapped = exception;
  413.         }

  415.         /** Get the wrapped exception.
  416.          * @return wrapped exception
  417.          */
  418.         public MaxCountExceededException getException() {
  419.             return wrapped;
  420.         }
  421.     }
  422. }
Created with JaCoCo 0.8.10.202304240956Code Coverage Report for Apache Commons Math 4.0-SNAPSHOT