#ifndef SPLINESJKPGE
#define SPLINESJKPGE

#include <string>
#include <iostream>
#include <math.h>

using namespace std;
#include "TObject.h"
#include "TH1.h"
#include "TF1.h"

#include "Table.hh"
#include "foreach.hh"

struct SplineKnot
{
  double x, f, df;
  double dx;
  double a, b, c, d; // coefficients for spline to the right of knot

  SplineKnot(double x_, double f_, double df_) : x(x_), f(f_), df(df_) {}
  SplineKnot() {}

  void set_coeffs(double x1, double x2, double f1, double f2, double df1, double df2)
  {
    dx = x2 - x1;
    df1 *= dx;
    df2 *= dx;

    const double A = df2 - df1;
    const double B = f2 - f1 - df1;
    a = f1;
    b = df1;
    d = A - 2 * B;
    c = B - d;
  }

  double df1() { return eval_df(x); }
  double df2() { return eval_df(x + dx); }

  double eval(double X) const
  {
    const double v = (X - x) / dx;
    return a + b * v + c * v * v + d * v * v * v;
  }

  double eval_df(double X) const
  {
    const double v = (X - x) / dx;
    return (b + 2 * c * v + 3 * d * v * v) / dx;
  }

  void eval_derivative(double X, double &f, double &df)
  {
    const double v = (X - x) / dx;
    const double v2 = v * v;
    f = a + b * v + c * v2 + d * v2 * v;
    df = (b + 2 * c * v + 3 * d * v2) / dx;
  }
};

inline double sign(double x) { return x / fabs(x); }

inline double interpolate_monotone2(double x1, double x2, double x3, double x4,
                                    double f1, double f2, double f3, double f4,
                                    double x, bool mono = true)
{
  //--------------------------------------------------
  // M. Steffen, A&A Vol 239, P443, 1990
  //--------------------------------------------------

  const double h1 = x2 - x1;
  const double h2 = x3 - x2;
  const double h3 = x4 - x3;

  const double s1 = (f2 - f1) / h1;
  const double s2 = (f3 - f2) / h2;
  const double s3 = (f4 - f3) / h3;

  double p2 = (s1 * h2 + s2 * h1) / (h1 + h2);
  double p3 = (s2 * h3 + s3 * h2) / (h2 + h3);

  if (mono)
  {
    double smin2 = min(fabs(s1), fabs(s2));
    if (fabs(p2) > 2 * smin2)
      p2 = 2 * smin2 * sign(p2);
    double smin3 = min(fabs(s2), fabs(s3));
    if (fabs(p3) > 2 * smin3)
      p3 = 2 * smin3 * sign(p3);
  }

  // evaluate the spline
  double dx = x3 - x2;
  const double df2 = p2 * dx;
  const double df3 = p3 * dx;
  const double A = df3 - df2;
  const double B = f3 - f2 - df2;
  double a = f2;
  double b = df2;
  double d = A - 2 * B;
  double c = B - d;
  const double v = (x - x2) / dx;
  return a + b * v + c * v * v + d * v * v * v;
}

inline double interpolate_monotone2(double xvals[],
                                    double fvals[],
                                    double x,
                                    bool mono = true)
{
  return interpolate_monotone2(xvals[0], xvals[1], xvals[2], xvals[3],
                               fvals[0], fvals[1], fvals[2], fvals[3],
                               x, mono);
}

inline double interpolate_monotone(double x1, double x2, double x3, double x4,
                                   double f1, double f2, double f3, double f4,
                                   double x, bool mono = true)
{

  const double d1 = (f2 - f1) / (x2 - x1);
  const double d2 = (f3 - f2) / (x3 - x2);
  const double d3 = (f4 - f3) / (x4 - x3);

  double m1 = d1;
  double m2 = (d2 + d1) / 2;
  double m3 = (d3 + d2) / 2;
  double m4 = d3;

  if (d1 == 0)
    m1 = m2 = 0;
  if (d2 == 0)
    m2 = m3 = 0;
  if (d3 == 0)
    m4 = m3 = 0;

  (void)m1;

  double a2 = m2 / d2;
  double b2 = m3 / d2;
  double a3 = m3 / d3;
  double b3 = m4 / d3;

  double alpha = a2;
  double beta = b2;
  double r = sqrt(alpha * alpha + beta * beta);
  double t = r > 3 ? 3 / r : 1;
  if (!mono)
    t = 1;
  double df2 = t * alpha * d2;
  double df3 = t * beta * d2;

  alpha = a3;
  beta = b3;
  double r2 = sqrt(alpha * alpha + beta * beta);
  t = r2 > 3 ? 3 / r2 : 1;
  if (!mono)
    t = 1;
  if (mono && r2 > r)
    df3 = t * alpha * d3;

  print(r, r2, df2, df3);

  // print ("d123, df2 df3 = ",d1, d2,d3,  df2, df3);

  // evaluate the spline
  double dx = x3 - x2;
  df2 *= dx;
  df3 *= dx;
  const double A = df3 - df2;
  const double B = f3 - f2 - df2;
  double a = f2;
  double b = df2;
  double d = A - 2 * B;
  double c = B - d;

  const double v = (x - x2) / dx;
  return a + b * v + c * v * v + d * v * v * v;
}

