source: src/main/java/agents/anac/y2019/harddealer/math3/ode/events/EventState.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: 17.6 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.analysis.UnivariateFunction;
21import agents.anac.y2019.harddealer.math3.analysis.solvers.AllowedSolution;
22import agents.anac.y2019.harddealer.math3.analysis.solvers.BracketedUnivariateSolver;
23import agents.anac.y2019.harddealer.math3.analysis.solvers.PegasusSolver;
24import agents.anac.y2019.harddealer.math3.analysis.solvers.UnivariateSolver;
25import agents.anac.y2019.harddealer.math3.analysis.solvers.UnivariateSolverUtils;
26import agents.anac.y2019.harddealer.math3.exception.MaxCountExceededException;
27import agents.anac.y2019.harddealer.math3.exception.NoBracketingException;
28import agents.anac.y2019.harddealer.math3.ode.EquationsMapper;
29import agents.anac.y2019.harddealer.math3.ode.ExpandableStatefulODE;
30import agents.anac.y2019.harddealer.math3.ode.sampling.StepInterpolator;
31import agents.anac.y2019.harddealer.math3.util.FastMath;
32
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 */
45public class EventState {
46
47 /** Event handler. */
48 private final EventHandler handler;
49
50 /** Maximal time interval between events handler checks. */
51 private final double maxCheckInterval;
52
53 /** Convergence threshold for event localization. */
54 private final double convergence;
55
56 /** Upper limit in the iteration count for event localization. */
57 private final int maxIterationCount;
58
59 /** Equation being integrated. */
60 private ExpandableStatefulODE expandable;
61
62 /** Time at the beginning of the step. */
63 private double t0;
64
65 /** Value of the events handler at the beginning of the step. */
66 private double g0;
67
68 /** Simulated sign of g0 (we cheat when crossing events). */
69 private boolean g0Positive;
70
71 /** Indicator of event expected during the step. */
72 private boolean pendingEvent;
73
74 /** Occurrence time of the pending event. */
75 private double pendingEventTime;
76
77 /** Occurrence time of the previous event. */
78 private double previousEventTime;
79
80 /** Integration direction. */
81 private boolean forward;
82
83 /** Variation direction around pending event.
84 * (this is considered with respect to the integration direction)
85 */
86 private boolean increasing;
87
88 /** Next action indicator. */
89 private EventHandler.Action nextAction;
90
91 /** Root-finding algorithm to use to detect state events. */
92 private final UnivariateSolver solver;
93
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 = FastMath.abs(convergence);
110 this.maxIterationCount = maxIterationCount;
111 this.solver = solver;
112
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
124 }
125
126 /** Get the underlying event handler.
127 * @return underlying event handler
128 */
129 public EventHandler getEventHandler() {
130 return handler;
131 }
132
133 /** Set the equation.
134 * @param expandable equation being integrated
135 */
136 public void setExpandable(final ExpandableStatefulODE expandable) {
137 this.expandable = expandable;
138 }
139
140 /** Get the maximal time interval between events handler checks.
141 * @return maximal time interval between events handler checks
142 */
143 public double getMaxCheckInterval() {
144 return maxCheckInterval;
145 }
146
147 /** Get the convergence threshold for event localization.
148 * @return convergence threshold for event localization
149 */
150 public double getConvergence() {
151 return convergence;
152 }
153
154 /** Get the upper limit in the iteration count for event localization.
155 * @return upper limit in the iteration count for event localization
156 */
157 public int getMaxIterationCount() {
158 return maxIterationCount;
159 }
160
161 /** Reinitialize the beginning of the step.
162 * @param interpolator valid for the current step
163 * @exception MaxCountExceededException if the interpolator throws one because
164 * the number of functions evaluations is exceeded
165 */
166 public void reinitializeBegin(final StepInterpolator interpolator)
167 throws MaxCountExceededException {
168
169 t0 = interpolator.getPreviousTime();
170 interpolator.setInterpolatedTime(t0);
171 g0 = handler.g(t0, getCompleteState(interpolator));
172 if (g0 == 0) {
173 // excerpt from MATH-421 issue:
174 // If an ODE solver is setup with an EventHandler that return STOP
175 // when the even is triggered, the integrator stops (which is exactly
176 // the expected behavior). If however the user wants to restart the
177 // solver from the final state reached at the event with the same
178 // configuration (expecting the event to be triggered again at a
179 // later time), then the integrator may fail to start. It can get stuck
180 // at the previous event. The use case for the bug MATH-421 is fairly
181 // general, so events occurring exactly at start in the first step should
182 // be ignored.
183
184 // extremely rare case: there is a zero EXACTLY at interval start
185 // we will use the sign slightly after step beginning to force ignoring this zero
186 final double epsilon = FastMath.max(solver.getAbsoluteAccuracy(),
187 FastMath.abs(solver.getRelativeAccuracy() * t0));
188 final double tStart = t0 + 0.5 * epsilon;
189 interpolator.setInterpolatedTime(tStart);
190 g0 = handler.g(tStart, getCompleteState(interpolator));
191 }
192 g0Positive = g0 >= 0;
193
194 }
195
196 /** Get the complete state (primary and secondary).
197 * @param interpolator interpolator to use
198 * @return complete state
199 */
200 private double[] getCompleteState(final StepInterpolator interpolator) {
201
202 final double[] complete = new double[expandable.getTotalDimension()];
203
204 expandable.getPrimaryMapper().insertEquationData(interpolator.getInterpolatedState(),
205 complete);
206 int index = 0;
207 for (EquationsMapper secondary : expandable.getSecondaryMappers()) {
208 secondary.insertEquationData(interpolator.getInterpolatedSecondaryState(index++),
209 complete);
210 }
211
212 return complete;
213
214 }
215
216 /** Evaluate the impact of the proposed step on the event handler.
217 * @param interpolator step interpolator for the proposed step
218 * @return true if the event handler triggers an event before
219 * the end of the proposed step
220 * @exception MaxCountExceededException if the interpolator throws one because
221 * the number of functions evaluations is exceeded
222 * @exception NoBracketingException if the event cannot be bracketed
223 */
224 public boolean evaluateStep(final StepInterpolator interpolator)
225 throws MaxCountExceededException, NoBracketingException {
226
227 try {
228 forward = interpolator.isForward();
229 final double t1 = interpolator.getCurrentTime();
230 final double dt = t1 - t0;
231 if (FastMath.abs(dt) < convergence) {
232 // we cannot do anything on such a small step, don't trigger any events
233 return false;
234 }
235 final int n = FastMath.max(1, (int) FastMath.ceil(FastMath.abs(dt) / maxCheckInterval));
236 final double h = dt / n;
237
238 final UnivariateFunction f = new UnivariateFunction() {
239 /** {@inheritDoc} */
240 public double value(final double t) throws LocalMaxCountExceededException {
241 try {
242 interpolator.setInterpolatedTime(t);
243 return handler.g(t, getCompleteState(interpolator));
244 } catch (MaxCountExceededException mcee) {
245 throw new LocalMaxCountExceededException(mcee);
246 }
247 }
248 };
249
250 double ta = t0;
251 double ga = g0;
252 for (int i = 0; i < n; ++i) {
253
254 // evaluate handler value at the end of the substep
255 final double tb = (i == n - 1) ? t1 : t0 + (i + 1) * h;
256 interpolator.setInterpolatedTime(tb);
257 final double gb = handler.g(tb, getCompleteState(interpolator));
258
259 // check events occurrence
260 if (g0Positive ^ (gb >= 0)) {
261 // there is a sign change: an event is expected during this step
262
263 // variation direction, with respect to the integration direction
264 increasing = gb >= ga;
265
266 // find the event time making sure we select a solution just at or past the exact root
267 final double root;
268 if (solver instanceof BracketedUnivariateSolver<?>) {
269 @SuppressWarnings("unchecked")
270 BracketedUnivariateSolver<UnivariateFunction> bracketing =
271 (BracketedUnivariateSolver<UnivariateFunction>) solver;
272 root = forward ?
273 bracketing.solve(maxIterationCount, f, ta, tb, AllowedSolution.RIGHT_SIDE) :
274 bracketing.solve(maxIterationCount, f, tb, ta, AllowedSolution.LEFT_SIDE);
275 } else {
276 final double baseRoot = forward ?
277 solver.solve(maxIterationCount, f, ta, tb) :
278 solver.solve(maxIterationCount, f, tb, ta);
279 final int remainingEval = maxIterationCount - solver.getEvaluations();
280 BracketedUnivariateSolver<UnivariateFunction> bracketing =
281 new PegasusSolver(solver.getRelativeAccuracy(), solver.getAbsoluteAccuracy());
282 root = forward ?
283 UnivariateSolverUtils.forceSide(remainingEval, f, bracketing,
284 baseRoot, ta, tb, AllowedSolution.RIGHT_SIDE) :
285 UnivariateSolverUtils.forceSide(remainingEval, f, bracketing,
286 baseRoot, tb, ta, AllowedSolution.LEFT_SIDE);
287 }
288
289 if ((!Double.isNaN(previousEventTime)) &&
290 (FastMath.abs(root - ta) <= convergence) &&
291 (FastMath.abs(root - previousEventTime) <= convergence)) {
292 // we have either found nothing or found (again ?) a past event,
293 // retry the substep excluding this value, and taking care to have the
294 // required sign in case the g function is noisy around its zero and
295 // crosses the axis several times
296 do {
297 ta = forward ? ta + convergence : ta - convergence;
298 ga = f.value(ta);
299 } while ((g0Positive ^ (ga >= 0)) && (forward ^ (ta >= tb)));
300
301 if (forward ^ (ta >= tb)) {
302 // we were able to skip this spurious root
303 --i;
304 } else {
305 // we can't avoid this root before the end of the step,
306 // we have to handle it despite it is close to the former one
307 // maybe we have two very close roots
308 pendingEventTime = root;
309 pendingEvent = true;
310 return true;
311 }
312 } else if (Double.isNaN(previousEventTime) ||
313 (FastMath.abs(previousEventTime - root) > convergence)) {
314 pendingEventTime = root;
315 pendingEvent = true;
316 return true;
317 } else {
318 // no sign change: there is no event for now
319 ta = tb;
320 ga = gb;
321 }
322
323 } else {
324 // no sign change: there is no event for now
325 ta = tb;
326 ga = gb;
327 }
328
329 }
330
331 // no event during the whole step
332 pendingEvent = false;
333 pendingEventTime = Double.NaN;
334 return false;
335
336 } catch (LocalMaxCountExceededException lmcee) {
337 throw lmcee.getException();
338 }
339
340 }
341
342 /** Get the occurrence time of the event triggered in the current step.
343 * @return occurrence time of the event triggered in the current
344 * step or infinity if no events are triggered
345 */
346 public double getEventTime() {
347 return pendingEvent ?
348 pendingEventTime :
349 (forward ? Double.POSITIVE_INFINITY : Double.NEGATIVE_INFINITY);
350 }
351
352 /** Acknowledge the fact the step has been accepted by the integrator.
353 * @param t value of the independent <i>time</i> variable at the
354 * end of the step
355 * @param y array containing the current value of the state vector
356 * at the end of the step
357 */
358 public void stepAccepted(final double t, final double[] y) {
359
360 t0 = t;
361 g0 = handler.g(t, y);
362
363 if (pendingEvent && (FastMath.abs(pendingEventTime - t) <= convergence)) {
364 // force the sign to its value "just after the event"
365 previousEventTime = t;
366 g0Positive = increasing;
367 nextAction = handler.eventOccurred(t, y, !(increasing ^ forward));
368 } else {
369 g0Positive = g0 >= 0;
370 nextAction = EventHandler.Action.CONTINUE;
371 }
372 }
373
374 /** Check if the integration should be stopped at the end of the
375 * current step.
376 * @return true if the integration should be stopped
377 */
378 public boolean stop() {
379 return nextAction == EventHandler.Action.STOP;
380 }
381
382 /** Let the event handler reset the state if it wants.
383 * @param t value of the independent <i>time</i> variable at the
384 * beginning of the next step
385 * @param y array were to put the desired state vector at the beginning
386 * of the next step
387 * @return true if the integrator should reset the derivatives too
388 */
389 public boolean reset(final double t, final double[] y) {
390
391 if (!(pendingEvent && (FastMath.abs(pendingEventTime - t) <= convergence))) {
392 return false;
393 }
394
395 if (nextAction == EventHandler.Action.RESET_STATE) {
396 handler.resetState(t, y);
397 }
398 pendingEvent = false;
399 pendingEventTime = Double.NaN;
400
401 return (nextAction == EventHandler.Action.RESET_STATE) ||
402 (nextAction == EventHandler.Action.RESET_DERIVATIVES);
403
404 }
405
406 /** Local wrapper to propagate exceptions. */
407 private static class LocalMaxCountExceededException extends RuntimeException {
408
409 /** Serializable UID. */
410 private static final long serialVersionUID = 20120901L;
411
412 /** Wrapped exception. */
413 private final MaxCountExceededException wrapped;
414
415 /** Simple constructor.
416 * @param exception exception to wrap
417 */
418 LocalMaxCountExceededException(final MaxCountExceededException exception) {
419 wrapped = exception;
420 }
421
422 /** Get the wrapped exception.
423 * @return wrapped exception
424 */
425 public MaxCountExceededException getException() {
426 return wrapped;
427 }
428
429 }
430
431}
Note: See TracBrowser for help on using the repository browser.