View Javadoc
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 }