/**
 * @file Interpolation.js
 *
 * Permission is hereby granted to any person obtaining a copy of this software 
 * and associated documentation files (the "Software"), to use it for personal 
 * or non-commercial purposes, with the following restrictions:
 *
 * 1. **No Copying or Redistribution**: The Software or any of its parts may not 
 *    be copied, merged, distributed, sublicensed, or sold without explicit 
 *    prior written permission from the author.
 * 
 * 2. **Commercial Use**: Any use of the Software for commercial purposes requires 
 *    a valid license, obtainable only with the explicit consent of the author.
 * 
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR 
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, 
 * FITNESS FOR A PARTICULAR PURPOSE, AND NONINFRINGEMENT. IN NO EVENT SHALL THE 
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES, OR OTHER 
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT, OR OTHERWISE, ARISING FROM, 
 * OUT OF, OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE 
 * SOFTWARE.
 *
 * Ownership of this code remains solely with the original author. Unauthorized 
 * use of this Software is strictly prohibited.
 *
 * Author:
 * - Rene De Ren
 * Email:
 * - rene@thegoldenbasket.nl
 *
/*
    Interpolate using cubic Hermite splines. The breakpoints in arrays xbp and ybp are assumed to be sorted.
    Evaluate the function in all points of the array xeval.
    Methods:
        "Linear"                yuck
        "FiniteDifference"      classic cubic interpolation, no tension parameter
        "Cardinal"              cubic cardinal splines, uses tension parameter which must be between [0,1]
        "FritschCarlson"        monotonic - tangents are first initialized, then adjusted if they are not monotonic
        "FritschButland"        monotonic - faster algorithm () but somewhat higher apparent "tension"
        "Steffen"               monotonic - also only one pass, results usualonly requires one passly between FritschCarlson and FritschButland
    Sources:
        Fritsch & Carlson (1980), "Monotone Piecewise Cubic Interpolation", doi:10.1137/0717021.
        Fritsch & Butland (1984), "A Method for Constructing Local Monotone Piecewise Cubic Interpolants", doi:10.1137/0905021.
        Steffen (1990), "A Simple Method for Monotonic Interpolation in One Dimension", http://adsabs.harvard.edu/abs/1990A%26A...239..443S
    
    Year : (c) 2023 
    Author : Rene De Ren
    Contact details : zn375ix3@gmail.com 
    Location : The Netherlands

*/

class Interpolation {
  constructor(config = {}) {
    this.input_xdata = [];
    this.input_ydata = [];
    this.y2 = [];
    this.n = 0;
    this.error = 0;
    this.interpolationtype = config.type || "monotone_cubic_spline";
    this.tension = config.tension || 0.5;
  }

  load_spline(input_xdata, input_ydata, interpolationtype) {
    if (!Array.isArray(input_xdata) || !Array.isArray(input_ydata)) {
      throw new Error("Invalid input: x and y must be arrays");
    }

    if (input_xdata.length !== input_ydata.length) {
      throw new Error("Arrays x and y must have the same length");
    }

    if (input_xdata.length < 2) {
      throw new Error("Arrays must contain at least 2 points for interpolation");
    }

    for (let i = 1; i < input_xdata.length; i++) {
      if (input_xdata[i] <= input_xdata[i - 1]) {
        throw new Error("X values must be strictly increasing");
      }
    }

    this.input_xdata = this.array_values(input_xdata);
    this.input_ydata = this.array_values(input_ydata);
    this.set_type(interpolationtype);
  }

  array_values(obj) {
    const new_array = [];
    for (let i in obj) {
      if (obj.hasOwnProperty(i)) {
        new_array.push(obj[i]);
      }
    }
    return new_array;
  }

  set_type(type) {
    if (type == "cubic_spline") {
      this.cubic_spline();
    } else if (type == "monotone_cubic_spline") {
      this.monotonic_cubic_spline();
    } else if (type == "linear") {
    } else {
      this.error = 1000;
    }
    this.interpolationtype = type;
  }

