// @(#)root/hist:$Id$
// Author: Federico Carminati   28/02/2000

/*************************************************************************
 * Copyright (C) 1995-2022, Rene Brun and Fons Rademakers.               *
 * All rights reserved.                                                  *
 *                                                                       *
 * For the licensing terms see $ROOTSYS/LICENSE.                         *
 * For the list of contributors see $ROOTSYS/README/CREDITS.             *
 *************************************************************************/

#ifndef ROOT_TSpline
#define ROOT_TSpline

#ifdef R__LESS_INCLUDES
#include "TNamed.h"
#include "TAttLine.h"
#include "TAttFill.h"
#include "TAttMarker.h"
#else
#include "TGraph.h"
#endif

class TH1;
class TH1F;
class TF1;
class TGraph;

class TSpline : public TNamed, public TAttLine,
                public TAttFill, public TAttMarker
{
protected:
   Double_t fDelta = -1.;      ///< Distance between equidistant knots
   Double_t fXmin = 0.;        ///< Minimum value of abscissa
   Double_t fXmax = 0.;        ///< Maximum value of abscissa
   Int_t fNp = 0;              ///< Number of knots
   Bool_t fKstep = kFALSE;     ///< True of equidistant knots
   TH1F *fHistogram = nullptr; ///< Temporary histogram
   TGraph *fGraph = nullptr;   ///< Graph for drawing the knots
   Int_t fNpx = 100;           ///< Number of points used for graphical representation

   TSpline(const TSpline &);
   TSpline &operator=(const TSpline &);
   virtual void BuildCoeff() = 0;

public:
   TSpline() {}
   TSpline(const char *title, Double_t delta, Double_t xmin, Double_t xmax, Int_t np, Bool_t step)
      : TNamed("Spline", title), TAttFill(0, 1), fDelta(delta), fXmin(xmin), fXmax(xmax), fNp(np), fKstep(step)
   {
   }
   ~TSpline() override;

   virtual void     GetKnot(Int_t i, Double_t &x, Double_t &y) const =0;
   Int_t    DistancetoPrimitive(Int_t px, Int_t py) override;
   void     Draw(Option_t *option="") override;
   void     ExecuteEvent(Int_t event, Int_t px, Int_t py) override;
   virtual Double_t GetDelta() const {return fDelta;}
   TH1F            *GetHistogram() const {return fHistogram;}
   virtual Int_t    GetNp()    const {return fNp;}
   virtual Int_t    GetNpx()   const {return fNpx;}
   virtual Double_t GetXmin()  const {return fXmin;}
   virtual Double_t GetXmax()  const {return fXmax;}
   void     Paint(Option_t *option="") override;
   virtual Double_t Eval(Double_t x) const=0;
   void     SaveAs(const char * /*filename*/ = "",Option_t * /*option*/ = "") const override {}
   void             SetNpx(Int_t n) {fNpx=n;}

   ClassDefOverride(TSpline,2) // Spline base class
};


//______________________________________________________________________________
class TSplinePoly : public TObject
{
protected:
   Double_t fX = 0.;     ///< Abscissa
   Double_t fY = 0.;     ///< Constant term

public:
   TSplinePoly() {}
   TSplinePoly(Double_t x, Double_t y) : fX(x), fY(y) {}
   TSplinePoly(TSplinePoly const &other) : TObject(other), fX(0), fY(0) { CopyPoly(other); }

   TSplinePoly &operator=(TSplinePoly const &other);

   Double_t &X() { return fX; }
   Double_t &Y() { return fY; }
   void GetKnot(Double_t &x, Double_t &y) const
   {
      x = fX;
      y = fY;
   }

   virtual Double_t Eval(Double_t) const { return fY; }

private:
   void CopyPoly(TSplinePoly const &other);

   ClassDefOverride(TSplinePoly,2) // Spline polynomial terms
};

//______________________________________________________________________________
class TSplinePoly3 : public TSplinePoly
{
private:
   Double_t fB = 0.; ///< First order expansion coefficient :  fB*1! is the first derivative at x
   Double_t fC = 0.; ///< Second order expansion coefficient : fC*2! is the second derivative at x
   Double_t fD = 0.; ///< Third order expansion coefficient :  fD*3! is the third derivative at x

public:
   TSplinePoly3() {}
   TSplinePoly3(Double_t x, Double_t y, Double_t b, Double_t c, Double_t d) : TSplinePoly(x, y), fB(b), fC(c), fD(d) {}
   TSplinePoly3(TSplinePoly3 const &other) : TSplinePoly(other) { CopyPoly(other); }

