source: src/main/java/agents/anac/y2019/harddealer/math3/ode/AbstractFieldIntegrator.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;
19
20import java.util.ArrayList;
21import java.util.Collection;
22import java.util.Collections;
23import java.util.Comparator;
24import java.util.Iterator;
25import java.util.List;
26import java.util.SortedSet;
27import java.util.TreeSet;
28
29import agents.anac.y2019.harddealer.math3.Field;
30import agents.anac.y2019.harddealer.math3.RealFieldElement;
31import agents.anac.y2019.harddealer.math3.analysis.solvers.BracketedRealFieldUnivariateSolver;
32import agents.anac.y2019.harddealer.math3.analysis.solvers.FieldBracketingNthOrderBrentSolver;
33import agents.anac.y2019.harddealer.math3.exception.DimensionMismatchException;
34import agents.anac.y2019.harddealer.math3.exception.MaxCountExceededException;
35import agents.anac.y2019.harddealer.math3.exception.NoBracketingException;
36import agents.anac.y2019.harddealer.math3.exception.NumberIsTooSmallException;
37import agents.anac.y2019.harddealer.math3.exception.util.LocalizedFormats;
38import agents.anac.y2019.harddealer.math3.ode.events.FieldEventHandler;
39import agents.anac.y2019.harddealer.math3.ode.events.FieldEventState;
40import agents.anac.y2019.harddealer.math3.ode.sampling.AbstractFieldStepInterpolator;
41import agents.anac.y2019.harddealer.math3.ode.sampling.FieldStepHandler;
42import agents.anac.y2019.harddealer.math3.util.FastMath;
43import agents.anac.y2019.harddealer.math3.util.IntegerSequence;
44
45/**
46 * Base class managing common boilerplate for all integrators.
47 * @param <T> the type of the field elements
48 * @since 3.6
49 */
50public abstract class AbstractFieldIntegrator<T extends RealFieldElement<T>> implements FirstOrderFieldIntegrator<T> {
51
52 /** Default relative accuracy. */
53 private static final double DEFAULT_RELATIVE_ACCURACY = 1e-14;
54
55 /** Default function value accuracy. */
56 private static final double DEFAULT_FUNCTION_VALUE_ACCURACY = 1e-15;
57
58 /** Step handler. */
59 private Collection<FieldStepHandler<T>> stepHandlers;
60
61 /** Current step start. */
62 private FieldODEStateAndDerivative<T> stepStart;
63
64 /** Current stepsize. */
65 private T stepSize;
66
67 /** Indicator for last step. */
68 private boolean isLastStep;
69
70 /** Indicator that a state or derivative reset was triggered by some event. */
71 private boolean resetOccurred;
72
73 /** Field to which the time and state vector elements belong. */
74 private final Field<T> field;
75
76 /** Events states. */
77 private Collection<FieldEventState<T>> eventsStates;
78
79 /** Initialization indicator of events states. */
80 private boolean statesInitialized;
81
82 /** Name of the method. */
83 private final String name;
84
85 /** Counter for number of evaluations. */
86 private IntegerSequence.Incrementor evaluations;
87
88 /** Differential equations to integrate. */
89 private transient FieldExpandableODE<T> equations;
90
91 /** Build an instance.
92 * @param field field to which the time and state vector elements belong
93 * @param name name of the method
94 */
95 protected AbstractFieldIntegrator(final Field<T> field, final String name) {
96 this.field = field;
97 this.name = name;
98 stepHandlers = new ArrayList<FieldStepHandler<T>>();
99 stepStart = null;
100 stepSize = null;
101 eventsStates = new ArrayList<FieldEventState<T>>();
102 statesInitialized = false;
103 evaluations = IntegerSequence.Incrementor.create().withMaximalCount(Integer.MAX_VALUE);
104 }
105
106 /** Get the field to which state vector elements belong.
107 * @return field to which state vector elements belong
108 */
109 public Field<T> getField() {
110 return field;
111 }
112
113 /** {@inheritDoc} */
114 public String getName() {
115 return name;
116 }
117
118 /** {@inheritDoc} */
119 public void addStepHandler(final FieldStepHandler<T> handler) {
120 stepHandlers.add(handler);
121 }
122
123 /** {@inheritDoc} */
124 public Collection<FieldStepHandler<T>> getStepHandlers() {
125 return Collections.unmodifiableCollection(stepHandlers);
126 }
127
128 /** {@inheritDoc} */
129 public void clearStepHandlers() {
130 stepHandlers.clear();
131 }
132
133 /** {@inheritDoc} */
134 public void addEventHandler(final FieldEventHandler<T> handler,
135 final double maxCheckInterval,
136 final double convergence,
137 final int maxIterationCount) {
138 addEventHandler(handler, maxCheckInterval, convergence,
139 maxIterationCount,
140 new FieldBracketingNthOrderBrentSolver<T>(field.getZero().add(DEFAULT_RELATIVE_ACCURACY),
141 field.getZero().add(convergence),
142 field.getZero().add(DEFAULT_FUNCTION_VALUE_ACCURACY),
143 5));
144 }
145
146 /** {@inheritDoc} */
147 public void addEventHandler(final FieldEventHandler<T> handler,
148 final double maxCheckInterval,
149 final double convergence,
150 final int maxIterationCount,
151 final BracketedRealFieldUnivariateSolver<T> solver) {
152 eventsStates.add(new FieldEventState<T>(handler, maxCheckInterval, field.getZero().add(convergence),
153 maxIterationCount, solver));
154 }
155
156 /** {@inheritDoc} */
157 public Collection<FieldEventHandler<T>> getEventHandlers() {
158 final List<FieldEventHandler<T>> list = new ArrayList<FieldEventHandler<T>>(eventsStates.size());
159 for (FieldEventState<T> state : eventsStates) {
160 list.add(state.getEventHandler());
161 }
162 return Collections.unmodifiableCollection(list);
163 }
164
165 /** {@inheritDoc} */
166 public void clearEventHandlers() {
167 eventsStates.clear();
168 }
169
170 /** {@inheritDoc} */
171 public FieldODEStateAndDerivative<T> getCurrentStepStart() {
172 return stepStart;
173 }
174
175 /** {@inheritDoc} */
176 public T getCurrentSignedStepsize() {
177 return stepSize;
178 }
179
180 /** {@inheritDoc} */
181 public void setMaxEvaluations(int maxEvaluations) {
182 evaluations = evaluations.withMaximalCount((maxEvaluations < 0) ? Integer.MAX_VALUE : maxEvaluations);
183 }
184
185 /** {@inheritDoc} */
186 public int getMaxEvaluations() {
187 return evaluations.getMaximalCount();
188 }
189
190 /** {@inheritDoc} */
191 public int getEvaluations() {
192 return evaluations.getCount();
193 }
194
195 /** Prepare the start of an integration.
196 * @param eqn equations to integrate
197 * @param t0 start value of the independent <i>time</i> variable
198 * @param y0 array containing the start value of the state vector
199 * @param t target time for the integration
200 * @return initial state with derivatives added
201 */
202 protected FieldODEStateAndDerivative<T> initIntegration(final FieldExpandableODE<T> eqn,
203 final T t0, final T[] y0, final T t) {
204
205 this.equations = eqn;
206 evaluations = evaluations.withStart(0);
207
208 // initialize ODE
209 eqn.init(t0, y0, t);
210
211 // set up derivatives of initial state
212 final T[] y0Dot = computeDerivatives(t0, y0);
213 final FieldODEStateAndDerivative<T> state0 = new FieldODEStateAndDerivative<T>(t0, y0, y0Dot);
214
215 // initialize event handlers
216 for (final FieldEventState<T> state : eventsStates) {
217 state.getEventHandler().init(state0, t);
218 }
219
220 // initialize step handlers
221 for (FieldStepHandler<T> handler : stepHandlers) {
222 handler.init(state0, t);
223 }
224
225 setStateInitialized(false);
226
227 return state0;
228
229 }
230
231 /** Get the differential equations to integrate.
232 * @return differential equations to integrate
233 */
234 protected FieldExpandableODE<T> getEquations() {
235 return equations;
236 }
237
238 /** Get the evaluations counter.
239 * @return evaluations counter
240 */
241 protected IntegerSequence.Incrementor getEvaluationsCounter() {
242 return evaluations;
243 }
244
245 /** Compute the derivatives and check the number of evaluations.
246 * @param t current value of the independent <I>time</I> variable
247 * @param y array containing the current value of the state vector
248 * @return state completed with derivatives
249 * @exception DimensionMismatchException if arrays dimensions do not match equations settings
250 * @exception MaxCountExceededException if the number of functions evaluations is exceeded
251 * @exception NullPointerException if the ODE equations have not been set (i.e. if this method
252 * is called outside of a call to {@link #integrate(FieldExpandableODE, FieldODEState,
253 * RealFieldElement) integrate}
254 */
255 public T[] computeDerivatives(final T t, final T[] y)
256 throws DimensionMismatchException, MaxCountExceededException, NullPointerException {
257 evaluations.increment();
258 return equations.computeDerivatives(t, y);
259 }
260
261 /** Set the stateInitialized flag.
262 * <p>This method must be called by integrators with the value
263 * {@code false} before they start integration, so a proper lazy
264 * initialization is done automatically on the first step.</p>
265 * @param stateInitialized new value for the flag
266 */
267 protected void setStateInitialized(final boolean stateInitialized) {
268 this.statesInitialized = stateInitialized;
269 }
270
271 /** Accept a step, triggering events and step handlers.
272 * @param interpolator step interpolator
273 * @param tEnd final integration time
274 * @return state at end of step
275 * @exception MaxCountExceededException if the interpolator throws one because
276 * the number of functions evaluations is exceeded
277 * @exception NoBracketingException if the location of an event cannot be bracketed
278 * @exception DimensionMismatchException if arrays dimensions do not match equations settings
279 */
280 protected FieldODEStateAndDerivative<T> acceptStep(final AbstractFieldStepInterpolator<T> interpolator,
281 final T tEnd)
282 throws MaxCountExceededException, DimensionMismatchException, NoBracketingException {
283
284 FieldODEStateAndDerivative<T> previousState = interpolator.getGlobalPreviousState();
285 final FieldODEStateAndDerivative<T> currentState = interpolator.getGlobalCurrentState();
286
287 // initialize the events states if needed
288 if (! statesInitialized) {
289 for (FieldEventState<T> state : eventsStates) {
290 state.reinitializeBegin(interpolator);
291 }
292 statesInitialized = true;
293 }
294
295 // search for next events that may occur during the step
296 final int orderingSign = interpolator.isForward() ? +1 : -1;
297 SortedSet<FieldEventState<T>> occurringEvents = new TreeSet<FieldEventState<T>>(new Comparator<FieldEventState<T>>() {
298
299 /** {@inheritDoc} */
300 public int compare(FieldEventState<T> es0, FieldEventState<T> es1) {
301 return orderingSign * Double.compare(es0.getEventTime().getReal(), es1.getEventTime().getReal());
302 }
303
304 });
305
306 for (final FieldEventState<T> state : eventsStates) {
307 if (state.evaluateStep(interpolator)) {
308 // the event occurs during the current step
309 occurringEvents.add(state);
310 }
311 }
312
313 AbstractFieldStepInterpolator<T> restricted = interpolator;
314 while (!occurringEvents.isEmpty()) {
315
316 // handle the chronologically first event
317 final Iterator<FieldEventState<T>> iterator = occurringEvents.iterator();
318 final FieldEventState<T> currentEvent = iterator.next();
319 iterator.remove();
320
321 // get state at event time
322 final FieldODEStateAndDerivative<T> eventState = restricted.getInterpolatedState(currentEvent.getEventTime());
323
324 // restrict the interpolator to the first part of the step, up to the event
325 restricted = restricted.restrictStep(previousState, eventState);
326
327 // advance all event states to current time
328 for (final FieldEventState<T> state : eventsStates) {
329 state.stepAccepted(eventState);
330 isLastStep = isLastStep || state.stop();
331 }
332
333 // handle the first part of the step, up to the event
334 for (final FieldStepHandler<T> handler : stepHandlers) {
335 handler.handleStep(restricted, isLastStep);
336 }
337
338 if (isLastStep) {
339 // the event asked to stop integration
340 return eventState;
341 }
342
343 FieldODEState<T> newState = null;
344 resetOccurred = false;
345 for (final FieldEventState<T> state : eventsStates) {
346 newState = state.reset(eventState);
347 if (newState != null) {
348 // some event handler has triggered changes that
349 // invalidate the derivatives, we need to recompute them
350 final T[] y = equations.getMapper().mapState(newState);
351 final T[] yDot = computeDerivatives(newState.getTime(), y);
352 resetOccurred = true;
353 return equations.getMapper().mapStateAndDerivative(newState.getTime(), y, yDot);
354 }
355 }
356
357 // prepare handling of the remaining part of the step
358 previousState = eventState;
359 restricted = restricted.restrictStep(eventState, currentState);
360
361 // check if the same event occurs again in the remaining part of the step
362 if (currentEvent.evaluateStep(restricted)) {
363 // the event occurs during the current step
364 occurringEvents.add(currentEvent);
365 }
366
367 }
368
369 // last part of the step, after the last event
370 for (final FieldEventState<T> state : eventsStates) {
371 state.stepAccepted(currentState);
372 isLastStep = isLastStep || state.stop();
373 }
374 isLastStep = isLastStep || currentState.getTime().subtract(tEnd).abs().getReal() <= FastMath.ulp(tEnd.getReal());
375
376 // handle the remaining part of the step, after all events if any
377 for (FieldStepHandler<T> handler : stepHandlers) {
378 handler.handleStep(restricted, isLastStep);
379 }
380
381 return currentState;
382
383 }
384
385 /** Check the integration span.
386 * @param eqn set of differential equations
387 * @param t target time for the integration
388 * @exception NumberIsTooSmallException if integration span is too small
389 * @exception DimensionMismatchException if adaptive step size integrators
390 * tolerance arrays dimensions are not compatible with equations settings
391 */
392 protected void sanityChecks(final FieldODEState<T> eqn, final T t)
393 throws NumberIsTooSmallException, DimensionMismatchException {
394
395 final double threshold = 1000 * FastMath.ulp(FastMath.max(FastMath.abs(eqn.getTime().getReal()),
396 FastMath.abs(t.getReal())));
397 final double dt = eqn.getTime().subtract(t).abs().getReal();
398 if (dt <= threshold) {
399 throw new NumberIsTooSmallException(LocalizedFormats.TOO_SMALL_INTEGRATION_INTERVAL,
400 dt, threshold, false);
401 }
402
403 }
404
405 /** Check if a reset occurred while last step was accepted.
406 * @return true if a reset occurred while last step was accepted
407 */
408 protected boolean resetOccurred() {
409 return resetOccurred;
410 }
411
412 /** Set the current step size.
413 * @param stepSize step size to set
414 */
415 protected void setStepSize(final T stepSize) {
416 this.stepSize = stepSize;
417 }
418
419 /** Get the current step size.
420 * @return current step size
421 */
422 protected T getStepSize() {
423 return stepSize;
424 }
425 /** Set current step start.
426 * @param stepStart step start
427 */
428 protected void setStepStart(final FieldODEStateAndDerivative<T> stepStart) {
429 this.stepStart = stepStart;
430 }
431
432 /** Getcurrent step start.
433 * @return current step start
434 */
435 protected FieldODEStateAndDerivative<T> getStepStart() {
436 return stepStart;
437 }
438
439 /** Set the last state flag.
440 * @param isLastStep if true, this step is the last one
441 */
442 protected void setIsLastStep(final boolean isLastStep) {
443 this.isLastStep = isLastStep;
444 }
445
446 /** Check if this step is the last one.
447 * @return true if this step is the last one
448 */
449 protected boolean isLastStep() {
450 return isLastStep;
451 }
452
453}
Note: See TracBrowser for help on using the repository browser.