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 // search for next events that may occur during the step
310 final int orderingSign = interpolator.isForward() ? +1 : -1;
311 final Queue<EventState> occurringEvents = new PriorityQueue<>(new Comparator<EventState>() {
312 /** {@inheritDoc} */
313 @Override
314 public int compare(final EventState es0, final EventState es1) {
315 return orderingSign * Double.compare(es0.getEventTime(), es1.getEventTime());
316 }
317 });
318
319 resetOccurred = false;
320 boolean doneWithStep = false;
321 resetEvents:
322 do {
323
324 // Evaluate all event detectors and end steps for events
325 occurringEvents.clear();
326 final ODEStateInterpolator finalRestricted = restricted;
327 Stream.concat(detectorBasedEventsStates.stream(), stepEndEventsStates.stream()).
328 forEach(s -> { if (s.evaluateStep(finalRestricted)) {
329 // the event occurs during the current step
330 occurringEvents.add(s);
331 }
332 });
333
334 do {
335
336 eventLoop:
337 while (!occurringEvents.isEmpty()) {
338
339 // handle the chronologically first event
340 final EventState currentEvent = occurringEvents.poll();
341
342 // get state at event time
343 ODEStateAndDerivative eventState = restricted.getInterpolatedState(currentEvent.getEventTime());
344
345 // restrict the interpolator to the first part of the step, up to the event
346 restricted = restricted.restrictStep(previousState, eventState);
347
348 // try to advance all event states related to detectors to current time
349 for (final DetectorBasedEventState state : detectorBasedEventsStates) {
350 if (state != currentEvent && state.tryAdvance(eventState, interpolator)) {
351 // we need to handle another event first
352 // remove event we just updated to prevent heap corruption
353 occurringEvents.remove(state);
354 // add it back to update its position in the heap
355 occurringEvents.add(state);
356 // re-queue the event we were processing
357 occurringEvents.add(currentEvent);
358 continue eventLoop;
359 }
360 }
361 // all event detectors agree we can advance to the current event time
362
363 // handle the first part of the step, up to the event
364 for (final ODEStepHandler handler : stepHandlers) {
365 handler.handleStep(restricted);
366 }
367
368 // acknowledge event occurrence
369 final EventOccurrence occurrence = currentEvent.doEvent(eventState);
370 final Action action = occurrence.getAction();
371 isLastStep = action == Action.STOP;
372
373 if (isLastStep) {
374
375 // ensure the event is after the root if it is returned STOP
376 // this lets the user integrate to a STOP event and then restart
377 // integration from the same time.
378 final ODEStateAndDerivative savedState = eventState;
379 eventState = interpolator.getInterpolatedState(occurrence.getStopTime());
380 restricted = interpolator.restrictStep(savedState, eventState);
381
382 // handle the almost zero size last part of the final step, at event time
383 for (final ODEStepHandler handler : stepHandlers) {
384 handler.handleStep(restricted);
385 handler.finish(restricted.getCurrentState());
386 }
387
388 }
389
390 if (isLastStep) {
391 // the event asked to stop integration
392 return eventState;
393 }
394
395 if (action == Action.RESET_DERIVATIVES || action == Action.RESET_STATE) {
396 // some event handler has triggered changes that
397 // invalidate the derivatives, we need to recompute them
398 final ODEState newState = occurrence.getNewState();
399 final double[] y = newState.getCompleteState();
400 final double[] yDot = computeDerivatives(newState.getTime(), y);
401 resetOccurred = true;
402 final ODEStateAndDerivative newStateAndDerivative = equations.getMapper().mapStateAndDerivative(newState.getTime(),
403 y, yDot);
404 detectorBasedEventsStates.forEach(e -> e.getEventDetector().reset(newStateAndDerivative, tEnd));
405 return newStateAndDerivative;
406 }
407 // at this point action == Action.CONTINUE or Action.RESET_EVENTS
408
409 // prepare handling of the remaining part of the step
410 previousState = eventState;
411 restricted = restricted.restrictStep(eventState, currentState);
412
413 if (action == Action.RESET_EVENTS) {
414 continue resetEvents;
415 }
416
417 // at this point action == Action.CONTINUE
418 // check if the same event occurs again in the remaining part of the step
419 if (currentEvent.evaluateStep(restricted)) {
420 // the event occurs during the current step
421 occurringEvents.add(currentEvent);
422 }
423
424 }
425
426 // last part of the step, after the last event. Advance all event
427 // detectors to the end of the step. Should never find new events unless
428 // a previous event modified the g function of another event detector when
429 // it returned Action.CONTINUE. Detecting such events here is unreliable
430 // and RESET_EVENTS should be used instead. Other option is to replace
431 // tryAdvance(...) with a doAdvance(...) that throws an exception when
432 // the g function sign is not as expected.
433 for (final DetectorBasedEventState state : detectorBasedEventsStates) {
434 if (state.tryAdvance(currentState, interpolator)) {
435 occurringEvents.add(state);
436 }
437 }
438
439 } while (!occurringEvents.isEmpty());
440
441 doneWithStep = true;
442 } while (!doneWithStep);
443
444 isLastStep = isLastStep || FastMath.abs(currentState.getTime() - tEnd) < FastMath.ulp(tEnd);
445
446 // handle the remaining part of the step, after all events if any
447 for (ODEStepHandler handler : stepHandlers) {
448 handler.handleStep(restricted);
449 if (isLastStep) {
450 handler.finish(restricted.getCurrentState());
451 }
452 }
453
454 return currentState;
455
456 }
457
458 /** Check the integration span.
459 * @param initialState initial state
460 * @param t target time for the integration
461 * @exception MathIllegalArgumentException if integration span is too small
462 * @exception MathIllegalArgumentException if adaptive step size integrators
463 * tolerance arrays dimensions are not compatible with equations settings
464 */
465 protected void sanityChecks(final ODEState initialState, final double t)
466 throws MathIllegalArgumentException {
467
468 final double threshold = 1000 * FastMath.ulp(FastMath.max(FastMath.abs(initialState.getTime()),
469 FastMath.abs(t)));
470 final double dt = FastMath.abs(initialState.getTime() - t);
471 if (dt < threshold) {
472 throw new MathIllegalArgumentException(LocalizedODEFormats.TOO_SMALL_INTEGRATION_INTERVAL,
473 dt, threshold, false);
474 }
475
476 }
477
478 /** Check if a reset occurred while last step was accepted.
479 * @return true if a reset occurred while last step was accepted
480 */
481 protected boolean resetOccurred() {
482 return resetOccurred;
483 }
484
485 /** Set the current step size.
486 * @param stepSize step size to set
487 */
488 protected void setStepSize(final double stepSize) {
489 this.stepSize = stepSize;
490 }
491
492 /** Get the current step size.
493 * @return current step size
494 */
495 protected double getStepSize() {
496 return stepSize;
497 }
498 /** Set current step start.
499 * @param stepStart step start
500 */
501 protected void setStepStart(final ODEStateAndDerivative stepStart) {
502 this.stepStart = stepStart;
503 }
504
505 /** {@inheritDoc} */
506 @Override
507 public ODEStateAndDerivative getStepStart() {
508 return stepStart;
509 }
510
511 /** Set the last state flag.
512 * @param isLastStep if true, this step is the last one
513 */
514 protected void setIsLastStep(final boolean isLastStep) {
515 this.isLastStep = isLastStep;
516 }
517
518 /** Check if this step is the last one.
519 * @return true if this step is the last one
520 */
521 protected boolean isLastStep() {
522 return isLastStep;
523 }
524
525 }