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    *      http://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  package org.apache.commons.math4.legacy.ode;
19  
20  import java.util.ArrayList;
21  import java.util.Collection;
22  import java.util.Collections;
23  import java.util.Comparator;
24  import java.util.Iterator;
25  import java.util.List;
26  import java.util.SortedSet;
27  import java.util.TreeSet;
28  
29  import org.apache.commons.math4.legacy.core.Field;
30  import org.apache.commons.math4.legacy.core.RealFieldElement;
31  import org.apache.commons.math4.legacy.analysis.solvers.BracketedRealFieldUnivariateSolver;
32  import org.apache.commons.math4.legacy.analysis.solvers.FieldBracketingNthOrderBrentSolver;
33  import org.apache.commons.math4.legacy.exception.DimensionMismatchException;
34  import org.apache.commons.math4.legacy.exception.MaxCountExceededException;
35  import org.apache.commons.math4.legacy.exception.NoBracketingException;
36  import org.apache.commons.math4.legacy.exception.NumberIsTooSmallException;
37  import org.apache.commons.math4.legacy.exception.util.LocalizedFormats;
38  import org.apache.commons.math4.legacy.ode.events.FieldEventHandler;
39  import org.apache.commons.math4.legacy.ode.events.FieldEventState;
40  import org.apache.commons.math4.legacy.ode.sampling.AbstractFieldStepInterpolator;
41  import org.apache.commons.math4.legacy.ode.sampling.FieldStepHandler;
42  import org.apache.commons.math4.core.jdkmath.JdkMath;
43  import org.apache.commons.math4.legacy.core.IntegerSequence;
44  
45  /**
46   * Base class managing common boilerplate for all integrators.
47   * @param <T> the type of the field elements
48   * @since 3.6
49   */
50  public abstract class AbstractFieldIntegrator<T extends RealFieldElement<T>> implements FirstOrderFieldIntegrator<T> {
51  
52      /** Default relative accuracy. */
53      private static final double DEFAULT_RELATIVE_ACCURACY = 1e-14;
54  
55      /** Default function value accuracy. */
56      private static final double DEFAULT_FUNCTION_VALUE_ACCURACY = 1e-15;
57  
58      /** Step handler. */
59      private Collection<FieldStepHandler<T>> stepHandlers;
60  
61      /** Current step start. */
62      private FieldODEStateAndDerivative<T> stepStart;
63  
64      /** Current stepsize. */
65      private T stepSize;
66  
67      /** Indicator for last step. */
68      private boolean isLastStep;
69  
70      /** Indicator that a state or derivative reset was triggered by some event. */
71      private boolean resetOccurred;
72  
73      /** Field to which the time and state vector elements belong. */
74      private final Field<T> field;
75  
76      /** Events states. */
77      private Collection<FieldEventState<T>> eventsStates;
78  
79      /** Initialization indicator of events states. */
80      private boolean statesInitialized;
81  
82      /** Name of the method. */
83      private final String name;
84  
85      /** Counter for number of evaluations. */
86      private IntegerSequence.Incrementor evaluations;
87  
88      /** Differential equations to integrate. */
89      private transient FieldExpandableODE<T> equations;
90  
91      /** Build an instance.
92       * @param field field to which the time and state vector elements belong
93       * @param name name of the method
94       */
95      protected AbstractFieldIntegrator(final Field<T> field, final String name) {
96          this.field        = field;
97          this.name         = name;
98          stepHandlers      = new ArrayList<>();
99          stepStart         = null;
100         stepSize          = null;
101         eventsStates      = new ArrayList<>();
102         statesInitialized = false;
103         evaluations       = IntegerSequence.Incrementor.create().withMaximalCount(Integer.MAX_VALUE);
104     }
105 
106     /** Get the field to which state vector elements belong.
107      * @return field to which state vector elements belong
108      */
109     public Field<T> getField() {
110         return field;
111     }
112 
113     /** {@inheritDoc} */
114     @Override
115     public String getName() {
116         return name;
117     }
118 
119     /** {@inheritDoc} */
120     @Override
121     public void addStepHandler(final FieldStepHandler<T> handler) {
122         stepHandlers.add(handler);
123     }
124 
125     /** {@inheritDoc} */
126     @Override
127     public Collection<FieldStepHandler<T>> getStepHandlers() {
128         return Collections.unmodifiableCollection(stepHandlers);
129     }
130 
131     /** {@inheritDoc} */
132     @Override
133     public void clearStepHandlers() {
134         stepHandlers.clear();
135     }
136 
137     /** {@inheritDoc} */
138     @Override
139     public void addEventHandler(final FieldEventHandler<T> handler,
140                                 final double maxCheckInterval,
141                                 final double convergence,
142                                 final int maxIterationCount) {
143         addEventHandler(handler, maxCheckInterval, convergence,
144                         maxIterationCount,
145                         new FieldBracketingNthOrderBrentSolver<>(field.getZero().add(DEFAULT_RELATIVE_ACCURACY),
146                                                                   field.getZero().add(convergence),
147                                                                   field.getZero().add(DEFAULT_FUNCTION_VALUE_ACCURACY),
148                                                                   5));
149     }
150 
151     /** {@inheritDoc} */
152     @Override
153     public void addEventHandler(final FieldEventHandler<T> handler,
154                                 final double maxCheckInterval,
155                                 final double convergence,
156                                 final int maxIterationCount,
157                                 final BracketedRealFieldUnivariateSolver<T> solver) {
158         eventsStates.add(new FieldEventState<>(handler, maxCheckInterval, field.getZero().add(convergence),
159                                                 maxIterationCount, solver));
160     }
161 
162     /** {@inheritDoc} */
163     @Override
164     public Collection<FieldEventHandler<T>> getEventHandlers() {
165         final List<FieldEventHandler<T>> list = new ArrayList<>(eventsStates.size());
166         for (FieldEventState<T> state : eventsStates) {
167             list.add(state.getEventHandler());
168         }
169         return Collections.unmodifiableCollection(list);
170     }
171 
172     /** {@inheritDoc} */
173     @Override
174     public void clearEventHandlers() {
175         eventsStates.clear();
176     }
177 
178     /** {@inheritDoc} */
179     @Override
180     public FieldODEStateAndDerivative<T> getCurrentStepStart() {
181         return stepStart;
182     }
183 
184     /** {@inheritDoc} */
185     @Override
186     public T getCurrentSignedStepsize() {
187         return stepSize;
188     }
189 
190     /** {@inheritDoc} */
191     @Override
192     public void setMaxEvaluations(int maxEvaluations) {
193         evaluations = evaluations.withMaximalCount((maxEvaluations < 0) ? Integer.MAX_VALUE : maxEvaluations);
194     }
195 
196     /** {@inheritDoc} */
197     @Override
198     public int getMaxEvaluations() {
199         return evaluations.getMaximalCount();
200     }
201 
202     /** {@inheritDoc} */
203     @Override
204     public int getEvaluations() {
205         return evaluations.getCount();
206     }
207 
208     /** Prepare the start of an integration.
209      * @param eqn equations to integrate
210      * @param t0 start value of the independent <i>time</i> variable
211      * @param y0 array containing the start value of the state vector
212      * @param t target time for the integration
213      * @return initial state with derivatives added
214      */
215     protected FieldODEStateAndDerivative<T> initIntegration(final FieldExpandableODE<T> eqn,
216                                                             final T t0, final T[] y0, final T t) {
217 
218         this.equations = eqn;
219         evaluations    = evaluations.withStart(0);
220 
221         // initialize ODE
222         eqn.init(t0, y0, t);
223 
224         // set up derivatives of initial state
225         final T[] y0Dot = computeDerivatives(t0, y0);
226         final FieldODEStateAndDerivative<T> state0 = new FieldODEStateAndDerivative<>(t0, y0, y0Dot);
227 
228         // initialize event handlers
229         for (final FieldEventState<T> state : eventsStates) {
230             state.getEventHandler().init(state0, t);
231         }
232 
233         // initialize step handlers
234         for (FieldStepHandler<T> handler : stepHandlers) {
235             handler.init(state0, t);
236         }
237 
238         setStateInitialized(false);
239 
240         return state0;
241     }
242 
243     /** Get the differential equations to integrate.
244      * @return differential equations to integrate
245      */
246     protected FieldExpandableODE<T> getEquations() {
247         return equations;
248     }
249 
250     /** Get the evaluations counter.
251      * @return evaluations counter
252      */
253     protected IntegerSequence.Incrementor getEvaluationsCounter() {
254         return evaluations;
255     }
256 
257     /** Compute the derivatives and check the number of evaluations.
258      * @param t current value of the independent <I>time</I> variable
259      * @param y array containing the current value of the state vector
260      * @return state completed with derivatives
261      * @exception DimensionMismatchException if arrays dimensions do not match equations settings
262      * @exception MaxCountExceededException if the number of functions evaluations is exceeded
263      * @exception NullPointerException if the ODE equations have not been set (i.e. if this method
264      * is called outside of a call to {@link #integrate(FieldExpandableODE, FieldODEState,
265      * RealFieldElement) integrate}
266      */
267     public T[] computeDerivatives(final T t, final T[] y)
268         throws DimensionMismatchException, MaxCountExceededException, NullPointerException {
269         evaluations.increment();
270         return equations.computeDerivatives(t, y);
271     }
272 
273     /** Set the stateInitialized flag.
274      * <p>This method must be called by integrators with the value
275      * {@code false} before they start integration, so a proper lazy
276      * initialization is done automatically on the first step.</p>
277      * @param stateInitialized new value for the flag
278      */
279     protected void setStateInitialized(final boolean stateInitialized) {
280         this.statesInitialized = stateInitialized;
281     }
282 
283     /** Accept a step, triggering events and step handlers.
284      * @param interpolator step interpolator
285      * @param tEnd final integration time
286      * @return state at end of step
287      * @exception MaxCountExceededException if the interpolator throws one because
288      * the number of functions evaluations is exceeded
289      * @exception NoBracketingException if the location of an event cannot be bracketed
290      * @exception DimensionMismatchException if arrays dimensions do not match equations settings
291      */
292     protected FieldODEStateAndDerivative<T> acceptStep(final AbstractFieldStepInterpolator<T> interpolator,
293                                                        final T tEnd)
294         throws MaxCountExceededException, DimensionMismatchException, NoBracketingException {
295 
296             FieldODEStateAndDerivative<T> previousState = interpolator.getGlobalPreviousState();
297             final FieldODEStateAndDerivative<T> currentState = interpolator.getGlobalCurrentState();
298 
299             // initialize the events states if needed
300             if (! statesInitialized) {
301                 for (FieldEventState<T> state : eventsStates) {
302                     state.reinitializeBegin(interpolator);
303                 }
304                 statesInitialized = true;
305             }
306 
307             // search for next events that may occur during the step
308             final int orderingSign = interpolator.isForward() ? +1 : -1;
309             SortedSet<FieldEventState<T>> occurringEvents = new TreeSet<>(new Comparator<FieldEventState<T>>() {
310 
311                 /** {@inheritDoc} */
312                 @Override
313                 public int compare(FieldEventState<T> es0, FieldEventState<T> es1) {
314                     return orderingSign * Double.compare(es0.getEventTime().getReal(), es1.getEventTime().getReal());
315                 }
316             });
317 
318             for (final FieldEventState<T> state : eventsStates) {
319                 if (state.evaluateStep(interpolator)) {
320                     // the event occurs during the current step
321                     occurringEvents.add(state);
322                 }
323             }
324 
325             AbstractFieldStepInterpolator<T> restricted = interpolator;
326             while (!occurringEvents.isEmpty()) {
327 
328                 // handle the chronologically first event
329                 final Iterator<FieldEventState<T>> iterator = occurringEvents.iterator();
330                 final FieldEventState<T> currentEvent = iterator.next();
331                 iterator.remove();
332 
333                 // get state at event time
334                 final FieldODEStateAndDerivative<T> eventState = restricted.getInterpolatedState(currentEvent.getEventTime());
335 
336                 // restrict the interpolator to the first part of the step, up to the event
337                 restricted = restricted.restrictStep(previousState, eventState);
338 
339                 // advance all event states to current time
340                 for (final FieldEventState<T> state : eventsStates) {
341                     state.stepAccepted(eventState);
342                     isLastStep = isLastStep || state.stop();
343                 }
344 
345                 // handle the first part of the step, up to the event
346                 for (final FieldStepHandler<T> handler : stepHandlers) {
347                     handler.handleStep(restricted, isLastStep);
348                 }
349 
350                 if (isLastStep) {
351                     // the event asked to stop integration
352                     return eventState;
353                 }
354 
355                 FieldODEState<T> newState = null;
356                 resetOccurred = false;
357                 for (final FieldEventState<T> state : eventsStates) {
358                     newState = state.reset(eventState);
359                     if (newState != null) {
360                         // some event handler has triggered changes that
361                         // invalidate the derivatives, we need to recompute them
362                         final T[] y    = equations.getMapper().mapState(newState);
363                         final T[] yDot = computeDerivatives(newState.getTime(), y);
364                         resetOccurred = true;
365                         return equations.getMapper().mapStateAndDerivative(newState.getTime(), y, yDot);
366                     }
367                 }
368 
369                 // prepare handling of the remaining part of the step
370                 previousState = eventState;
371                 restricted = restricted.restrictStep(eventState, currentState);
372 
373                 // check if the same event occurs again in the remaining part of the step
374                 if (currentEvent.evaluateStep(restricted)) {
375                     // the event occurs during the current step
376                     occurringEvents.add(currentEvent);
377                 }
378             }
379 
380             // last part of the step, after the last event
381             for (final FieldEventState<T> state : eventsStates) {
382                 state.stepAccepted(currentState);
383                 isLastStep = isLastStep || state.stop();
384             }
385             isLastStep = isLastStep || currentState.getTime().subtract(tEnd).abs().getReal() <= JdkMath.ulp(tEnd.getReal());
386 
387             // handle the remaining part of the step, after all events if any
388             for (FieldStepHandler<T> handler : stepHandlers) {
389                 handler.handleStep(restricted, isLastStep);
390             }
391 
392             return currentState;
393     }
394 
395     /** Check the integration span.
396      * @param eqn set of differential equations
397      * @param t target time for the integration
398      * @exception NumberIsTooSmallException if integration span is too small
399      * @exception DimensionMismatchException if adaptive step size integrators
400      * tolerance arrays dimensions are not compatible with equations settings
401      */
402     protected void sanityChecks(final FieldODEState<T> eqn, final T t)
403         throws NumberIsTooSmallException, DimensionMismatchException {
404 
405         final double threshold = 1000 * JdkMath.ulp(JdkMath.max(JdkMath.abs(eqn.getTime().getReal()),
406                                                                   JdkMath.abs(t.getReal())));
407         final double dt = eqn.getTime().subtract(t).abs().getReal();
408         if (dt <= threshold) {
409             throw new NumberIsTooSmallException(LocalizedFormats.TOO_SMALL_INTEGRATION_INTERVAL,
410                                                 dt, threshold, false);
411         }
412     }
413 
414     /** Check if a reset occurred while last step was accepted.
415      * @return true if a reset occurred while last step was accepted
416      */
417     protected boolean resetOccurred() {
418         return resetOccurred;
419     }
420 
421     /** Set the current step size.
422      * @param stepSize step size to set
423      */
424     protected void setStepSize(final T stepSize) {
425         this.stepSize = stepSize;
426     }
427 
428     /** Get the current step size.
429      * @return current step size
430      */
431     protected T getStepSize() {
432         return stepSize;
433     }
434     /** Set current step start.
435      * @param stepStart step start
436      */
437     protected void setStepStart(final FieldODEStateAndDerivative<T> stepStart) {
438         this.stepStart = stepStart;
439     }
440 
441     /** Getcurrent step start.
442      * @return current step start
443      */
444     protected FieldODEStateAndDerivative<T> getStepStart() {
445         return stepStart;
446     }
447 
448     /** Set the last state flag.
449      * @param isLastStep if true, this step is the last one
450      */
451     protected void setIsLastStep(final boolean isLastStep) {
452         this.isLastStep = isLastStep;
453     }
454 
455     /** Check if this step is the last one.
456      * @return true if this step is the last one
457      */
458     protected boolean isLastStep() {
459         return isLastStep;
460     }
461 }