FieldEventState.java
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math4.legacy.ode.events;
import org.apache.commons.math4.legacy.core.RealFieldElement;
import org.apache.commons.math4.legacy.analysis.RealFieldUnivariateFunction;
import org.apache.commons.math4.legacy.analysis.solvers.AllowedSolution;
import org.apache.commons.math4.legacy.analysis.solvers.BracketedRealFieldUnivariateSolver;
import org.apache.commons.math4.legacy.exception.MaxCountExceededException;
import org.apache.commons.math4.legacy.exception.NoBracketingException;
import org.apache.commons.math4.legacy.ode.FieldODEState;
import org.apache.commons.math4.legacy.ode.FieldODEStateAndDerivative;
import org.apache.commons.math4.legacy.ode.sampling.FieldStepInterpolator;
import org.apache.commons.math4.core.jdkmath.JdkMath;
/** This class handles the state for one {@link EventHandler
* event handler} during integration steps.
*
* <p>Each time the integrator proposes a step, the event handler
* switching function should be checked. This class handles the state
* of one handler during one integration step, with references to the
* state at the end of the preceding step. This information is used to
* decide if the handler should trigger an event or not during the
* proposed step.</p>
*
* @param <T> the type of the field elements
* @since 3.6
*/
public class FieldEventState<T extends RealFieldElement<T>> {
/** Event handler. */
private final FieldEventHandler<T> handler;
/** Maximal time interval between events handler checks. */
private final double maxCheckInterval;
/** Convergence threshold for event localization. */
private final T convergence;
/** Upper limit in the iteration count for event localization. */
private final int maxIterationCount;
/** Time at the beginning of the step. */
private T t0;
/** Value of the events handler at the beginning of the step. */
private T g0;
/** Simulated sign of g0 (we cheat when crossing events). */
private boolean g0Positive;
/** Indicator of event expected during the step. */
private boolean pendingEvent;
/** Occurrence time of the pending event. */
private T pendingEventTime;
/** Occurrence time of the previous event. */
private T previousEventTime;
/** Integration direction. */
private boolean forward;
/** Variation direction around pending event.
* (this is considered with respect to the integration direction)
*/
private boolean increasing;
/** Next action indicator. */
private Action nextAction;
/** Root-finding algorithm to use to detect state events. */
private final BracketedRealFieldUnivariateSolver<T> solver;
/** Simple constructor.
* @param handler event handler
* @param maxCheckInterval maximal time interval between switching
* function checks (this interval prevents missing sign changes in
* case the integration steps becomes very large)
* @param convergence convergence threshold in the event time search
* @param maxIterationCount upper limit of the iteration count in
* the event time search
* @param solver Root-finding algorithm to use to detect state events
*/
public FieldEventState(final FieldEventHandler<T> handler, final double maxCheckInterval,
final T convergence, final int maxIterationCount,
final BracketedRealFieldUnivariateSolver<T> solver) {
this.handler = handler;
this.maxCheckInterval = maxCheckInterval;
this.convergence = convergence.abs();
this.maxIterationCount = maxIterationCount;
this.solver = solver;
// some dummy values ...
t0 = null;
g0 = null;
g0Positive = true;
pendingEvent = false;
pendingEventTime = null;
previousEventTime = null;
increasing = true;
nextAction = Action.CONTINUE;
}
/** Get the underlying event handler.
* @return underlying event handler
*/
public FieldEventHandler<T> getEventHandler() {
return handler;
}
/** Get the maximal time interval between events handler checks.
* @return maximal time interval between events handler checks
*/
public double getMaxCheckInterval() {
return maxCheckInterval;
}
/** Get the convergence threshold for event localization.
* @return convergence threshold for event localization
*/
public T getConvergence() {
return convergence;
}
/** Get the upper limit in the iteration count for event localization.
* @return upper limit in the iteration count for event localization
*/
public int getMaxIterationCount() {
return maxIterationCount;
}
/** Reinitialize the beginning of the step.
* @param interpolator valid for the current step
* @exception MaxCountExceededException if the interpolator throws one because
* the number of functions evaluations is exceeded
*/
public void reinitializeBegin(final FieldStepInterpolator<T> interpolator)
throws MaxCountExceededException {
final FieldODEStateAndDerivative<T> s0 = interpolator.getPreviousState();
t0 = s0.getTime();
g0 = handler.g(s0);
if (g0.getReal() == 0) {
// excerpt from MATH-421 issue:
// If an ODE solver is setup with an EventHandler that return STOP
// when the even is triggered, the integrator stops (which is exactly
// the expected behavior). If however the user wants to restart the
// solver from the final state reached at the event with the same
// configuration (expecting the event to be triggered again at a
// later time), then the integrator may fail to start. It can get stuck
// at the previous event. The use case for the bug MATH-421 is fairly
// general, so events occurring exactly at start in the first step should
// be ignored.
// extremely rare case: there is a zero EXACTLY at interval start
// we will use the sign slightly after step beginning to force ignoring this zero
final double epsilon = JdkMath.max(solver.getAbsoluteAccuracy().getReal(),
JdkMath.abs(solver.getRelativeAccuracy().multiply(t0).getReal()));
final T tStart = t0.add(0.5 * epsilon);
g0 = handler.g(interpolator.getInterpolatedState(tStart));
}
g0Positive = g0.getReal() >= 0;
}
/** Evaluate the impact of the proposed step on the event handler.
* @param interpolator step interpolator for the proposed step
* @return true if the event handler triggers an event before
* the end of the proposed step
* @exception MaxCountExceededException if the interpolator throws one because
* the number of functions evaluations is exceeded
* @exception NoBracketingException if the event cannot be bracketed
*/
public boolean evaluateStep(final FieldStepInterpolator<T> interpolator)
throws MaxCountExceededException, NoBracketingException {
forward = interpolator.isForward();
final FieldODEStateAndDerivative<T> s1 = interpolator.getCurrentState();
final T t1 = s1.getTime();
final T dt = t1.subtract(t0);
if (dt.abs().subtract(convergence).getReal() < 0) {
// we cannot do anything on such a small step, don't trigger any events
return false;
}
final int n = JdkMath.max(1, (int) JdkMath.ceil(JdkMath.abs(dt.getReal()) / maxCheckInterval));
final T h = dt.divide(n);
final RealFieldUnivariateFunction<T> f = new RealFieldUnivariateFunction<T>() {
/** {@inheritDoc} */
@Override
public T value(final T t) {
return handler.g(interpolator.getInterpolatedState(t));
}
};
T ta = t0;
T ga = g0;
for (int i = 0; i < n; ++i) {
// evaluate handler value at the end of the substep
final T tb = (i == n - 1) ? t1 : t0.add(h.multiply(i + 1));
final T gb = handler.g(interpolator.getInterpolatedState(tb));
// check events occurrence
if (g0Positive ^ (gb.getReal() >= 0)) {
// there is a sign change: an event is expected during this step
// variation direction, with respect to the integration direction
increasing = gb.subtract(ga).getReal() >= 0;
// find the event time making sure we select a solution just at or past the exact root
final T root = forward ?
solver.solve(maxIterationCount, f, ta, tb, AllowedSolution.RIGHT_SIDE) :
solver.solve(maxIterationCount, f, tb, ta, AllowedSolution.LEFT_SIDE);
if (previousEventTime != null &&
root.subtract(ta).abs().subtract(convergence).getReal() <= 0 &&
root.subtract(previousEventTime).abs().subtract(convergence).getReal() <= 0) {
// we have either found nothing or found (again ?) a past event,
// retry the substep excluding this value, and taking care to have the
// required sign in case the g function is noisy around its zero and
// crosses the axis several times
do {
ta = forward ? ta.add(convergence) : ta.subtract(convergence);
ga = f.value(ta);
} while ((g0Positive ^ (ga.getReal() >= 0)) && (forward ^ (ta.subtract(tb).getReal() >= 0)));
if (forward ^ (ta.subtract(tb).getReal() >= 0)) {
// we were able to skip this spurious root
--i;
} else {
// we can't avoid this root before the end of the step,
// we have to handle it despite it is close to the former one
// maybe we have two very close roots
pendingEventTime = root;
pendingEvent = true;
return true;
}
} else if (previousEventTime == null ||
previousEventTime.subtract(root).abs().subtract(convergence).getReal() > 0) {
pendingEventTime = root;
pendingEvent = true;
return true;
} else {
// no sign change: there is no event for now
ta = tb;
ga = gb;
}
} else {
// no sign change: there is no event for now
ta = tb;
ga = gb;
}
}
// no event during the whole step
pendingEvent = false;
pendingEventTime = null;
return false;
}
/** Get the occurrence time of the event triggered in the current step.
* @return occurrence time of the event triggered in the current
* step or infinity if no events are triggered
*/
public T getEventTime() {
return pendingEvent ?
pendingEventTime :
t0.getField().getZero().add(forward ? Double.POSITIVE_INFINITY : Double.NEGATIVE_INFINITY);
}
/** Acknowledge the fact the step has been accepted by the integrator.
* @param state state at the end of the step
*/
public void stepAccepted(final FieldODEStateAndDerivative<T> state) {
t0 = state.getTime();
g0 = handler.g(state);
if (pendingEvent && pendingEventTime.subtract(state.getTime()).abs().subtract(convergence).getReal() <= 0) {
// force the sign to its value "just after the event"
previousEventTime = state.getTime();
g0Positive = increasing;
nextAction = handler.eventOccurred(state, increasing == forward);
} else {
g0Positive = g0.getReal() >= 0;
nextAction = Action.CONTINUE;
}
}
/** Check if the integration should be stopped at the end of the
* current step.
* @return true if the integration should be stopped
*/
public boolean stop() {
return nextAction == Action.STOP;
}
/** Let the event handler reset the state if it wants.
* @param state state at the beginning of the next step
* @return reset state (may by the same as initial state if only
* derivatives should be reset), or null if nothing is reset
*/
public FieldODEState<T> reset(final FieldODEStateAndDerivative<T> state) {
if (!(pendingEvent && pendingEventTime.subtract(state.getTime()).abs().subtract(convergence).getReal() <= 0)) {
return null;
}
final FieldODEState<T> newState;
if (nextAction == Action.RESET_STATE) {
newState = handler.resetState(state);
} else if (nextAction == Action.RESET_DERIVATIVES) {
newState = state;
} else {
newState = null;
}
pendingEvent = false;
pendingEventTime = null;
return newState;
}
}