   TSplinePoly3 &operator=(TSplinePoly3 const &other);

   Double_t &B() { return fB; }
   Double_t &C() { return fC; }
   Double_t &D() { return fD; }
   Double_t Eval(Double_t x) const override
   {
      Double_t dx = x - fX;
      return (fY + dx * (fB + dx * (fC + dx * fD)));
   }
   Double_t Derivative(Double_t x) const
   {
      Double_t dx = x - fX;
      return (fB + dx * (2 * fC + 3 * fD * dx));
   }

private:
   void CopyPoly(TSplinePoly3 const &other);

   ClassDefOverride(TSplinePoly3,1)  // Third spline polynomial terms
};

//______________________________________________________________________________
class TSplinePoly5 : public TSplinePoly
{
private:
   Double_t fB = 0.; ///< First order expansion coefficient :  fB*1! is the first derivative at x
   Double_t fC = 0.; ///< Second order expansion coefficient : fC*2! is the second derivative at x
   Double_t fD = 0.; ///< Third order expansion coefficient :  fD*3! is the third derivative at x
   Double_t fE = 0.; ///< Fourth order expansion coefficient : fE*4! is the fourth derivative at x
   Double_t fF = 0.; ///< Fifth order expansion coefficient :  fF*5! is the fifth derivative at x

public:
   TSplinePoly5() {}
   TSplinePoly5(Double_t x, Double_t y, Double_t b, Double_t c, Double_t d, Double_t e, Double_t f)
      : TSplinePoly(x, y), fB(b), fC(c), fD(d), fE(e), fF(f)
   {
   }
   TSplinePoly5(TSplinePoly5 const &other) : TSplinePoly(other) { CopyPoly(other); }

   TSplinePoly5 &operator=(TSplinePoly5 const &other);

   Double_t &B() { return fB; }
   Double_t &C() { return fC; }
   Double_t &D() { return fD; }
   Double_t &E() { return fE; }
   Double_t &F() { return fF; }
   Double_t Eval(Double_t x) const override
   {
      Double_t dx=x-fX;
      return (fY+dx*(fB+dx*(fC+dx*(fD+dx*(fE+dx*fF)))));
   }
   Double_t Derivative(Double_t x) const{
      Double_t dx=x-fX;
      return (fB+dx*(2*fC+dx*(3*fD+dx*(4*fE+dx*(5*fF)))));
   }

private:
   void CopyPoly(TSplinePoly5 const &other);

   ClassDefOverride(TSplinePoly5,1)  // Quintic spline polynomial terms
};


//______________________________________________________________________________
class TSpline3 : public TSpline
{
protected:
   TSplinePoly3 *fPoly = nullptr; ///<[fNp] Array of polynomial terms
   Double_t fValBeg = 0.;         ///< Initial value of first or second derivative
   Double_t fValEnd = 0.;         ///< End value of first or second derivative
   Int_t fBegCond = -1;           ///< 0=no beg cond, 1=first derivative, 2=second derivative
   Int_t fEndCond = -1;           ///< 0=no end cond, 1=first derivative, 2=second derivative

   void   BuildCoeff() override;
   void   SetCond(const char *opt);

public:
   TSpline3() {}
   TSpline3(const char *title,
            Double_t x[], Double_t y[], Int_t n, const char *opt=nullptr,
            Double_t valbeg=0, Double_t valend=0);
   TSpline3(const char *title,
            Double_t xmin, Double_t xmax,
            Double_t y[], Int_t n, const char *opt=nullptr,
            Double_t valbeg=0, Double_t valend=0);
   TSpline3(const char *title,
            Double_t x[], const TF1 *func, Int_t n, const char *opt=nullptr,
            Double_t valbeg=0, Double_t valend=0);
   TSpline3(const char *title,
            Double_t xmin, Double_t xmax,
            const TF1 *func, Int_t n, const char *opt=nullptr,
            Double_t valbeg=0, Double_t valend=0);
   TSpline3(const char *title,
            const TGraph *g, const char *opt=nullptr,
            Double_t valbeg=0, Double_t valend=0);
   TSpline3(const TH1 *h, const char *opt=nullptr,
            Double_t valbeg=0, Double_t valend=0);
   TSpline3(const TSpline3&);
   TSpline3& operator=(const TSpline3&);
   Int_t    FindX(Double_t x) const;
   Double_t Eval(Double_t x) const override;
   Double_t Derivative(Double_t x) const;
   ~TSpline3() override {if (fPoly) delete [] fPoly;}
   void GetCoeff(Int_t i, Double_t &x, Double_t &y, Double_t &b,
                 Double_t &c, Double_t &d) const
      {x=fPoly[i].X();y=fPoly[i].Y();
      b=fPoly[i].B();c=fPoly[i].C();d=fPoly[i].D();}
   void GetKnot(Int_t i, Double_t &x, Double_t &y) const override
      {x=fPoly[i].X(); y=fPoly[i].Y();}
    void     SaveAs(const char *filename="",Option_t *option="") const override;
    void     SavePrimitive(std::ostream &out, Option_t *option = "") override;
   virtual  void     SetPoint(Int_t i, Double_t x, Double_t y);
   virtual  void     SetPointCoeff(Int_t i, Double_t b, Double_t c, Double_t d);
   static void Test();

