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 * https://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
18 /*
19 * This is not the original file distributed by the Apache Software Foundation
20 * It has been modified by the Hipparchus project
21 */
22
23 package org.hipparchus.ode;
24
25 import java.util.ArrayList;
26 import java.util.Collections;
27 import java.util.Comparator;
28 import java.util.List;
29 import java.util.PriorityQueue;
30 import java.util.Queue;
31 import java.util.stream.Collectors;
32 import java.util.stream.Stream;
33
34 import org.hipparchus.exception.MathIllegalArgumentException;
35 import org.hipparchus.exception.MathIllegalStateException;
36 import org.hipparchus.ode.events.Action;
37 import org.hipparchus.ode.events.DetectorBasedEventState;
38 import org.hipparchus.ode.events.EventOccurrence;
39 import org.hipparchus.ode.events.EventState;
40 import org.hipparchus.ode.events.ODEEventDetector;
41 import org.hipparchus.ode.events.ODEStepEndHandler;
42 import org.hipparchus.ode.events.StepEndEventState;
43 import org.hipparchus.ode.sampling.AbstractODEStateInterpolator;
44 import org.hipparchus.ode.sampling.ODEStateInterpolator;
45 import org.hipparchus.ode.sampling.ODEStepHandler;
46 import org.hipparchus.util.FastMath;
47 import org.hipparchus.util.Incrementor;
48
49 /**
50 * Base class managing common boilerplate for all integrators.
51 */
52 public abstract class AbstractIntegrator implements ODEIntegrator {
53
54 /** Step handler. */
55 private List<ODEStepHandler> stepHandlers;
56
57 /** Current step start time. */
58 private ODEStateAndDerivative stepStart;
59
60 /** Current stepsize. */
61 private double stepSize;
62
63 /** Indicator for last step. */
64 private boolean isLastStep;
65
66 /** Indicator that a state or derivative reset was triggered by some event. */
67 private boolean resetOccurred;
68
69 /** Events states related to event detectors. */
70 private List<DetectorBasedEventState> detectorBasedEventsStates;
71
72 /** Events states related to step end. */
73 private List<StepEndEventState> stepEndEventsStates;
74
75 /** Initialization indicator of events states. */
76 private boolean statesInitialized;
77
78 /** Name of the method. */
79 private final String name;
80
81 /** Counter for number of evaluations. */
82 private Incrementor evaluations;
83
84 /** Differential equations to integrate. */
85 private ExpandableODE equations;
86
87 /** Build an instance.
88 * @param name name of the method
89 */
90 protected AbstractIntegrator(final String name) {
91 this.name = name;
92 stepHandlers = new ArrayList<>();
93 stepStart = null;
94 stepSize = Double.NaN;
95 detectorBasedEventsStates = new ArrayList<>();
96 stepEndEventsStates = new ArrayList<>();
97 statesInitialized = false;
98 evaluations = new Incrementor();
99 }
100
101 /** {@inheritDoc} */
102 @Override
103 public String getName() {
104 return name;
105 }
106
107 /** {@inheritDoc} */
108 @Override
109 public void addStepHandler(final ODEStepHandler handler) {
110 stepHandlers.add(handler);
111 }
112
113 /** {@inheritDoc} */
114 @Override
115 public List<ODEStepHandler> getStepHandlers() {
116 return Collections.unmodifiableList(stepHandlers);
117 }
118
119 /** {@inheritDoc} */
120 @Override
121 public void clearStepHandlers() {
122 stepHandlers.clear();
123 }
124
125 /** {@inheritDoc} */
126 @Override
127 public void addEventDetector(final ODEEventDetector detector) {
128 detectorBasedEventsStates.add(new DetectorBasedEventState(detector));
129 }
130
131 /** {@inheritDoc} */
132 @Override
133 public List<ODEEventDetector> getEventDetectors() {
134 return detectorBasedEventsStates.stream().map(DetectorBasedEventState::getEventDetector).collect(Collectors.toList());
135 }
136
137 /** {@inheritDoc} */
138 @Override
139 public void clearEventDetectors() {
140 detectorBasedEventsStates.clear();
141 }
142
143 /** {@inheritDoc} */
144 @Override
145 public void addStepEndHandler(ODEStepEndHandler handler) {
146 stepEndEventsStates.add(new StepEndEventState(handler));
147 }
148
149 /** {@inheritDoc} */
150 @Override
151 public List<ODEStepEndHandler> getStepEndHandlers() {
152 return stepEndEventsStates.stream().map(StepEndEventState::getHandler).collect(Collectors.toList());
153 }
154
155 /** {@inheritDoc} */
156 @Override
157 public void clearStepEndHandlers() {
158 stepEndEventsStates.clear();
159 }
160
161 /** {@inheritDoc} */
162 @Override
163 public double getCurrentSignedStepsize() {
164 return stepSize;
165 }
166
167 /** {@inheritDoc} */
168 @Override
169 public void setMaxEvaluations(int maxEvaluations) {
170 evaluations = evaluations.withMaximalCount((maxEvaluations < 0) ? Integer.MAX_VALUE : maxEvaluations);
171 }
172
173 /** {@inheritDoc} */
174 @Override
175 public int getMaxEvaluations() {
176 return evaluations.getMaximalCount();
177 }
178
179 /** {@inheritDoc} */
180 @Override
181 public int getEvaluations() {
182 return evaluations.getCount();
183 }
184
185 /**
186 * Prepare the start of an integration.
187 *
188 * @param eqn equations to integrate
189 * @param s0 initial state vector
190 * @param t target time for the integration
191 * @return Initial state with computed derivatives.
192 */
193 protected ODEStateAndDerivative initIntegration(final ExpandableODE eqn,
194 final ODEState s0, final double t) {
195
196 this.equations = eqn;
197 evaluations = evaluations.withCount(0);
198
199 // initialize ODE
200 eqn.init(s0, t);
201
202 // set up derivatives of initial state (including primary and secondary components)
203 final double t0 = s0.getTime();
204 final double[] y0 = s0.getCompleteState();
205 final double[] y0Dot = computeDerivatives(t0, y0);
206
207 // built the state
208 final ODEStateAndDerivative s0WithDerivatives =
209 eqn.getMapper().mapStateAndDerivative(t0, y0, y0Dot);
210
211 // initialize detector based event states (both and step end based)
212 detectorBasedEventsStates.forEach(s -> {
213 s.init(s0WithDerivatives, t);
214 s.getEventDetector().getHandler().init(s0WithDerivatives, t, s.getEventDetector());
215 });
216
217 // initialize step end based event states
218 stepEndEventsStates.forEach(s -> {
219 s.init(s0WithDerivatives, t);
220 s.getHandler().init(s0WithDerivatives, t);
221 });
222
223 // initialize step handlers
224 for (ODEStepHandler handler : stepHandlers) {
225 handler.init(s0WithDerivatives, t);
226 }
227
228 setStateInitialized(false);
229
230 return s0WithDerivatives;
231
232 }
233
234 /** Get the differential equations to integrate.
235 * @return differential equations to integrate
236 */
237 protected ExpandableODE getEquations() {
238 return equations;
239 }
240
241 /** Get the evaluations counter.
242 * @return evaluations counter
243 */
244 protected Incrementor getEvaluationsCounter() {
245 return evaluations;
246 }
247
248 /** Compute the derivatives and check the number of evaluations.
249 * @param t current value of the independent <I>time</I> variable
250 * @param y array containing the current value of the state vector
251 * @return state completed with derivatives
252 * @exception MathIllegalArgumentException if arrays dimensions do not match equations settings
253 * @exception MathIllegalStateException if the number of functions evaluations is exceeded
254 * @exception NullPointerException if the ODE equations have not been set (i.e. if this method
255 * is called outside of a call to {@link #integrate(ExpandableODE, ODEState, double) integrate}
256 */
257 public double[] computeDerivatives(final double t, final double[] y)
258 throws MathIllegalArgumentException, MathIllegalStateException, NullPointerException {
259 evaluations.increment();
260 return equations.computeDerivatives(t, y);
261 }
262
263 /** Increment evaluations of derivatives.
264 *
265 * @param nTimes number of evaluations to increment
266 */
267 protected void incrementEvaluations(final int nTimes) {
268 evaluations.increment(nTimes);
269 }
270
271 /** Set the stateInitialized flag.
272 * <p>This method must be called by integrators with the value
273 * {@code false} before they start integration, so a proper lazy
274 * initialization is done automatically on the first step.</p>
275 * @param stateInitialized new value for the flag
276 */
277 protected void setStateInitialized(final boolean stateInitialized) {
278 this.statesInitialized = stateInitialized;
279 }
280
281 /** Accept a step, triggering events and step handlers.
282 * @param interpolator step interpolator
283 * @param tEnd final integration time
284 * @return state at end of step
285 * @exception MathIllegalStateException if the interpolator throws one because
286 * the number of functions evaluations is exceeded
287 * @exception MathIllegalArgumentException if the location of an event cannot be bracketed
288 * @exception MathIllegalArgumentException if arrays dimensions do not match equations settings
289 */
290 protected ODEStateAndDerivative acceptStep(final AbstractODEStateInterpolator interpolator,
291 final double tEnd)
292 throws MathIllegalArgumentException, MathIllegalStateException {
293
294 ODEStateAndDerivative previousState = interpolator.getGlobalPreviousState();
295 final ODEStateAndDerivative currentState = interpolator.getGlobalCurrentState();
296 ODEStateInterpolator restricted = interpolator;
297
298
299 // initialize the events states if needed
300 if (!statesInitialized) {
301 // initialize event states
302 detectorBasedEventsStates.forEach(s -> s.reinitializeBegin(interpolator));
303 statesInitialized = true;
304 }
305
306 // set end of step
307 stepEndEventsStates.forEach(s -> s.setStepEnd(currentState.getTime()));
308
309 for (final ODEStepHandler handler : stepHandlers) {
310 handler.updateOnStep(interpolator);
311 }
312
313 // search for next events that may occur during the step
314 final int orderingSign = interpolator.isForward() ? +1 : -1;
315 final Queue<EventState> occurringEvents = new PriorityQueue<>(new Comparator<EventState>() {
316 /** {@inheritDoc} */
317 @Override
318 public int compare(final EventState es0, final EventState es1) {
319 return orderingSign * Double.compare(es0.getEventTime(), es1.getEventTime());
320 }
321 });
322
323 resetOccurred = false;
324 boolean doneWithStep = false;
325 resetEvents:
326 do {
327
328 // Evaluate all event detectors and end steps for events
329 occurringEvents.clear();
330 final ODEStateInterpolator finalRestricted = restricted;
331 Stream.concat(detectorBasedEventsStates.stream(), stepEndEventsStates.stream()).
332 forEach(s -> { if (s.evaluateStep(finalRestricted)) {
333 // the event occurs during the current step
334 occurringEvents.add(s);
335 }
336 });
337
338 do {
339
340 eventLoop:
341 while (!occurringEvents.isEmpty()) {
342
343 // handle the chronologically first event
344 final EventState currentEvent = occurringEvents.poll();
345
346 // get state at event time
347 ODEStateAndDerivative eventState = restricted.getInterpolatedState(currentEvent.getEventTime());
348
349 // restrict the interpolator to the first part of the step, up to the event
350 restricted = restricted.restrictStep(previousState, eventState);
351
352 // try to advance all event states related to detectors to current time
353 for (final DetectorBasedEventState state : detectorBasedEventsStates) {
354 if (state != currentEvent && state.tryAdvance(eventState, interpolator)) {
355 // we need to handle another event first
356 // remove event we just updated to prevent heap corruption
357 occurringEvents.remove(state);
358 // add it back to update its position in the heap
359 occurringEvents.add(state);
360 // re-queue the event we were processing
361 occurringEvents.add(currentEvent);
362 continue eventLoop;
363 }
364 }
365 // all event detectors agree we can advance to the current event time
366
367 // handle the first part of the step, up to the event
368 for (final ODEStepHandler handler : stepHandlers) {
369 handler.handleStep(restricted);
370 }
371
372 // acknowledge event occurrence
373 final EventOccurrence occurrence = currentEvent.doEvent(eventState);
374 final Action action = occurrence.getAction();
375 isLastStep = action == Action.STOP;
376
377 if (isLastStep) {
378
379 // ensure the event is after the root if it is returned STOP
380 // this lets the user integrate to a STOP event and then restart
381 // integration from the same time.
382 final ODEStateAndDerivative savedState = eventState;
383 eventState = interpolator.getInterpolatedState(occurrence.getStopTime());
384 restricted = interpolator.restrictStep(savedState, eventState);
385
386 // handle the almost zero size last part of the final step, at event time
387 for (final ODEStepHandler handler : stepHandlers) {
388 handler.handleStep(restricted);
389 handler.finish(restricted.getCurrentState());
390 }
391
392 }
393
394 if (isLastStep) {
395 // the event asked to stop integration
396 return eventState;
397 }
398
399 if (action == Action.RESET_DERIVATIVES || action == Action.RESET_STATE) {
400 // some event handler has triggered changes that
401 // invalidate the derivatives, we need to recompute them
402 final ODEState newState = occurrence.getNewState();
403 final double[] y = newState.getCompleteState();
404 final double[] yDot = computeDerivatives(newState.getTime(), y);
405 resetOccurred = true;
406 final ODEStateAndDerivative newStateAndDerivative = equations.getMapper().mapStateAndDerivative(newState.getTime(),
407 y, yDot);
408 detectorBasedEventsStates.forEach(e -> e.getEventDetector().reset(newStateAndDerivative, tEnd));
409 return newStateAndDerivative;
410 }
411 // at this point action == Action.CONTINUE or Action.RESET_EVENTS
412
413 // prepare handling of the remaining part of the step
414 previousState = eventState;
415 restricted = restricted.restrictStep(eventState, currentState);
416
417 if (action == Action.RESET_EVENTS) {
418 continue resetEvents;
419 }
420
421 // at this point action == Action.CONTINUE
422 // check if the same event occurs again in the remaining part of the step
423 if (currentEvent.evaluateStep(restricted)) {
424 // the event occurs during the current step
425 occurringEvents.add(currentEvent);
426 }
427
428 }
429
430 // last part of the step, after the last event. Advance all event
431 // detectors to the end of the step. Should never find new events unless
432 // a previous event modified the g function of another event detector when
433 // it returned Action.CONTINUE. Detecting such events here is unreliable
434 // and RESET_EVENTS should be used instead. Other option is to replace
435 // tryAdvance(...) with a doAdvance(...) that throws an exception when
436 // the g function sign is not as expected.
437 for (final DetectorBasedEventState state : detectorBasedEventsStates) {
438 if (state.tryAdvance(currentState, interpolator)) {
439 occurringEvents.add(state);
440 }
441 }
442
443 } while (!occurringEvents.isEmpty());
444
445 doneWithStep = true;
446 } while (!doneWithStep);
447
448 isLastStep = isLastStep || FastMath.abs(currentState.getTime() - tEnd) < FastMath.ulp(tEnd);
449
450 // handle the remaining part of the step, after all events if any
451 for (ODEStepHandler handler : stepHandlers) {
452 handler.handleStep(restricted);
453 if (isLastStep) {
454 handler.finish(restricted.getCurrentState());
455 }
456 }
457
458 return currentState;
459
460 }
461
462 /** Check the integration span.
463 * @param initialState initial state
464 * @param t target time for the integration
465 * @exception MathIllegalArgumentException if integration span is too small
466 * @exception MathIllegalArgumentException if adaptive step size integrators
467 * tolerance arrays dimensions are not compatible with equations settings
468 */
469 protected void sanityChecks(final ODEState initialState, final double t)
470 throws MathIllegalArgumentException {
471
472 final double threshold = 1000 * FastMath.ulp(FastMath.max(FastMath.abs(initialState.getTime()),
473 FastMath.abs(t)));
474 final double dt = FastMath.abs(initialState.getTime() - t);
475 if (dt < threshold) {
476 throw new MathIllegalArgumentException(LocalizedODEFormats.TOO_SMALL_INTEGRATION_INTERVAL,
477 dt, threshold, false);
478 }
479
480 }
481
482 /** Check if a reset occurred while last step was accepted.
483 * @return true if a reset occurred while last step was accepted
484 */
485 protected boolean resetOccurred() {
486 return resetOccurred;
487 }
488
489 /** Set the current step size.
490 * @param stepSize step size to set
491 */
492 protected void setStepSize(final double stepSize) {
493 this.stepSize = stepSize;
494 }
495
496 /** Get the current step size.
497 * @return current step size
498 */
499 protected double getStepSize() {
500 return stepSize;
501 }
502 /** Set current step start.
503 * @param stepStart step start
504 */
505 protected void setStepStart(final ODEStateAndDerivative stepStart) {
506 this.stepStart = stepStart;
507 }
508
509 /** {@inheritDoc} */
510 @Override
511 public ODEStateAndDerivative getStepStart() {
512 return stepStart;
513 }
514
515 /** Set the last state flag.
516 * @param isLastStep if true, this step is the last one
517 */
518 protected void setIsLastStep(final boolean isLastStep) {
519 this.isLastStep = isLastStep;
520 }
521
522 /** Check if this step is the last one.
523 * @return true if this step is the last one
524 */
525 protected boolean isLastStep() {
526 return isLastStep;
527 }
528
529 }