//Example of the various methods of interpolation provided by the 
// ROOT::Math::Interpolator class
// 
//Example can also be run in ROOT by doing : 
// 
// root> .x testInterpolation.cxx
//
#include "TGraph.h"
#include "TAxis.h"
#include "TCanvas.h"
#include "TLegend.h"
#include "TApplication.h"
#include "Math/Interpolator.h"
#include <iostream>
#include <string>
#include <cstdlib>

#include <cmath>

bool showGraphics = true;

TGraph *grorig = 0;

void interpolate( const  ROOT::Math::Interpolator & itp ) { 

   std::string type = itp.Type(); 
   std::cout << "\n" << type << "  interpolation:" << std::endl;

   std::cout << "x[i]     y[i]     deriv[i]     deriv2[i]    integral[i] \n" << std::endl; 
   // print result of interpolation
   const Int_t n = 51;
   Int_t i = 0;
   Float_t xcoord[n], ycoord[n];
   for ( double xi = 0; xi < 10; xi += 0.2) { 
      xcoord[i] = xi;
      ycoord[i] = itp.Eval(xi);
      double dyi = itp.Deriv(xi); 
      double dy2i = itp.Deriv2(xi);
      double igyi = itp.Integ(0, xi); 
      std::cout << xcoord[i]  << "  " << ycoord[i] << "  " << dyi << "  " << dy2i << "  " <<  igyi << std::endl;
      i++; 
   }

   TGraph *gr = new TGraph(n,xcoord,ycoord);
   gr->SetMarkerColor(kBlue);
   gr->SetMarkerStyle(7);
   gr->Draw("CP");
   TLegend * l = new TLegend(0.1,0.7,0.4,0.9); 
   l->AddEntry(grorig,"original data"); 
   l->AddEntry(gr,type.c_str());
   l->Draw();

   return;
}


void testInterpolation() {

   // Create data
   int n = 10; 
   Float_t xorig[11], yorig[11];
   std::vector<double> x(n+1); 
   std::vector<double> y(n+1); 
   for (int i = 0; i < n+1; ++i) {  
      x[i] = i + 0.5 * std::sin(i+0.0); 
      y[i] = i + std::cos(i * i+0.0); 
      xorig[i] = x[i];
      yorig[i] = y[i];
   } 

   TCanvas *c1 = new TCanvas("c1","Original (red), Linear (upper left), Polynomial (upper right), Spline , Spline periodic, Akima (lower left) and Akima Periodic (lower right) Interpolation",10,10,1000,800);
   c1->Divide(2,3);
   c1->cd(1);

   grorig = new TGraph(n+1,xorig,yorig);
   grorig->SetMarkerColor(kRed);
   grorig->SetMarkerStyle(20);
   grorig->GetYaxis()->SetRange(0,40);
   grorig->Draw("AP"); 

   //ROOT::Math::Interpolator itp1(x, y, ROOT::Math::Interpolation::kLINEAR); 
   ROOT::Math::Interpolator itp1(x.size(), ROOT::Math::Interpolation::kLINEAR); 
   itp1.SetData(x,y);
   interpolate(itp1);


   c1->cd(2);
   grorig->Draw("AP"); 

   ROOT::Math::Interpolator itp2(x, y, ROOT::Math::Interpolation::kPOLYNOMIAL); 
   interpolate(itp2);


   c1->cd(3);
   grorig->Draw("AP"); 

   //std::cout << "Cubic Spline Interpolation: " << std::endl;
   ROOT::Math::Interpolator itp3( 2*x.size(), ROOT::Math::Interpolation::kCSPLINE); 
   itp3.SetData(x.size(), &x[0], &y[0]);
   interpolate(itp3);

   c1->cd(4);
   grorig->Draw("AP"); 

   //std::cout << "Akima  Interpolation: " << std::endl;
   ROOT::Math::Interpolator itp4(x, y, ROOT::Math::Interpolation::kAKIMA); 
   interpolate(itp4);

   c1->cd(5);
   grorig->Draw("AP"); 

   //std::cout << "Cubic Spline Periodic Interpolation: " << std::endl;
   ROOT::Math::Interpolator itp5(x, y, ROOT::Math::Interpolation::kCSPLINE_PERIODIC); 
   interpolate(itp5);



   c1->cd(6);
   grorig->Draw("AP"); 

   //std::cout << "Akima Periodic Interpolation: " << std::endl;
   ROOT::Math::Interpolator itp6(x, y, ROOT::Math::Interpolation::kAKIMA_PERIODIC); 
   interpolate(itp6);

   std::cout << "\n***********************************" << std::endl;
   std::cout << "Using default Interpolation type: " << std::endl;

   ROOT::Math::Interpolator itp7; 
   itp7.SetData(x,y);
   interpolate(itp7);


}

#ifndef __CINT__
int main(int argc, char **argv)
{
   using std::cerr;
   using std::cout;
   using std::endl;

   if ( argc > 1 && argc != 2 )
   {
      cerr << "Usage: " << argv[0] << " [-ng]\n";
      cerr << "  where:\n";
      cerr << "     -ng : no graphics mode";
      cerr << endl;
      exit(1);
   }

   if ( argc == 2 && strcmp( argv[1], "-ng") == 0 ) 
   {
      showGraphics = false;
   }

   TApplication* theApp = 0;

   if ( showGraphics )
      theApp = new TApplication("App",&argc,argv);

   testInterpolation();

   if ( showGraphics )
   {
      theApp->Run();
      delete theApp;
      theApp = 0;
   }

   return 0;
}
#endif