AbstractFieldIntegrator.java

  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.  */

  18. package org.apache.commons.math4.legacy.ode;

  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;

  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;

  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> {

  52.     /** Default relative accuracy. */
  53.     private static final double DEFAULT_RELATIVE_ACCURACY = 1e-14;

  55.     /** Default function value accuracy. */
  56.     private static final double DEFAULT_FUNCTION_VALUE_ACCURACY = 1e-15;

  58.     /** Step handler. */
  59.     private Collection<FieldStepHandler<T>> stepHandlers;

  61.     /** Current step start. */
  62.     private FieldODEStateAndDerivative<T> stepStart;

  64.     /** Current stepsize. */
  65.     private T stepSize;

  67.     /** Indicator for last step. */
  68.     private boolean isLastStep;

  70.     /** Indicator that a state or derivative reset was triggered by some event. */
  71.     private boolean resetOccurred;

  73.     /** Field to which the time and state vector elements belong. */
  74.     private final Field<T> field;

  76.     /** Events states. */
  77.     private Collection<FieldEventState<T>> eventsStates;

  79.     /** Initialization indicator of events states. */
  80.     private boolean statesInitialized;

  82.     /** Name of the method. */
  83.     private final String name;

  85.     /** Counter for number of evaluations. */
  86.     private IntegerSequence.Incrementor evaluations;

  88.     /** Differential equations to integrate. */
  89.     private transient FieldExpandableODE<T> equations;

  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.     }

  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.     }

  113.     /** {@inheritDoc} */
  114.     @Override
  115.     public String getName() {
  116.         return name;
  117.     }

  119.     /** {@inheritDoc} */
  120.     @Override
  121.     public void addStepHandler(final FieldStepHandler<T> handler) {
  122.         stepHandlers.add(handler);
  123.     }

  125.     /** {@inheritDoc} */
  126.     @Override
  127.     public Collection<FieldStepHandler<T>> getStepHandlers() {
  128.         return Collections.unmodifiableCollection(stepHandlers);
  129.     }

  131.     /** {@inheritDoc} */
  132.     @Override
  133.     public void clearStepHandlers() {
  134.         stepHandlers.clear();
  135.     }

  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.     }

  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.     }

  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.     }

  172.     /** {@inheritDoc} */
  173.     @Override
  174.     public void clearEventHandlers() {
  175.         eventsStates.clear();
  176.     }

  178.     /** {@inheritDoc} */
  179.     @Override
  180.     public FieldODEStateAndDerivative<T> getCurrentStepStart() {
  181.         return stepStart;
  182.     }

  184.     /** {@inheritDoc} */
  185.     @Override
  186.     public T getCurrentSignedStepsize() {
  187.         return stepSize;
  188.     }

  190.     /** {@inheritDoc} */
  191.     @Override
  192.     public void setMaxEvaluations(int maxEvaluations) {
  193.         evaluations = evaluations.withMaximalCount((maxEvaluations < 0) ? Integer.MAX_VALUE : maxEvaluations);
  194.     }

  196.     /** {@inheritDoc} */
  197.     @Override
  198.     public int getMaxEvaluations() {
  199.         return evaluations.getMaximalCount();
  200.     }

  202.     /** {@inheritDoc} */
  203.     @Override
  204.     public int getEvaluations() {
  205.         return evaluations.getCount();
  206.     }

  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) {

  218.         this.equations = eqn;
  219.         evaluations    = evaluations.withStart(0);

  221.         // initialize ODE
  222.         eqn.init(t0, y0, t);

  224.         // set up derivatives of initial state
  225.         final T[] y0Dot = computeDerivatives(t0, y0);
  226.         final FieldODEStateAndDerivative<T> state0 = new FieldODEStateAndDerivative<>(t0, y0, y0Dot);

  228.         // initialize event handlers
  229.         for (final FieldEventState<T> state : eventsStates) {
  230.             state.getEventHandler().init(state0, t);
  231.         }

  233.         // initialize step handlers
  234.         for (FieldStepHandler<T> handler : stepHandlers) {
  235.             handler.init(state0, t);
  236.         }

  238.         setStateInitialized(false);

  240.         return state0;
  241.     }

  243.     /** Get the differential equations to integrate.
  244.      * @return differential equations to integrate
  245.      */
  246.     protected FieldExpandableODE<T> getEquations() {
  247.         return equations;
  248.     }

  250.     /** Get the evaluations counter.
  251.      * @return evaluations counter
  252.      */
  253.     protected IntegerSequence.Incrementor getEvaluationsCounter() {
  254.         return evaluations;
  255.     }

  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.     }

  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.     }

  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 {

  296.             FieldODEStateAndDerivative<T> previousState = interpolator.getGlobalPreviousState();
  297.             final FieldODEStateAndDerivative<T> currentState = interpolator.getGlobalCurrentState();

  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.             }

  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>>() {

  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.             });

  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.             }

  325.             AbstractFieldStepInterpolator<T> restricted = interpolator;
  326.             while (!occurringEvents.isEmpty()) {

  328.                 // handle the chronologically first event
  329.                 final Iterator<FieldEventState<T>> iterator = occurringEvents.iterator();
  330.                 final FieldEventState<T> currentEvent = iterator.next();
  331.                 iterator.remove();

  333.                 // get state at event time
  334.                 final FieldODEStateAndDerivative<T> eventState = restricted.getInterpolatedState(currentEvent.getEventTime());

  336.                 // restrict the interpolator to the first part of the step, up to the event
  337.                 restricted = restricted.restrictStep(previousState, eventState);

  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.                 }

  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.                 }

  350.                 if (isLastStep) {
  351.                     // the event asked to stop integration
  352.                     return eventState;
  353.                 }

  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.                 }

  369.                 // prepare handling of the remaining part of the step
  370.                 previousState = eventState;
  371.                 restricted = restricted.restrictStep(eventState, currentState);

  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.             }

  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());

  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.             }

  392.             return currentState;
  393.     }

  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 {

  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.     }

  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.     }

  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.     }

  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.     }

  441.     /** Getcurrent step start.
  442.      * @return current step start
  443.      */
  444.     protected FieldODEStateAndDerivative<T> getStepStart() {
  445.         return stepStart;
  446.     }

  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.     }

  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. }
Created with JaCoCo 0.8.10.202304240956Code Coverage Report for Apache Commons Math 4.0-SNAPSHOT