  interpolate(xpoint) {
    if (!this.input_xdata || !this.input_ydata || this.input_xdata.length < 2) {
      throw new Error("Spline not properly initialized");
    }

    if (xpoint <= this.input_xdata[0]) return this.input_ydata[0];
    if (xpoint >= this.input_xdata[this.input_xdata.length - 1]) return this.input_ydata[this.input_ydata.length - 1];

    let interpolatedval = 0;

    if (this.interpolationtype == "cubic_spline") {
      interpolatedval = this.interpolate_cubic(xpoint);
    } else if (this.interpolationtype == "monotone_cubic_spline") {
      interpolatedval = this.interpolate_cubic_monotonic(xpoint);
    } else if (this.interpolationtype == "linear") {
      interpolatedval = this.linear(xpoint);
    } else {
      console.log(this.interpolationtype);
      interpolatedval = "Unknown type";
    }
    return interpolatedval;
  }

  cubic_spline() {
    var xdata = this.input_xdata;
    var ydata = this.input_ydata;
    var delta = [];

    var n = ydata.length;
    this.n = n;

    if (n !== xdata.length) {
      this.error = 1;
    }

    this.y2[0] = 0.0;
    this.y2[n - 1] = 0.0;
    delta[0] = 0.0;

    for (let i = 1; i < n - 1; ++i) {
      let d = xdata[i + 1] - xdata[i - 1];

      if (d == 0) {
        this.error = 2;
      }

      let s = (xdata[i] - xdata[i - 1]) / d;

      let p = s * this.y2[i - 1] + 2.0;

      this.y2[i] = (s - 1.0) / p;

      delta[i] = (ydata[i + 1] - ydata[i]) / (xdata[i + 1] - xdata[i]) - (ydata[i] - ydata[i - 1]) / (xdata[i] - xdata[i - 1]);

      delta[i] = (6.0 * delta[i]) / (xdata[i + 1] - xdata[i - 1]) - (s * delta[i - 1]) / p;
    }

    for (let j = n - 2; j >= 0; --j) {
      this.y2[j] = this.y2[j] * this.y2[j + 1] + delta[j];
    }
  }

  linear(xpoint) {
    var i_min = 0;
    var i_max = 0;
    var o_min = 0;
    var o_max = 0;

    for (let i = 0; i < this.input_xdata.length; i++) {
      if (xpoint >= this.input_xdata[i] && xpoint < this.input_xdata[i + 1]) {
        i_min = this.input_xdata[i];
        i_max = this.input_xdata[i + 1];

        o_min = this.input_ydata[i];
        o_max = this.input_ydata[i + 1];

        break;
      }
    }

    let o_number;

    if (i_min < i_max) {
      o_number = o_min + ((xpoint - i_min) * (o_max - o_min)) / (i_max - i_min);
    } else {
      o_number = xpoint;
    }

    return o_number;
  }

  interpolate_cubic(xpoint) {
    let xdata = this.input_xdata;
    let ydata = this.input_ydata;

    let max = this.n - 1;
    let min = 0;

    while (max - min > 1) {
      let k = Math.floor((max + min) / 2);

      if (xdata[k] > xpoint) max = k;
      else min = k;
    }

    let h = xdata[max] - xdata[min];

    if (h == 0) {
      this.error = 3;
    }

    let a = (xdata[max] - xpoint) / h;
    let b = (xpoint - xdata[min]) / h;

    let interpolatedvalue = a * ydata[min] + b * ydata[max] + ((a * a * a - a) * this.y2[min] + (b * b * b - b) * this.y2[max]) * (h * h) / 6.0;

    return interpolatedvalue;
  }

  monotonic_cubic_spline() {
    let xdata = this.input_xdata;
    let ydata = this.input_ydata;

    let interpolationtype = this.interpolationtype;
    let tension = this.tension;

    let n = ydata.length;
    this.n = n;

    if (this.n !== xdata.length) {
      this.error = 1;
    }

    let obj = this.calc_tangents(xdata, ydata, tension);
    this.y1 = obj[0];
    this.delta = obj[1];
  }

  interpolate_cubic_monotonic(xpoint) {
    let xdata = this.input_xdata;
    let ydata = this.input_ydata;
    let xinterval = 0;

    let y1 = this.y1;
    let delta = this.delta;
    let c = [];
    let d = [];
    let n = this.n;

    for (let k = 0; k < n - 1; k++) {
      xinterval = xdata[k + 1] - xdata[k];
      c[k] = (3 * delta[k] - 2 * y1[k] - y1[k + 1]) / xinterval;
      d[k] = (y1[k] + y1[k + 1] - 2 * delta[k]) / xinterval / xinterval;
    }

    let interpolatedvalues = [];
    let k = 0;

    if (xpoint < xdata[0] || xpoint > xdata[n - 1]) {
    }

    while (k < n - 1 && xpoint > xdata[k + 1] && !(xpoint < xdata[0] || xpoint > xdata[n - 1])) {
      k++;
    }

    let xdiffdown = xpoint - xdata[k];

    interpolatedvalues = ydata[k] + y1[k] * xdiffdown + c[k] * xdiffdown * xdiffdown + d[k] * xdiffdown * xdiffdown * xdiffdown;

    return interpolatedvalues;
  }

