testInterpolation.cxx

Go to the documentation of this file.

//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 "TList.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, bool drawSame = false ) {

   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++;
   }

   if (showGraphics) {
      TGraph *gr = new TGraph(n,xcoord,ycoord);
      gr->SetMarkerColor(kBlue);
      if (drawSame)     gr->SetMarkerColor(kGreen);
      gr->SetMarkerStyle(7);
      gr->Draw("CP");
      if (!drawSame) {
         TLegend * l = new TLegend(0.1,0.7,0.4,0.9);
         l->SetName("legend");
         l->AddEntry(grorig,"original data");
         l->AddEntry(gr,type.c_str());
         l->Draw();
      }
      else if (gPad) {
         TLegend * l = (TLegend*)  gPad->GetListOfPrimitives()->FindObject("legend");
         if (l) {
            l->AddEntry(gr, type.c_str() );
            l->Draw();
         }
      }
   }
   if (gPad) gPad->Update();

   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 = 0;
   if (showGraphics) {
      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);


   if (showGraphics) {
      c1->cd(2);
      grorig->Draw("AP");
   }

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


   if (showGraphics) {
      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);

   if (showGraphics) {
      c1->cd(4);
      grorig->Draw("AP");
   }

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

   if (showGraphics) {
      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);

   if (showGraphics) {
      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,true);


}

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

  // Parse command line arguments
  for (Int_t i=1 ;  i<argc ; i++) {
     std::string arg = argv[i] ;
     if (arg == "-g") {
      showGraphics = true;
     }
     if (arg == "-v") {
      showGraphics = true;
      //verbose = true;
     }
     if (arg == "-h") {
        cerr << "Usage: " << argv[0] << " [-g] [-v]\n";
        cerr << "  where:\n";
        cerr << "     -g : graphics mode\n";
        cerr << "     -v : verbose  mode";
        cerr << endl;
        return -1;
     }
   }

   TApplication* theApp = 0;

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

   testInterpolation();

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

   return 0;
}
#endif