/*
 * Licensed to the Apache Software Foundation (ASF) under one or more
 * contributor license agreements.  See the NOTICE file distributed with
 * this work for additional information regarding copyright ownership.
 * The ASF licenses this file to You under the Apache License, Version 2.0
 * (the "License"); you may not use this file except in compliance with
 * the License.  You may obtain a copy of the License at
 *
 *      http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */

package org.apache.commons.math3.ode;

import java.io.Serializable;
import java.util.ArrayList;
import java.util.List;

import org.apache.commons.math3.exception.DimensionMismatchException;
import org.apache.commons.math3.exception.MathIllegalArgumentException;
import org.apache.commons.math3.exception.MaxCountExceededException;
import org.apache.commons.math3.exception.util.LocalizedFormats;
import org.apache.commons.math3.ode.sampling.StepHandler;
import org.apache.commons.math3.ode.sampling.StepInterpolator;
import org.apache.commons.math3.util.FastMath;

This class stores all information provided by an ODE integrator
during the integration process and build a continuous model of the
solution from this.
This class act as a step handler from the integrator point of view. It is called iteratively during the integration process and stores a copy of all steps information in a sorted collection for later use. Once the integration process is over, the user can use the setInterpolatedTime and getInterpolatedState to retrieve this information at any time. It is important to wait for the integration to be over before attempting to call setInterpolatedTime because some internal variables are set only once the last step has been handled.
This is useful for example if the main loop of the user
application should remain independent from the integration process
or if one needs to mimic the behaviour of an analytical model
despite a numerical model is used (i.e. one needs the ability to
get the model value at any time or to navigate through the
data).
If problem modeling is done with several separate
integration phases for contiguous intervals, the same
ContinuousOutputModel can be used as step handler for all
integration phases as long as they are performed in order and in
the same direction. As an example, one can extrapolate the
trajectory of a satellite with one model (i.e. one set of
differential equations) up to the beginning of a maneuver, use
another more complex model including thrusters modeling and
accurate attitude control during the maneuver, and revert to the
first model after the end of the maneuver. If the same continuous
output model handles the steps of all integration phases, the user
do not need to bother when the maneuver begins or ends, he has all
the data available in a transparent manner.
An important feature of this class is that it implements the
Serializable interface. This means that the result of
an integration can be serialized and reused later (if stored into a
persistent medium like a filesystem or a database) or elsewhere (if
sent to another application). Only the result of the integration is
stored, there is no reference to the integrated problem by
itself.
One should be aware that the amount of data stored in a ContinuousOutputModel instance can be important if the state vector is large, if the integration interval is long or if the steps are small (which can result from small tolerance settings in adaptive
step size integrators).
See Also:
StepHandler
StepInterpolator
Since:
1.2/**
 * This class stores all information provided by an ODE integrator
 * during the integration process and build a continuous model of the
 * solution from this.
 *
 * <p>This class act as a step handler from the integrator point of
 * view. It is called iteratively during the integration process and
 * stores a copy of all steps information in a sorted collection for
 * later use. Once the integration process is over, the user can use
 * the {@link #setInterpolatedTime setInterpolatedTime} and {@link
 * #getInterpolatedState getInterpolatedState} to retrieve this
 * information at any time. It is important to wait for the
 * integration to be over before attempting to call {@link
 * #setInterpolatedTime setInterpolatedTime} because some internal
 * variables are set only once the last step has been handled.</p>
 *
 * <p>This is useful for example if the main loop of the user
 * application should remain independent from the integration process
 * or if one needs to mimic the behaviour of an analytical model
 * despite a numerical model is used (i.e. one needs the ability to
 * get the model value at any time or to navigate through the
 * data).</p>
 *
 * <p>If problem modeling is done with several separate
 * integration phases for contiguous intervals, the same
 * ContinuousOutputModel can be used as step handler for all
 * integration phases as long as they are performed in order and in
 * the same direction. As an example, one can extrapolate the
 * trajectory of a satellite with one model (i.e. one set of
 * differential equations) up to the beginning of a maneuver, use
 * another more complex model including thrusters modeling and
 * accurate attitude control during the maneuver, and revert to the
 * first model after the end of the maneuver. If the same continuous
 * output model handles the steps of all integration phases, the user
 * do not need to bother when the maneuver begins or ends, he has all
 * the data available in a transparent manner.</p>
 *
 * <p>An important feature of this class is that it implements the
 * <code>Serializable</code> interface. This means that the result of
 * an integration can be serialized and reused later (if stored into a
 * persistent medium like a filesystem or a database) or elsewhere (if
 * sent to another application). Only the result of the integration is
 * stored, there is no reference to the integrated problem by
 * itself.</p>
 *
 * <p>One should be aware that the amount of data stored in a
 * ContinuousOutputModel instance can be important if the state vector
 * is large, if the integration interval is long or if the steps are
 * small (which can result from small tolerance settings in {@link
 * org.apache.commons.math3.ode.nonstiff.AdaptiveStepsizeIntegrator adaptive
 * step size integrators}).</p>
 *
 * @see StepHandler
 * @see StepInterpolator
 * @since 1.2
 */