  calc_tangents(xdata, ydata, tension) {
    let method = this.interpolationtype;
    let n = xdata.length;
    let delta_array = [];
    let delta = 0;
    let y1 = [];

    for (let i = 0; i < n - 1; i++) {
      delta = (ydata[i + 1] - ydata[i]) / (xdata[i + 1] - xdata[i]);
      delta_array[i] = delta;

      if (i == 0) {
        y1[i] = delta;
      } else if (method == "cardinal") {
        y1[i] = (1 - tension) * (ydata[i + 1] - ydata[i - 1]) / (xdata[i + 1] - xdata[i - 1]);
      } else if (method == "fritschbutland") {
        let alpha = (1 + (xdata[i + 1] - xdata[i]) / (xdata[i + 1] - xdata[i - 1])) / 3;
        y1[i] = delta_array[i - 1] * delta <= 0 ? 0 : (delta_array[i - 1] * delta) / (alpha * delta + (1 - alpha) * delta_array[i - 1]);
      } else if (method == "fritschcarlson") {
        y1[i] = delta_array[i - 1] * delta < 0 ? 0 : (delta_array[i - 1] + delta) / 2;
      } else if (method == "steffen") {
        let p = ((xdata[i + 1] - xdata[i]) * delta_array[i - 1] + (xdata[i] - xdata[i - 1]) * delta) / (xdata[i + 1] - xdata[i - 1]);
        y1[i] = (Math.sign(delta_array[i - 1]) + Math.sign(delta)) * Math.min(Math.abs(delta_array[i - 1]), Math.abs(delta), 0.5 * Math.abs(p));
      } else {
        y1[i] = (delta_array[i - 1] + delta) / 2;
      }
    }

    y1[n - 1] = delta_array[n - 2];
    if (method != "fritschcarlson") {
      return [y1, delta_array];
    }

    for (let i = 0; i < n - 1; i++) {
      let delta = delta_array[i];
      if (delta == 0) {
        y1[i] = 0;
        y1[i + 1] = 0;
        continue;
      }
      let alpha = y1[i] / delta;
      let beta = y1[i + 1] / delta;
      let tau = 3 / Math.sqrt(Math.pow(alpha, 2) + Math.pow(beta, 2));
      if (tau < 1) {
        y1[i] = tau * alpha * delta;
        y1[i + 1] = tau * beta * delta;
      }
    }
    return [y1, delta_array];
  }

  interpolate_lin_curve_points(i_curve, o_min, o_max) {
    if (!Array.isArray(i_curve)) {
      throw new Error("xArray must be an array");
    }

    let o_curve = {};
    let i_min = 0;
    let i_max = 0;

    i_min = Math.min(...Object.values(i_curve));
    i_max = Math.max(...Object.values(i_curve));

    i_curve.forEach((val, index) => {
      o_curve[index] = this.interpolate_lin_single_point(val, i_min, i_max, o_min, o_max);
    });

    o_curve = Object.values(o_curve);

    return o_curve;
  }

  interpolate_lin_single_point(i_number, i_min, i_max, o_min, o_max) {
    if (typeof i_number !== "number" || typeof i_min !== "number" || typeof i_max !== "number" || typeof o_min !== "number" || typeof o_max !== "number") {
      throw new Error("All parameters must be numbers");
    }

    if (i_max === i_min) {
      return o_min;
    }

    let o_number;
    //i_number = this.limit_input(i_number, i_min, i_max);

    o_number = o_min + ((i_number - i_min) * (o_max - o_min)) / (i_max - i_min);

    o_number = this.limit_input(o_number, o_min, o_max);

    return o_number;
  }

  limit_input(input, min, max) {
    let output;

    if (input < min) {
      output = min;
    } else if (input > max) {
      output = max;
    } else {
      output = input;
    }
    return output;
  }
}

module.exports = Interpolation;