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