/*
 * Copyright (C) 2013 - present by OpenGamma Inc. and the OpenGamma group of companies
 *
 * Please see distribution for license.
 */
package com.opengamma.strata.math.impl.interpolation;

import java.util.Arrays;

import com.opengamma.strata.collect.array.DoubleArray;
import com.opengamma.strata.collect.array.DoubleMatrix;

/**
 * Hermite interpolation is determined if one specifies first derivatives for a cubic
 * interpolant and first and second derivatives for a quintic interpolant.
 */
public class HermiteCoefficientsProvider {

  /**
   * @param values (yValues_i)
   * @param intervals (xValues_{i+1} - xValues_{i})
   * @param slopes (yValues_{i+1} - yValues_{i})/(xValues_{i+1} - xValues_{i})
   * @param first First derivatives at xValues_i
   * @return Coefficient matrix whose i-th row vector is { a_n, a_{n-1}, ...} for the i-th interval,
   *   where a_n, a_{n-1},... are coefficients of f(x) = a_n (x-x_i)^n + a_{n-1} (x-x_i)^{n-1} + .... with n=3
   */
  public double[][] solve(double[] values, double[] intervals, double[] slopes, double[] first) {
    int nInt = intervals.length;
    double[][] res = new double[nInt][4];
    for (int i = 0; i < nInt; ++i) {
      Arrays.fill(res[i], 0.);
    }

    for (int i = 0; i < nInt; ++i) {
      res[i][3] = values[i];
      res[i][2] = first[i];
      res[i][1] = (3. * slopes[i] - first[i + 1] - 2. * first[i]) / intervals[i];
      res[i][0] = -(2. * slopes[i] - first[i + 1] - first[i]) / intervals[i] / intervals[i];
    }

    return res;
  }

  /**
   * @param values Y values of data
   * @param intervals (xValues_{i+1} - xValues_{i})
   * @param slopes (yValues_{i+1} - yValues_{i})/(xValues_{i+1} - xValues_{i})
   * @param slopeSensitivity Derivative values of slope with respect to yValues
   * @param firstWithSensitivity First derivative values at xValues_i and their yValues dependencies
   * @return Coefficient matrix and its node dependencies
   */
  public DoubleMatrix[] solveWithSensitivity(
      double[] values,
      double[] intervals,
      double[] slopes,
      double[][] slopeSensitivity,
      DoubleArray[] firstWithSensitivity) {

    int nData = values.length;
    double[] first = firstWithSensitivity[0].toArray();
    DoubleMatrix[] res = new DoubleMatrix[nData];

    double[][] coef = solve(values, intervals, slopes, first);
    res[0] = DoubleMatrix.copyOf(coef);

    for (int i = 0; i < nData - 1; ++i) {
      double[][] coefSense = new double[4][nData];
      Arrays.fill(coefSense[3], 0.);
      coefSense[3][i] = 1.;
      for (int k = 0; k < nData; ++k) {
        coefSense[0][k] =
            -(2. * slopeSensitivity[i][k] - firstWithSensitivity[i + 2].get(k) - firstWithSensitivity[i + 1].get(k)) /
                intervals[i] / intervals[i];
        coefSense[1][k] =
            (3. * slopeSensitivity[i][k] - firstWithSensitivity[i + 2].get(k) - 2. * firstWithSensitivity[i + 1].get(k)) /
                intervals[i];
        coefSense[2][k] = firstWithSensitivity[i + 1].get(k);
      }
      res[i + 1] = DoubleMatrix.copyOf(coefSense);
    }

    return res;
  }

  /**
   * @param values (yValues_i)
   * @param intervals (xValues_{i+1} - xValues_{i})
   * @param slopes (yValues_{i+1} - yValues_{i})/(xValues_{i+1} - xValues_{i})
   * @param first First derivatives at xValues_i
   * @param second Second derivatives at xValues_i
   * @return Coefficient matrix whose i-th row vector is { a_n, a_{n-1}, ...} for the i-th interval,
   *   where a_n, a_{n-1},... are coefficients of f(x) = a_n (x-x_i)^n + a_{n-1} (x-x_i)^{n-1} + .... with n=5
   */
  public double[][] solve(double[] values, double[] intervals, double[] slopes, double[] first, double[] second) {
    int nInt = intervals.length;
    double[][] res = new double[nInt][6];
    for (int i = 0; i < nInt; ++i) {
      Arrays.fill(res[i], 0.);
    }

    for (int i = 0; i < nInt; ++i) {
      res[i][5] = values[i];
      res[i][4] = first[i];
      res[i][3] = 0.5 * second[i];
      res[i][2] = 0.5 * (second[i + 1] - 3. * second[i]) / intervals[i] +
          2. * (5. * slopes[i] - 3. * first[i] - 2. * first[i + 1]) / intervals[i] / intervals[i];
      res[i][1] = 0.5 * (3. * second[i] - 2. * second[i + 1]) / intervals[i] / intervals[i] +
          (8. * first[i] + 7. * first[i + 1] - 15. * slopes[i]) / intervals[i] / intervals[i] / intervals[i];
      res[i][0] = 0.5 * (second[i + 1] - second[i]) / intervals[i] / intervals[i] / intervals[i] +
          3. * (2. * slopes[i] - first[i + 1] - first[i]) / intervals[i] / intervals[i] / intervals[i] / intervals[i];
    }

    return res;
  }

