source: src/main/java/agents/anac/y2019/harddealer/math3/ode/events/FieldEventState.java

Last change on this file was 204, checked in by Katsuhide Fujita, 5 years ago

Fixed errors of ANAC2019 agents
  • Property svn:executable set to *
File size: 13.9 KB
Line 
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
18package agents.anac.y2019.harddealer.math3.ode.events;
19
20import agents.anac.y2019.harddealer.math3.RealFieldElement;
21import agents.anac.y2019.harddealer.math3.analysis.RealFieldUnivariateFunction;
22import agents.anac.y2019.harddealer.math3.analysis.solvers.AllowedSolution;
23import agents.anac.y2019.harddealer.math3.analysis.solvers.BracketedRealFieldUnivariateSolver;
24import agents.anac.y2019.harddealer.math3.exception.MaxCountExceededException;
25import agents.anac.y2019.harddealer.math3.exception.NoBracketingException;
26import agents.anac.y2019.harddealer.math3.ode.FieldODEState;
27import agents.anac.y2019.harddealer.math3.ode.FieldODEStateAndDerivative;
28import agents.anac.y2019.harddealer.math3.ode.sampling.FieldStepInterpolator;
29import agents.anac.y2019.harddealer.math3.util.FastMath;
30
31/** This class handles the state for one {@link EventHandler
32 * event handler} during integration steps.
33 *
34 * <p>Each time the integrator proposes a step, the event handler
35 * switching function should be checked. This class handles the state
36 * of one handler during one integration step, with references to the
37 * state at the end of the preceding step. This information is used to
38 * decide if the handler should trigger an event or not during the
39 * proposed step.</p>
40 *
41 * @param <T> the type of the field elements
42 * @since 3.6
43 */
44public class FieldEventState<T extends RealFieldElement<T>> {
45
46 /** Event handler. */
47 private final FieldEventHandler<T> handler;
48
49 /** Maximal time interval between events handler checks. */
50 private final double maxCheckInterval;
51
52 /** Convergence threshold for event localization. */
53 private final T convergence;
54
55 /** Upper limit in the iteration count for event localization. */
56 private final int maxIterationCount;
57
58 /** Time at the beginning of the step. */
59 private T t0;
60
61 /** Value of the events handler at the beginning of the step. */
62 private T g0;
63
64 /** Simulated sign of g0 (we cheat when crossing events). */
65 private boolean g0Positive;
66
67 /** Indicator of event expected during the step. */
68 private boolean pendingEvent;
69
70 /** Occurrence time of the pending event. */
71 private T pendingEventTime;
72
73 /** Occurrence time of the previous event. */
74 private T previousEventTime;
75
76 /** Integration direction. */
77 private boolean forward;
78
79 /** Variation direction around pending event.
80 * (this is considered with respect to the integration direction)
81 */
82 private boolean increasing;
83
84 /** Next action indicator. */
85 private Action nextAction;
86
87 /** Root-finding algorithm to use to detect state events. */
88 private final BracketedRealFieldUnivariateSolver<T> solver;
89
90 /** Simple constructor.
91 * @param handler event handler
92 * @param maxCheckInterval maximal time interval between switching
93 * function checks (this interval prevents missing sign changes in
94 * case the integration steps becomes very large)
95 * @param convergence convergence threshold in the event time search
96 * @param maxIterationCount upper limit of the iteration count in
97 * the event time search
98 * @param solver Root-finding algorithm to use to detect state events
99 */
100 public FieldEventState(final FieldEventHandler<T> handler, final double maxCheckInterval,
101 final T convergence, final int maxIterationCount,
102 final BracketedRealFieldUnivariateSolver<T> solver) {
103 this.handler = handler;
104 this.maxCheckInterval = maxCheckInterval;
105 this.convergence = convergence.abs();
106 this.maxIterationCount = maxIterationCount;
107 this.solver = solver;
108
109 // some dummy values ...
110 t0 = null;
111 g0 = null;
112 g0Positive = true;
113 pendingEvent = false;
114 pendingEventTime = null;
115 previousEventTime = null;
116 increasing = true;
117 nextAction = Action.CONTINUE;
118
119 }
120
121 /** Get the underlying event handler.
122 * @return underlying event handler
123 */
124 public FieldEventHandler<T> getEventHandler() {
125 return handler;
126 }
127
128 /** Get the maximal time interval between events handler checks.
129 * @return maximal time interval between events handler checks
130 */
131 public double getMaxCheckInterval() {
132 return maxCheckInterval;
133 }
134
135 /** Get the convergence threshold for event localization.
136 * @return convergence threshold for event localization
137 */
138 public T getConvergence() {
139 return convergence;
140 }
141
142 /** Get the upper limit in the iteration count for event localization.
143 * @return upper limit in the iteration count for event localization
144 */
145 public int getMaxIterationCount() {
146 return maxIterationCount;
147 }
148
149 /** Reinitialize the beginning of the step.
150 * @param interpolator valid for the current step
151 * @exception MaxCountExceededException if the interpolator throws one because
152 * the number of functions evaluations is exceeded
153 */
154 public void reinitializeBegin(final FieldStepInterpolator<T> interpolator)
155 throws MaxCountExceededException {
156
157 final FieldODEStateAndDerivative<T> s0 = interpolator.getPreviousState();
158 t0 = s0.getTime();
159 g0 = handler.g(s0);
160 if (g0.getReal() == 0) {
161 // excerpt from MATH-421 issue:
162 // If an ODE solver is setup with an EventHandler that return STOP
163 // when the even is triggered, the integrator stops (which is exactly
164 // the expected behavior). If however the user wants to restart the
165 // solver from the final state reached at the event with the same
166 // configuration (expecting the event to be triggered again at a
167 // later time), then the integrator may fail to start. It can get stuck
168 // at the previous event. The use case for the bug MATH-421 is fairly
169 // general, so events occurring exactly at start in the first step should
170 // be ignored.
171
172 // extremely rare case: there is a zero EXACTLY at interval start
173 // we will use the sign slightly after step beginning to force ignoring this zero
174 final double epsilon = FastMath.max(solver.getAbsoluteAccuracy().getReal(),
175 FastMath.abs(solver.getRelativeAccuracy().multiply(t0).getReal()));
176 final T tStart = t0.add(0.5 * epsilon);
177 g0 = handler.g(interpolator.getInterpolatedState(tStart));
178 }
179 g0Positive = g0.getReal() >= 0;
180
181 }
182
183 /** Evaluate the impact of the proposed step on the event handler.
184 * @param interpolator step interpolator for the proposed step
185 * @return true if the event handler triggers an event before
186 * the end of the proposed step
187 * @exception MaxCountExceededException if the interpolator throws one because
188 * the number of functions evaluations is exceeded
189 * @exception NoBracketingException if the event cannot be bracketed
190 */
191 public boolean evaluateStep(final FieldStepInterpolator<T> interpolator)
192 throws MaxCountExceededException, NoBracketingException {
193
194 forward = interpolator.isForward();
195 final FieldODEStateAndDerivative<T> s1 = interpolator.getCurrentState();
196 final T t1 = s1.getTime();
197 final T dt = t1.subtract(t0);
198 if (dt.abs().subtract(convergence).getReal() < 0) {
199 // we cannot do anything on such a small step, don't trigger any events
200 return false;
201 }
202 final int n = FastMath.max(1, (int) FastMath.ceil(FastMath.abs(dt.getReal()) / maxCheckInterval));
203 final T h = dt.divide(n);
204
205 final RealFieldUnivariateFunction<T> f = new RealFieldUnivariateFunction<T>() {
206 /** {@inheritDoc} */
207 public T value(final T t) {
208 return handler.g(interpolator.getInterpolatedState(t));
209 }
210 };
211
212 T ta = t0;
213 T ga = g0;
214 for (int i = 0; i < n; ++i) {
215
216 // evaluate handler value at the end of the substep
217 final T tb = (i == n - 1) ? t1 : t0.add(h.multiply(i + 1));
218 final T gb = handler.g(interpolator.getInterpolatedState(tb));
219
220 // check events occurrence
221 if (g0Positive ^ (gb.getReal() >= 0)) {
222 // there is a sign change: an event is expected during this step
223
224 // variation direction, with respect to the integration direction
225 increasing = gb.subtract(ga).getReal() >= 0;
226
227 // find the event time making sure we select a solution just at or past the exact root
228 final T root = forward ?
229 solver.solve(maxIterationCount, f, ta, tb, AllowedSolution.RIGHT_SIDE) :
230 solver.solve(maxIterationCount, f, tb, ta, AllowedSolution.LEFT_SIDE);
231
232 if (previousEventTime != null &&
233 root.subtract(ta).abs().subtract(convergence).getReal() <= 0 &&
234 root.subtract(previousEventTime).abs().subtract(convergence).getReal() <= 0) {
235 // we have either found nothing or found (again ?) a past event,
236 // retry the substep excluding this value, and taking care to have the
237 // required sign in case the g function is noisy around its zero and
238 // crosses the axis several times
239 do {
240 ta = forward ? ta.add(convergence) : ta.subtract(convergence);
241 ga = f.value(ta);
242 } while ((g0Positive ^ (ga.getReal() >= 0)) && (forward ^ (ta.subtract(tb).getReal() >= 0)));
243
244 if (forward ^ (ta.subtract(tb).getReal() >= 0)) {
245 // we were able to skip this spurious root
246 --i;
247 } else {
248 // we can't avoid this root before the end of the step,
249 // we have to handle it despite it is close to the former one
250 // maybe we have two very close roots
251 pendingEventTime = root;
252 pendingEvent = true;
253 return true;
254 }
255 } else if (previousEventTime == null ||
256 previousEventTime.subtract(root).abs().subtract(convergence).getReal() > 0) {
257 pendingEventTime = root;
258 pendingEvent = true;
259 return true;
260 } else {
261 // no sign change: there is no event for now
262 ta = tb;
263 ga = gb;
264 }
265
266 } else {
267 // no sign change: there is no event for now
268 ta = tb;
269 ga = gb;
270 }
271
272 }
273
274 // no event during the whole step
275 pendingEvent = false;
276 pendingEventTime = null;
277 return false;
278
279 }
280
281 /** Get the occurrence time of the event triggered in the current step.
282 * @return occurrence time of the event triggered in the current
283 * step or infinity if no events are triggered
284 */
285 public T getEventTime() {
286 return pendingEvent ?
287 pendingEventTime :
288 t0.getField().getZero().add(forward ? Double.POSITIVE_INFINITY : Double.NEGATIVE_INFINITY);
289 }
290
291 /** Acknowledge the fact the step has been accepted by the integrator.
292 * @param state state at the end of the step
293 */
294 public void stepAccepted(final FieldODEStateAndDerivative<T> state) {
295
296 t0 = state.getTime();
297 g0 = handler.g(state);
298
299 if (pendingEvent && pendingEventTime.subtract(state.getTime()).abs().subtract(convergence).getReal() <= 0) {
300 // force the sign to its value "just after the event"
301 previousEventTime = state.getTime();
302 g0Positive = increasing;
303 nextAction = handler.eventOccurred(state, !(increasing ^ forward));
304 } else {
305 g0Positive = g0.getReal() >= 0;
306 nextAction = Action.CONTINUE;
307 }
308 }
309
310 /** Check if the integration should be stopped at the end of the
311 * current step.
312 * @return true if the integration should be stopped
313 */
314 public boolean stop() {
315 return nextAction == Action.STOP;
316 }
317
318 /** Let the event handler reset the state if it wants.
319 * @param state state at the beginning of the next step
320 * @return reset state (may by the same as initial state if only
321 * derivatives should be reset), or null if nothing is reset
322 */
323 public FieldODEState<T> reset(final FieldODEStateAndDerivative<T> state) {
324
325 if (!(pendingEvent && pendingEventTime.subtract(state.getTime()).abs().subtract(convergence).getReal() <= 0)) {
326 return null;
327 }
328
329 final FieldODEState<T> newState;
330 if (nextAction == Action.RESET_STATE) {
331 newState = handler.resetState(state);
332 } else if (nextAction == Action.RESET_DERIVATIVES) {
333 newState = state;
334 } else {
335 newState = null;
336 }
337 pendingEvent = false;
338 pendingEventTime = null;
339
340 return newState;
341
342 }
343
344}
Note: See TracBrowser for help on using the repository browser.