/*
    * 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.math4.legacy.ode.sampling;
  
  import java.io.IOException;
  import java.io.ObjectInput;
  import java.io.ObjectOutput;
  import java.util.Arrays;
  
  import org.apache.commons.math4.legacy.exception.MaxCountExceededException;
  import org.apache.commons.math4.legacy.linear.Array2DRowRealMatrix;
  import org.apache.commons.math4.legacy.ode.EquationsMapper;
  import org.apache.commons.math4.core.jdkmath.JdkMath;
  
  /**
   * This class implements an interpolator for integrators using Nordsieck representation.
   *
   * <p>This interpolator computes dense output around the current point.
   * The interpolation equation is based on Taylor series formulas.
   *
   * @see org.apache.commons.math4.legacy.ode.nonstiff.AdamsBashforthIntegrator
   * @see org.apache.commons.math4.legacy.ode.nonstiff.AdamsMoultonIntegrator
   * @since 2.0
   */
  
  public class NordsieckStepInterpolator extends AbstractStepInterpolator {
  
      /** Serializable version identifier. */
      private static final long serialVersionUID = -7179861704951334960L;
  
      /** State variation. */
      protected double[] stateVariation;
  
      /** Step size used in the first scaled derivative and Nordsieck vector. */
      private double scalingH;
  
      /** Reference time for all arrays.
       * <p>Sometimes, the reference time is the same as previousTime,
       * sometimes it is the same as currentTime, so we use a separate
       * field to avoid any confusion.
       * </p>
       */
      private double referenceTime;
  
      /** First scaled derivative. */
      private double[] scaled;
  
      /** Nordsieck vector. */
      private Array2DRowRealMatrix nordsieck;
  
      /** Simple constructor.
       * This constructor builds an instance that is not usable yet, the
       * {@link AbstractStepInterpolator#reinitialize} method should be called
       * before using the instance in order to initialize the internal arrays. This
       * constructor is used only in order to delay the initialization in
       * some cases.
       */
      public NordsieckStepInterpolator() {
      }
  
      /** Copy constructor.
       * @param interpolator interpolator to copy from. The copy is a deep
       * copy: its arrays are separated from the original arrays of the
       * instance
       */
      public NordsieckStepInterpolator(final NordsieckStepInterpolator interpolator) {
          super(interpolator);
          scalingH      = interpolator.scalingH;
          referenceTime = interpolator.referenceTime;
          if (interpolator.scaled != null) {
              scaled = interpolator.scaled.clone();
          }
          if (interpolator.nordsieck != null) {
              nordsieck = new Array2DRowRealMatrix(interpolator.nordsieck.getDataRef(), true);
          }
          if (interpolator.stateVariation != null) {
              stateVariation = interpolator.stateVariation.clone();
          }
      }
  
      /** {@inheritDoc} */
      @Override
      protected StepInterpolator doCopy() {
          return new NordsieckStepInterpolator(this);
     }
 
     /** Reinitialize the instance.
      * <p>Beware that all arrays <em>must</em> be references to integrator
      * arrays, in order to ensure proper update without copy.</p>
      * @param y reference to the integrator array holding the state at
      * the end of the step
      * @param forward integration direction indicator
      * @param primaryMapper equations mapper for the primary equations set
      * @param secondaryMappers equations mappers for the secondary equations sets
      */
     @Override
     public void reinitialize(final double[] y, final boolean forward,
                              final EquationsMapper primaryMapper,
                              final EquationsMapper[] secondaryMappers) {
         super.reinitialize(y, forward, primaryMapper, secondaryMappers);
         stateVariation = new double[y.length];
     }
 
     /** Reinitialize the instance.
      * <p>Beware that all arrays <em>must</em> be references to integrator
      * arrays, in order to ensure proper update without copy.</p>
      * @param time time at which all arrays are defined
      * @param stepSize step size used in the scaled and Nordsieck arrays
      * @param scaledDerivative reference to the integrator array holding the first
      * scaled derivative
      * @param nordsieckVector reference to the integrator matrix holding the
      * Nordsieck vector
      */
     public void reinitialize(final double time, final double stepSize,
                              final double[] scaledDerivative,
                              final Array2DRowRealMatrix nordsieckVector) {
         this.referenceTime = time;
         this.scalingH      = stepSize;
         this.scaled        = scaledDerivative;
         this.nordsieck     = nordsieckVector;
 
         // make sure the state and derivatives will depend on the new arrays
         setInterpolatedTime(getInterpolatedTime());
     }
 
     /** Rescale the instance.
      * <p>Since the scaled and Nordsieck arrays are shared with the caller,
      * this method has the side effect of rescaling this arrays in the caller too.</p>
      * @param stepSize new step size to use in the scaled and Nordsieck arrays
      */
     public void rescale(final double stepSize) {
 
         final double ratio = stepSize / scalingH;
         for (int i = 0; i < scaled.length; ++i) {
             scaled[i] *= ratio;
         }
 
         final double[][] nData = nordsieck.getDataRef();
         double power = ratio;
         for (int i = 0; i < nData.length; ++i) {
             power *= ratio;
             final double[] nDataI = nData[i];
             for (int j = 0; j < nDataI.length; ++j) {
                 nDataI[j] *= power;
             }
         }
 
         scalingH = stepSize;
     }
 
     /**
      * Get the state vector variation from current to interpolated state.
      * <p>This method is aimed at computing y(t<sub>interpolation</sub>)
      * -y(t<sub>current</sub>) accurately by avoiding the cancellation errors
      * that would occur if the subtraction were performed explicitly.</p>
      * <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.</p>
      * @return state vector at time {@link #getInterpolatedTime}
      * @see #getInterpolatedDerivatives()
      * @exception MaxCountExceededException if the number of functions evaluations is exceeded
      */
     public double[] getInterpolatedStateVariation() throws MaxCountExceededException {
         // compute and ignore interpolated state
         // to make sure state variation is computed as a side effect
         getInterpolatedState();
         return stateVariation;
     }
 
     /** {@inheritDoc} */
     @Override
     protected void computeInterpolatedStateAndDerivatives(final double theta, final double oneMinusThetaH) {
 
         final double x = interpolatedTime - referenceTime;
         final double normalizedAbscissa = x / scalingH;
 
         Arrays.fill(stateVariation, 0.0);
         Arrays.fill(interpolatedDerivatives, 0.0);
 
         // apply Taylor formula from high order to low order,
         // for the sake of numerical accuracy
         final double[][] nData = nordsieck.getDataRef();
         for (int i = nData.length - 1; i >= 0; --i) {
             final int order = i + 2;
             final double[] nDataI = nData[i];
             final double power = JdkMath.pow(normalizedAbscissa, order);
             for (int j = 0; j < nDataI.length; ++j) {
                 final double d = nDataI[j] * power;
                 stateVariation[j]          += d;
                 interpolatedDerivatives[j] += order * d;
             }
         }
 
         for (int j = 0; j < currentState.length; ++j) {
             stateVariation[j] += scaled[j] * normalizedAbscissa;
             interpolatedState[j] = currentState[j] + stateVariation[j];
             interpolatedDerivatives[j] =
                 (interpolatedDerivatives[j] + scaled[j] * normalizedAbscissa) / x;
         }
     }
 
     /** {@inheritDoc} */
     @Override
     public void writeExternal(final ObjectOutput out)
         throws IOException {
 
         // save the state of the base class
         writeBaseExternal(out);
 
         // save the local attributes
         out.writeDouble(scalingH);
         out.writeDouble(referenceTime);
 
         final int n = (currentState == null) ? -1 : currentState.length;
         if (scaled == null) {
             out.writeBoolean(false);
         } else {
             out.writeBoolean(true);
             for (int j = 0; j < n; ++j) {
                 out.writeDouble(scaled[j]);
             }
         }
 
         if (nordsieck == null) {
             out.writeBoolean(false);
         } else {
             out.writeBoolean(true);
             out.writeObject(nordsieck);
         }
 
         // we don't save state variation, it will be recomputed
     }
 
     /** {@inheritDoc} */
     @Override
     public void readExternal(final ObjectInput in)
         throws IOException, ClassNotFoundException {
 
         // read the base class
         final double t = readBaseExternal(in);
 
         // read the local attributes
         scalingH      = in.readDouble();
         referenceTime = in.readDouble();
 
         final int n = (currentState == null) ? -1 : currentState.length;
         final boolean hasScaled = in.readBoolean();
         if (hasScaled) {
             scaled = new double[n];
             for (int j = 0; j < n; ++j) {
                 scaled[j] = in.readDouble();
             }
         } else {
             scaled = null;
         }
 
         final boolean hasNordsieck = in.readBoolean();
         if (hasNordsieck) {
             nordsieck = (Array2DRowRealMatrix) in.readObject();
         } else {
             nordsieck = null;
         }
 
         if (hasScaled && hasNordsieck) {
             // we can now set the interpolated time and state
             stateVariation = new double[n];
             setInterpolatedTime(t);
         } else {
             stateVariation = null;
         }
     }
 }