struct Spline : TObject
{
  vector<SplineKnot> knots;

  double xmin() { return knots[0].x; }
  double xmax() { return knots[knots.size() - 1].x; }
  Spline copy() { return *this; }

  Spline() {}
  Spline(vector<double> x, vector<double> y, string option)
  {
    for (int i = 0; i < (int)x.size(); i++)
      add_point(x[i], y[i]);
    init(option);
  }

  Spline(TH1 &h, string option = "mono")
  {
    int N = h.GetNbinsX();
    add_point(h.GetBinCenter(0), h.GetBinContent(1));

    for (int b = 1; b <= N; b++)
    {
      add_point(h.GetBinCenter(b), h.GetBinContent(b));
    }

    add_point(h.GetBinCenter(N + 1), h.GetBinContent(N));
    init(option);
  }

  void add_point(double x, double f, double df = 0)
  {
    knots.push_back(SplineKnot(x, f, df));
  }

  void init(string opt = "mono")
  {
    init_derivatives();
    if (opt == "mono")
    {
      // cout << "--------------------make monotone " << endl;
      make_monotone();
    }
    set_coeffs();
  }

  void init_derivatives()
  {
    for (int i = 1; i < (int)knots.size() - 1; i++)
    {
      const double d1 = (knots[i + 1].f - knots[i].f) / (knots[i + 1].x - knots[i].x);
      const double d2 = (knots[i].f - knots[i - 1].f) / (knots[i].x - knots[i - 1].x);

      if (d1 * d2 < 0)
      {
        knots[i].df = 0;
      }
      else
      {
        knots[i].df = (d1 + d2) / 2;
      }

      if (d1 == 0)
        knots[i].df = 0;
      if (d2 == 0)
        knots[i].df = knots[i - 1].df = 0;
    }
  }

  bool make_monotone()
  {
    for (int i = 0; i < (int)knots.size() - 1; i++)
    {
      SplineKnot &k1 = knots[i];
      SplineKnot &k2 = knots[i + 1];

      const double delta_k = (k2.f - k1.f) / (k2.x - k1.x);

      if (delta_k == 0)
        continue;

      const double alpha = k1.df / delta_k;
      const double beta = k2.df / delta_k;

      // print ("knot ",i,k1.f,k1.df , alpha, delta_k, k1.f, k1.df);

      if (alpha < 0)
      {
        cout << " data not monotone in point " << i << endl;

        // for (int j=0; j< (int)knots.size()-1; j++)
        //  {
        // 	SplineKnot& k = knots[j];
        //  }
        fatal("");
        continue;
      }

      const double r = sqrt(alpha * alpha + beta * beta);
      const double t = r > 3 ? 3 / r : 1;
      k1.df = t * alpha * delta_k;
      k2.df = t * beta * delta_k;
    }

    set_coeffs();
    return true;
  }

  void set_coeffs()
  {
    for (int i = 0; i < (int)knots.size() - 1; i++)
    {
      knots[i].set_coeffs(knots[i].x, knots[i + 1].x,
                          knots[i].f, knots[i + 1].f,
                          knots[i].df, knots[i + 1].df);
    }
  }

  void dump()
  {
    for (int i = 0; i < (int)knots.size(); i++)
    {
      cout << i << "   x=" << knots[i].x << ", " << knots[i].f << " " << knots[i].df << endl;
    }
  }

  double eval(double x) const
  {
    if ((int)knots.size() == 0)
    {
      cout << "warning: trying to evaluate empty spline " << endl;
      return -1e100;
    }

    if (x < knots[0].x || x > knots[knots.size() - 1].x)
      return 0;

    int i = 0;
    while (x > knots[i].x)
      i++;
    return knots[i - 1].eval(x);
  }

  void eval_derivative(double x, double &f, double &df)
  {
    if ((int)knots.size() == 0)
    {
      cout << "warning: trying to evaluate empty spline " << endl;
    }

    if (x < knots[0].x || x > knots[knots.size() - 1].x)
    {
      x = df = 0;
    };

    int i = 0;
    while (x > knots[i].x)
      i++;
    knots[i - 1].eval_derivative(x, f, df);
  }

  double operator()(double x) const
  {
    return eval(x);
  }

  TH1F as_hist(int nbins = 50, string name = "")
  {
    TH1F r(("splh" + name).c_str(), "splh", nbins, xmin(), xmax());
    for (int i = 1; i <= r.GetNbinsX(); i++)
    {
      double x = r.GetBinCenter(i);
      r.SetBinContent(i, eval(x));
    }
    return r;
  }

  double operator()(double *x, double *dum) const
  {
    return eval(*x);
  }

  TF1 as_tf1()
  {
    return TF1("spline", *this, xmin(), xmax(), 0);
  }

  void Draw(string opt = "")
  {
    as_hist().Draw(opt.c_str());
  }

  string info();

  ClassDef(Spline, 2)
};

#endif