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