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