  /**
   *
   * @param values (yValues_i)
   * @param intervals (xValues_{i+1} - xValues_{i})
   * @param slopes (yValues_{i+1} - yValues_{i})/(xValues_{i+1} - xValues_{i})
   * @param slopeSensitivity Derivative values of slope with respect to yValues
   * @param firstWithSensitivity First derivative values at xValues_i and their yValues dependencies
   * @param secondWithSensitivity Second derivative values at xValues_i and their yValues dependencies
   * @return Coefficient matrix and its node dependencies
   */
  public DoubleMatrix[] solveWithSensitivity(
      double[] values,
      double[] intervals,
      double[] slopes,
      double[][] slopeSensitivity,
      DoubleArray[] firstWithSensitivity,
      DoubleArray[] secondWithSensitivity) {

    int nData = values.length;
    double[] first = firstWithSensitivity[0].toArray();
    double[] second = secondWithSensitivity[0].toArray();
    DoubleMatrix[] res = new DoubleMatrix[nData];

    double[][] coef = solve(values, intervals, slopes, first, second);
    res[0] = DoubleMatrix.copyOf(coef);

    for (int i = 0; i < nData - 1; ++i) {
      double interval = intervals[i];
      double[][] coefSense = new double[6][nData];
      Arrays.fill(coefSense[5], 0.);
      coefSense[5][i] = 1.;
      for (int k = 0; k < nData; ++k) {
        double cs0b = 2d * slopeSensitivity[i][k] - firstWithSensitivity[i + 2].get(k) - firstWithSensitivity[i + 1].get(k);
        double cs0a = secondWithSensitivity[i + 2].get(k) - secondWithSensitivity[i + 1].get(k);
        coefSense[0][k] = 0.5 * cs0a / interval / interval / interval + 3d * cs0b / interval / interval / interval / interval;

        double cs1a = 3d * secondWithSensitivity[i + 1].get(k) - 2d * secondWithSensitivity[i + 2].get(k);
        double cs1b =
            8d * firstWithSensitivity[i + 1].get(k) + 7d * firstWithSensitivity[i + 2].get(k) - 15d * slopeSensitivity[i][k];
        coefSense[1][k] = 0.5 * cs1a / interval / interval + cs1b / interval / interval / interval;

        double cs2a = secondWithSensitivity[i + 2].get(k) - 3d * secondWithSensitivity[i + 1].get(k);
        double cs2b = 5d * slopeSensitivity[i][k] - 3d * firstWithSensitivity[i + 1].get(k) - 2d * firstWithSensitivity[i + 2].get(k);
        coefSense[2][k] = 0.5 * cs2a / interval + 2. * cs2b / interval / interval;

        coefSense[3][k] = 0.5 * secondWithSensitivity[i + 1].get(k);
        coefSense[4][k] = firstWithSensitivity[i + 1].get(k);
      }
      res[i + 1] = DoubleMatrix.copyOf(coefSense);
    }

    return res;
  }

  /**
   * @param xValues The x values
   * @return Intervals of xValues, ( xValues_{i+1} - xValues_i )
   */
  public double[] intervalsCalculator(double[] xValues) {

    int nDataPts = xValues.length;
    double[] intervals = new double[nDataPts - 1];

    for (int i = 0; i < nDataPts - 1; ++i) {
      intervals[i] = xValues[i + 1] - xValues[i];
    }

    return intervals;
  }

  /**
   * @param yValues Y values of data
   * @param intervals Intervals of x data
   * @return ( yValues_{i+1} - yValues_i )/( xValues_{i+1} - xValues_i )
   */
  public double[] slopesCalculator(double[] yValues, double[] intervals) {

    int nDataPts = yValues.length;
    double[] slopes = new double[nDataPts - 1];

    for (int i = 0; i < nDataPts - 1; ++i) {
      slopes[i] = (yValues[i + 1] - yValues[i]) / intervals[i];
    }

    return slopes;
  }

  /**
   * Derivative values of slopes_i with respect to yValues_j, s_{ij}.
   * @param intervals Intervals of x data
   * @return The matrix s_{ij}
   */
  public double[][] slopeSensitivityCalculator(double[] intervals) {
    int nDataPts = intervals.length + 1;
    double[][] res = new double[nDataPts - 1][nDataPts];

    for (int i = 0; i < nDataPts - 1; ++i) {
      Arrays.fill(res[i], 0.);
      res[i][i] = -1. / intervals[i];
      res[i][i + 1] = 1. / intervals[i];
    }
    return res;
  }

  /**
   * @param ints1 The first interval
   * @param ints2 The second interval
   * @param slope1 The first gradient
   * @param slope2 The second gradient
   * @return Value of derivative at each endpoint
   */
  public double endpointDerivatives(double ints1, double ints2, double slope1, double slope2) {
    double val = (2. * ints1 + ints2) * slope1 / (ints1 + ints2) - ints1 * slope2 / (ints1 + ints2);

    if (Math.signum(val) != Math.signum(slope1)) {
      return 0.;
    }
    if (Math.signum(slope1) != Math.signum(slope2) && Math.abs(val) > 3. * Math.abs(slope1)) {
      return 3. * slope1;
    }
    return val;
  }
}