public class ContinuousOutputModel
  implements StepHandler, Serializable {

    
Serializable version identifier /** Serializable version identifier */
    private static final long serialVersionUID = -1417964919405031606L;

    
Initial integration time. /** Initial integration time. */
    private double initialTime;

    
Final integration time. /** Final integration time. */
    private double finalTime;

    
Integration direction indicator. /** Integration direction indicator. */
    private boolean forward;

    
Current interpolator index. /** Current interpolator index. */
    private int index;

    
Steps table. /** Steps table. */
    private List<StepInterpolator> steps;

  
Simple constructor.
Build an empty continuous output model.
/** Simple constructor.
   * Build an empty continuous output model.
   */
  public ContinuousOutputModel() {
    steps = new ArrayList<StepInterpolator>();
    initialTime = Double.NaN;
    finalTime   = Double.NaN;
    forward     = true;
    index       = 0;
  }

  
Append another model at the end of the instance.
Params:
model – model to add at the end of the instance
Throws:
MathIllegalArgumentException – if the model to append is not
compatible with the instance (dimension of the state vector,
propagation direction, hole between the dates)
MaxCountExceededException – if the number of functions evaluations is exceeded
during step finalization/** Append another model at the end of the instance.
   * @param model model to add at the end of the instance
   * @exception MathIllegalArgumentException if the model to append is not
   * compatible with the instance (dimension of the state vector,
   * propagation direction, hole between the dates)
   * @exception MaxCountExceededException if the number of functions evaluations is exceeded
   * during step finalization
   */
  public void append(final ContinuousOutputModel model)
    throws MathIllegalArgumentException, MaxCountExceededException {

    if (model.steps.size() == 0) {
      return;
    }

    if (steps.size() == 0) {
      initialTime = model.initialTime;
      forward     = model.forward;
    } else {

      if (getInterpolatedState().length != model.getInterpolatedState().length) {
          throw new DimensionMismatchException(model.getInterpolatedState().length,
                                               getInterpolatedState().length);
      }

      if (forward ^ model.forward) {
          throw new MathIllegalArgumentException(LocalizedFormats.PROPAGATION_DIRECTION_MISMATCH);
      }

      final StepInterpolator lastInterpolator = steps.get(index);
      final double current  = lastInterpolator.getCurrentTime();
      final double previous = lastInterpolator.getPreviousTime();
      final double step = current - previous;
      final double gap = model.getInitialTime() - current;
      if (FastMath.abs(gap) > 1.0e-3 * FastMath.abs(step)) {
        throw new MathIllegalArgumentException(LocalizedFormats.HOLE_BETWEEN_MODELS_TIME_RANGES,
                                               FastMath.abs(gap));
      }

    }

    for (StepInterpolator interpolator : model.steps) {
      steps.add(interpolator.copy());
    }

    index = steps.size() - 1;
    finalTime = (steps.get(index)).getCurrentTime();

  }

  
{@inheritDoc} /** {@inheritDoc} */
  public void init(double t0, double[] y0, double t) {
    initialTime = Double.NaN;
    finalTime   = Double.NaN;
    forward     = true;
    index       = 0;
    steps.clear();
  }

  
Handle the last accepted step.
A copy of the information provided by the last step is stored in
the instance for later use.
Params:
interpolator – interpolator for the last accepted step.
isLast – true if the step is the last one
Throws:
MaxCountExceededException – if the number of functions evaluations is exceeded
during step finalization/** Handle the last accepted step.
   * A copy of the information provided by the last step is stored in
   * the instance for later use.
   * @param interpolator interpolator for the last accepted step.
   * @param isLast true if the step is the last one
   * @exception MaxCountExceededException if the number of functions evaluations is exceeded
   * during step finalization
   */
  public void handleStep(final StepInterpolator interpolator, final boolean isLast)
      throws MaxCountExceededException {

    if (steps.size() == 0) {
      initialTime = interpolator.getPreviousTime();
      forward     = interpolator.isForward();
    }

    steps.add(interpolator.copy());

    if (isLast) {
      finalTime = interpolator.getCurrentTime();
      index     = steps.size() - 1;
    }

  }

  
Get the initial integration time.
Returns:
initial integration time/**
   * Get the initial integration time.
   * @return initial integration time
   */
  public double getInitialTime() {
    return initialTime;
  }

  
Get the final integration time.
Returns:
final integration time/**
   * Get the final integration time.
   * @return final integration time
   */
  public double getFinalTime() {
    return finalTime;
  }

  
Get the time of the interpolated point. If setInterpolatedTime has not been called, it returns the final integration time. 
Returns:
interpolation point time/**
   * Get the time of the interpolated point.
   * If {@link #setInterpolatedTime} has not been called, it returns
   * the final integration time.
   * @return interpolation point time
   */
  public double getInterpolatedTime() {
    return steps.get(index).getInterpolatedTime();
  }

  
Set the time of the interpolated point.
This method should not be called before the
integration is over because some internal variables are set only
once the last step has been handled.
Setting the time outside of the integration interval is now
allowed, but should be used with care since the accuracy of the
interpolator will probably be very poor far from this interval.
This allowance has been added to simplify implementation of search
algorithms near the interval endpoints.
Note that each time this method is called, the internal arrays returned in getInterpolatedState(), getInterpolatedDerivatives() and getInterpolatedSecondaryState(int) will be overwritten. So if their content must be preserved
across several calls, user must copy them.
Params:
time – time of the interpolated point
See Also:
getInterpolatedState()
getInterpolatedDerivatives()
getInterpolatedSecondaryState(int)/** Set the time of the interpolated point.
   * <p>This method should <strong>not</strong> be called before the
   * integration is over because some internal variables are set only
   * once the last step has been handled.</p>
   * <p>Setting the time outside of the integration interval is now
   * allowed, but should be used with care since the accuracy of the
   * interpolator will probably be very poor far from this interval.
   * This allowance has been added to simplify implementation of search
   * algorithms near the interval endpoints.</p>
   * <p>Note that each time this method is called, the internal arrays
   * returned in {@link #getInterpolatedState()}, {@link
   * #getInterpolatedDerivatives()} and {@link #getInterpolatedSecondaryState(int)}
   * <em>will</em> be overwritten. So if their content must be preserved
   * across several calls, user must copy them.</p>
   * @param time time of the interpolated point
   * @see #getInterpolatedState()
   * @see #getInterpolatedDerivatives()
   * @see #getInterpolatedSecondaryState(int)
   */
  public void setInterpolatedTime(final double time) {

      // initialize the search with the complete steps table
      int iMin = 0;
      final StepInterpolator sMin = steps.get(iMin);
      double tMin = 0.5 * (sMin.getPreviousTime() + sMin.getCurrentTime());

      int iMax = steps.size() - 1;
      final StepInterpolator sMax = steps.get(iMax);
      double tMax = 0.5 * (sMax.getPreviousTime() + sMax.getCurrentTime());

      // handle points outside of the integration interval
      // or in the first and last step
      if (locatePoint(time, sMin) <= 0) {
        index = iMin;
        sMin.setInterpolatedTime(time);
        return;
      }
      if (locatePoint(time, sMax) >= 0) {
        index = iMax;
        sMax.setInterpolatedTime(time);
        return;
      }

      // reduction of the table slice size
      while (iMax - iMin > 5) {

        // use the last estimated index as the splitting index
        final StepInterpolator si = steps.get(index);
        final int location = locatePoint(time, si);
        if (location < 0) {
          iMax = index;
          tMax = 0.5 * (si.getPreviousTime() + si.getCurrentTime());
        } else if (location > 0) {
          iMin = index;
          tMin = 0.5 * (si.getPreviousTime() + si.getCurrentTime());
        } else {
          // we have found the target step, no need to continue searching
          si.setInterpolatedTime(time);
          return;
        }

        // compute a new estimate of the index in the reduced table slice
        final int iMed = (iMin + iMax) / 2;
        final StepInterpolator sMed = steps.get(iMed);
        final double tMed = 0.5 * (sMed.getPreviousTime() + sMed.getCurrentTime());

        if ((FastMath.abs(tMed - tMin) < 1e-6) || (FastMath.abs(tMax - tMed) < 1e-6)) {
          // too close to the bounds, we estimate using a simple dichotomy
          index = iMed;
        } else {
          // estimate the index using a reverse quadratic polynom
          // (reverse means we have i = P(t), thus allowing to simply
          // compute index = P(time) rather than solving a quadratic equation)
          final double d12 = tMax - tMed;
          final double d23 = tMed - tMin;
          final double d13 = tMax - tMin;
          final double dt1 = time - tMax;
          final double dt2 = time - tMed;
          final double dt3 = time - tMin;
          final double iLagrange = ((dt2 * dt3 * d23) * iMax -
                                    (dt1 * dt3 * d13) * iMed +
                                    (dt1 * dt2 * d12) * iMin) /
                                   (d12 * d23 * d13);
          index = (int) FastMath.rint(iLagrange);
        }

        // force the next size reduction to be at least one tenth
        final int low  = FastMath.max(iMin + 1, (9 * iMin + iMax) / 10);
        final int high = FastMath.min(iMax - 1, (iMin + 9 * iMax) / 10);
        if (index < low) {
          index = low;
        } else if (index > high) {
          index = high;
        }

      }

      // now the table slice is very small, we perform an iterative search
      index = iMin;
      while ((index <= iMax) && (locatePoint(time, steps.get(index)) > 0)) {
        ++index;
      }

      steps.get(index).setInterpolatedTime(time);

  }

  
Get the state vector of the interpolated point.
The returned vector is a reference to a reused array, so it should not be modified and it should be copied if it needs to be preserved across several calls to the associated setInterpolatedTime(double) method.
Throws:
MaxCountExceededException – if the number of functions evaluations is exceeded
See Also:
setInterpolatedTime(double)
getInterpolatedDerivatives()
getInterpolatedSecondaryState(int)
getInterpolatedSecondaryDerivatives(int)
Returns:
state vector at time getInterpolatedTime/**
   * Get the state vector of the interpolated point.
   * <p>The returned vector is a reference to a reused array, so
   * it should not be modified and it should be copied if it needs
   * to be preserved across several calls to the associated
   * {@link #setInterpolatedTime(double)} method.</p>
   * @return state vector at time {@link #getInterpolatedTime}
   * @exception MaxCountExceededException if the number of functions evaluations is exceeded
   * @see #setInterpolatedTime(double)
   * @see #getInterpolatedDerivatives()
   * @see #getInterpolatedSecondaryState(int)
   * @see #getInterpolatedSecondaryDerivatives(int)
   */
  public double[] getInterpolatedState() throws MaxCountExceededException {
    return steps.get(index).getInterpolatedState();
  }

  
Get the derivatives of the state vector of the interpolated point.
The returned vector is a reference to a reused array, so it should not be modified and it should be copied if it needs to be preserved across several calls to the associated setInterpolatedTime(double) method.
Throws:
MaxCountExceededException – if the number of functions evaluations is exceeded
See Also:
setInterpolatedTime(double)
getInterpolatedState()
getInterpolatedSecondaryState(int)
getInterpolatedSecondaryDerivatives(int)
Returns:
derivatives of the state vector at time getInterpolatedTime
Since:
3.4/**
   * Get the derivatives of the state vector of the interpolated point.
   * <p>The returned vector is a reference to a reused array, so
   * it should not be modified and it should be copied if it needs
   * to be preserved across several calls to the associated
   * {@link #setInterpolatedTime(double)} method.</p>
   * @return derivatives of the state vector at time {@link #getInterpolatedTime}
   * @exception MaxCountExceededException if the number of functions evaluations is exceeded
   * @see #setInterpolatedTime(double)
   * @see #getInterpolatedState()
   * @see #getInterpolatedSecondaryState(int)
   * @see #getInterpolatedSecondaryDerivatives(int)
   * @since 3.4
   */
  public double[] getInterpolatedDerivatives() throws MaxCountExceededException {
    return steps.get(index).getInterpolatedDerivatives();
  }

  
Get the interpolated secondary state corresponding to the secondary equations.
The returned vector is a reference to a reused array, so it should not be modified and it should be copied if it needs to be preserved across several calls to the associated setInterpolatedTime(double) method.
Params:
secondaryStateIndex – index of the secondary set, as returned by 
ExpandableStatefulODE.addSecondaryEquations(SecondaryEquations)
Throws:
MaxCountExceededException – if the number of functions evaluations is exceeded
See Also:
setInterpolatedTime(double)
getInterpolatedState()
getInterpolatedDerivatives()
getInterpolatedSecondaryDerivatives(int)
Returns:
interpolated secondary state at the current interpolation date
Since:
3.2/** Get the interpolated secondary state corresponding to the secondary equations.
   * <p>The returned vector is a reference to a reused array, so
   * it should not be modified and it should be copied if it needs
   * to be preserved across several calls to the associated
   * {@link #setInterpolatedTime(double)} method.</p>
   * @param secondaryStateIndex index of the secondary set, as returned by {@link
   * org.apache.commons.math3.ode.ExpandableStatefulODE#addSecondaryEquations(
   * org.apache.commons.math3.ode.SecondaryEquations)
   * ExpandableStatefulODE.addSecondaryEquations(SecondaryEquations)}
   * @return interpolated secondary state at the current interpolation date
   * @see #setInterpolatedTime(double)
   * @see #getInterpolatedState()
   * @see #getInterpolatedDerivatives()
   * @see #getInterpolatedSecondaryDerivatives(int)
   * @since 3.2
   * @exception MaxCountExceededException if the number of functions evaluations is exceeded
   */
  public double[] getInterpolatedSecondaryState(final int secondaryStateIndex)
    throws MaxCountExceededException {
    return steps.get(index).getInterpolatedSecondaryState(secondaryStateIndex);
  }

  
Get the interpolated secondary derivatives corresponding to the secondary equations.
The returned vector is a reference to a reused array, so it should not be modified and it should be copied if it needs to be preserved across several calls to the associated setInterpolatedTime(double) method.
Params:
secondaryStateIndex – index of the secondary set, as returned by 
ExpandableStatefulODE.addSecondaryEquations(SecondaryEquations)
Throws:
MaxCountExceededException – if the number of functions evaluations is exceeded
See Also:
setInterpolatedTime(double)
getInterpolatedState()
getInterpolatedDerivatives()
getInterpolatedSecondaryState(int)
Returns:
interpolated secondary derivatives at the current interpolation date
Since:
3.4/** Get the interpolated secondary derivatives corresponding to the secondary equations.
   * <p>The returned vector is a reference to a reused array, so
   * it should not be modified and it should be copied if it needs
   * to be preserved across several calls to the associated
   * {@link #setInterpolatedTime(double)} method.</p>
   * @param secondaryStateIndex index of the secondary set, as returned by {@link
   * org.apache.commons.math3.ode.ExpandableStatefulODE#addSecondaryEquations(
   * org.apache.commons.math3.ode.SecondaryEquations)
   * ExpandableStatefulODE.addSecondaryEquations(SecondaryEquations)}
   * @return interpolated secondary derivatives at the current interpolation date
   * @see #setInterpolatedTime(double)
   * @see #getInterpolatedState()
   * @see #getInterpolatedDerivatives()
   * @see #getInterpolatedSecondaryState(int)
   * @since 3.4
   * @exception MaxCountExceededException if the number of functions evaluations is exceeded
   */
  public double[] getInterpolatedSecondaryDerivatives(final int secondaryStateIndex)
    throws MaxCountExceededException {
    return steps.get(index).getInterpolatedSecondaryDerivatives(secondaryStateIndex);
  }

  
Compare a step interval and a double.
Params:
time – point to locate
interval – step interval
Returns:
-1 if the double is before the interval, 0 if it is in
the interval, and +1 if it is after the interval, according to
the interval direction/** Compare a step interval and a double.
   * @param time point to locate
   * @param interval step interval
   * @return -1 if the double is before the interval, 0 if it is in
   * the interval, and +1 if it is after the interval, according to
   * the interval direction
   */
  private int locatePoint(final double time, final StepInterpolator interval) {
    if (forward) {
      if (time < interval.getPreviousTime()) {
        return -1;
      } else if (time > interval.getCurrentTime()) {
        return +1;
      } else {
        return 0;
      }
    }
    if (time > interval.getPreviousTime()) {
      return -1;
    } else if (time < interval.getCurrentTime()) {
      return +1;
    } else {
      return 0;
    }
  }

}