   ClassDefOverride(TSpline3,2)  // Class to create third natural splines
};


//______________________________________________________________________________
class TSpline5 : public TSpline
{
protected:
   TSplinePoly5  *fPoly = nullptr;     ///<[fNp] Array of polynomial terms

   void BuildCoeff() override;
   void BoundaryConditions(const char *opt, Int_t &beg, Int_t &end,
                           const char *&cb1, const char *&ce1, const char *&cb2,
                           const char *&ce2);
   void SetBoundaries(Double_t b1, Double_t e1, Double_t b2, Double_t e2,
                      const char *cb1, const char *ce1, const char *cb2,
                      const char *ce2);
public:
   TSpline5() {}
   TSpline5(const char *title,
            Double_t x[], Double_t y[], Int_t n,
            const char *opt=nullptr, Double_t b1=0, Double_t e1=0,
            Double_t b2=0, Double_t e2=0);
   TSpline5(const char *title,
            Double_t xmin, Double_t xmax,
            Double_t y[], Int_t n,
            const char *opt=nullptr, Double_t b1=0, Double_t e1=0,
            Double_t b2=0, Double_t e2=0);
   TSpline5(const char *title,
            Double_t x[], const TF1 *func, Int_t n,
            const char *opt=nullptr, Double_t b1=0, Double_t e1=0,
            Double_t b2=0, Double_t e2=0);
   TSpline5(const char *title,
            Double_t xmin, Double_t xmax,
            const TF1 *func, Int_t n,
            const char *opt=nullptr, Double_t b1=0, Double_t e1=0,
            Double_t b2=0, Double_t e2=0);
   TSpline5(const char *title,
            const TGraph *g,
            const char *opt=nullptr, Double_t b1=0, Double_t e1=0,
            Double_t b2=0, Double_t e2=0);
   TSpline5(const TH1 *h,
            const char *opt=nullptr, Double_t b1=0, Double_t e1=0,
            Double_t b2=0, Double_t e2=0);
   TSpline5(const TSpline5&);
   TSpline5& operator=(const TSpline5&);
   Int_t    FindX(Double_t x) const;
   Double_t Eval(Double_t x) const override;
   Double_t Derivative(Double_t x) const;
   ~TSpline5() override {if (fPoly) delete [] fPoly;}
   void GetCoeff(Int_t i, Double_t &x, Double_t &y, Double_t &b,
                 Double_t &c, Double_t &d, Double_t &e, Double_t &f) const
      {x=fPoly[i].X();y=fPoly[i].Y();b=fPoly[i].B();
      c=fPoly[i].C();d=fPoly[i].D();
      e=fPoly[i].E();f=fPoly[i].F();}
   void GetKnot(Int_t i, Double_t &x, Double_t &y) const override
      {x=fPoly[i].X(); y=fPoly[i].Y();}
    void     SaveAs(const char *filename="",Option_t *option="") const override;
    void     SavePrimitive(std::ostream &out, Option_t *option = "") override;
   virtual  void     SetPoint(Int_t i, Double_t x, Double_t y);
   virtual  void     SetPointCoeff(Int_t i, Double_t b, Double_t c, Double_t d,
                                   Double_t e, Double_t f);
   static void Test();

   ClassDefOverride(TSpline5,2) // Class to create quintic natural splines
};

#endif