// test unuran using the string interface to generate numbers according to the normal distributions
// compare CPU performancecwith TRandom::Gaus and opitonally GSL (using MathMore ) and CLHEP for 
// generating normal distributed random numbers 
//
// run within ROOT (.x unuranSimple.cxx+) or pass any extra parameter in the command line to get  
// a graphics output 
//
#include "TStopwatch.h"
#include "TUnuran.h"

#include "TH1.h"
#include "TF1.h"

#include "TRandom.h"
#include "TSystem.h"
#include "TApplication.h"
#include "TCanvas.h"
//#include "TRint.h"
#include "TVirtualFitter.h"
#include "TFitter.h"
#include "Math/DistFunc.h"

#include <iostream> 

#ifdef HAVE_MATHMORE
#include "Math/Random.h"
#include "Math/GSLRndmEngines.h"
#endif



#ifdef HAVE_CLHEP
#include "CLHEP/Random/RandFlat.h"
#include "CLHEP/Random/RandGauss.h"
#include "CLHEP/Random/MTwistEngine.h"
#include "CLHEP/Random/JamesRandom.h"
#include "CLHEP/Random/RanluxEngine.h"
#include "CLHEP/Random/Ranlux64Engine.h"
#include "CLHEP/Random/RanecuEngine.h"
#include "CLHEP/Random/Hurd160Engine.h"
#include "CLHEP/Random/Hurd288Engine.h"
#include "CLHEP/Random/RanshiEngine.h"
#include "CLHEP/Random/DualRand.h"
#include "CLHEP/Random/TripleRand.h"
#endif


using std::cout; 
using std::endl; 

int unuranSimple( ) { 

   // simple test of unuran


   std::cout << "Test Generation of Gaussian Numbers \n\n";

   TUnuran unr; 
   if (! unr.Init( "normal()", "method=arou") ) {
      cout << "Error initializing unuran" << endl;
      return -1;
   }

    // default is 10**7 but one should use 10**8 for serious timing
   int n = 10000000; 

   TStopwatch w;
   double time; 
   w.Start(); 

   for (int i = 0; i < n; ++i) 
      unr.Sample(); 

   w.Stop(); 
   time = w.CpuTime()*1.E9/n;
   cout << "Time using Unuran method arou =\t " <<   time << "\tns/call" << endl;


   if (! unr.Init( "normal()", "method=tdr") ) {
      cout << "Error initializing unuran" << endl;
      return -1;
   }
   w.Start(); 
   for (int i = 0; i < n; ++i) 
      unr.Sample(); 

   w.Stop(); 
   time = w.CpuTime()*1.E9/n;
   cout << "Time using Unuran method tdr =\t " <<   time << "\tns/call" << endl;

   if (! unr.Init( "normal()", "method=hinv") ) {
      cout << "Error initializing unuran" << endl;
      return -1;
   }
   w.Start(); 
   for (int i = 0; i < n; ++i) 
      unr.Sample(); 

   w.Stop(); 
   time = w.CpuTime()*1.E9/n;
   cout << "Time using Unuran method hinv =\t " <<   time << "\tns/call" << endl;

   w.Start();
   for (int i = 0; i < n; ++i) 
      gRandom->Gaus(0,1); 
   w.Stop(); 
   time = w.CpuTime()*1.E9/n;
   cout << "Time using TRandom::Gaus  =\t " <<   time << "\tns/call" << endl;

   // using Rannor
   w.Start();
   double x1,x2; 
   for (int i = 0; i < n/2; ++i) 
      gRandom->Rannor(x1,x2); 
   w.Stop(); 
   time = w.CpuTime()*1.E9/n;
   cout << "Time using TRandom::Rannor  =\t " <<   time << "\tns/call" << endl;

#ifdef HAVE_MATHMORE
   // using GSL Ziggurat
   ROOT::Math::Random<ROOT::Math::GSLRngMT>     rgsl;
   w.Start();
   for (int i = 0; i < n; ++i) 
      rgsl.Gaus(0,1); 
   w.Stop(); 
   time = w.CpuTime()*1.E9/n;
   cout << "Time using GSL::Gaus  =\t\t " <<   time << "\tns/call" << endl;

   // using GSL BoxMuller method
   w.Start();
   for (int i = 0; i < n; ++i) 
      rgsl.GausBM(0,1); 
   w.Stop(); 
   time = w.CpuTime()*1.E9/n;
   cout << "Time using GSL::GausBM  =  \t " <<   time << "\tns/call" << endl;

   // using GSL Ratio method
   w.Start();
   for (int i = 0; i < n; ++i) 
      rgsl.GausR(0,1); 
   w.Stop(); 
   time = w.CpuTime()*1.E9/n;
   cout << "Time using GSL::GausR  =\t " <<   time << "\tns/call" << endl;
#endif


   // run unuran standard Normal methods 2 (Polar from Marsaglia)

   if (! unr.Init( "normal()", "method=cstd;variant=2") ) {
      cout << "Error initializing unuran" << endl;
      return -1;
   }
   w.Start(); 
   for (int i = 0; i < n; ++i) 
      unr.Sample(); 

   w.Stop(); 
   time = w.CpuTime()*1.E9/n;
   cout << "Time using Unuran GausPolarR =\t " <<   time << "\tns/call" << endl;

   // run unuran standard Normal method 3 (KM)

   if (! unr.Init( "normal()", "method=cstd;variant=3") ) {
      cout << "Error initializing unuran" << endl;
      return -1;
   }
   w.Start(); 
   for (int i = 0; i < n; ++i) 
      unr.Sample(); 

   w.Stop(); 
   time = w.CpuTime()*1.E9/n;
   cout << "Time using Unuran Gaus  K-R   =\t " <<   time << "\tns/call" << endl;
  
   if (! unr.Init( "normal()", "method=cstd;variant=6") ) {
      cout << "Error initializing unuran" << endl;
      return -1;
   }
   w.Start(); 
   for (int i = 0; i < n; ++i) 
      unr.Sample(); 

   w.Stop(); 
   time = w.CpuTime()*1.E9/n;
   cout << "Time using Unuran Gaus exp6  =\t " <<   time << "\tns/call" << endl;


//    w.Start();
//    for (int i = 0; i < n; ++i) 
//       rgsl.GausRatio(0,1); 
//    w.Stop(); 
//    time = w.CpuTime()*1.E9/n;
//    cout << "Time using GSL::GausRatio  =\t\t " <<   time << "\tns/call" << endl;

#ifdef HAVE_CLHEP
   MTwistEngine eng(111);
   RandGauss rg(eng);
   w.Start();
   for (int i = 0; i < n; ++i) 
      rg(0,1); 
   w.Stop(); 
   time = w.CpuTime()*1.E9/n;
   cout << "Time using CLHEP::Gaus  =\t " <<   time << "\tns/call" << endl;
#endif

   std::cout << "\nTest uniform generator\n" << std::endl;
   w.Start();
   for (int i = 0; i < n; ++i) 
      gRandom->Rndm(); 
   w.Stop(); 
   time = w.CpuTime()*1.E9/n;
   cout << "Time using TRandom::Rndm  =\t " <<   time << "\tns/call" << endl;

#ifdef HAVE_CLHEP
   RandFlat rf(eng);
   w.Start();
   for (int i = 0; i < n; ++i) 
      eng.flat(); // use directly the engine (faster!)
   w.Stop(); 
   time = w.CpuTime()*1.E9/n;
   cout << "Time using CLHEP::MT  =\t\t " <<   time << "\tns/call" << endl;

   { 
      RanecuEngine e; 
      RandFlat rf2(e);
      w.Start();
      for (int i = 0; i < n; ++i) 
         e.flat(); 
      w.Stop(); 
      time = w.CpuTime()*1.E9/n;
      cout << "Time using CLHEP::Ranecu  =\t " <<   time << "\tns/call" << endl;
   }
   { 
      Hurd160Engine e; 
      RandFlat rf2(e);
      w.Start();
      for (int i = 0; i < n; ++i) 
         e.flat();
      w.Stop(); 
      time = w.CpuTime()*1.E9/n;
      cout << "Time using CLHEP::Hard160  =\t " <<   time << "\tns/call" << endl;
   }
   { 
      Hurd288Engine e; 
      RandFlat rf2(e);
      w.Start();
      for (int i = 0; i < n; ++i) 
         e.flat();
         //rf2(); 
      w.Stop(); 
      time = w.CpuTime()*1.E9/n;
      cout << "Time using CLHEP::Hard288  =\t " <<   time << "\tns/call" << endl;
   }
   { 
      DualRand e; 
      RandFlat rf2(e);
      w.Start();
      for (int i = 0; i < n; ++i) 
         e.flat();
      //rf2(); 
      w.Stop(); 
      time = w.CpuTime()*1.E9/n;
      cout << "Time using CLHEP::DualRand  =\t " <<   time << "\tns/call" << endl;
   }
   {
      TripleRand e; 
      RandFlat rf2(e);
      w.Start();
      for (int i = 0; i < n; ++i) 
         e.flat();
      //rf2(); 
      w.Stop(); 
      time = w.CpuTime()*1.E9/n;
      cout << "Time using CLHEP::TripleRand  =\t " <<   time << "\tns/call" << endl;
   }
   {
      RanshiEngine e; 
      RandFlat rf2(e);
      w.Start();
      for (int i = 0; i < n; ++i) 
         //rf2(); 
         e.flat();
      w.Stop(); 
      time = w.CpuTime()*1.E9/n;
      cout << "Time using CLHEP::Runshi  =\t " <<   time << "\tns/call" << endl;
   }
   {
      RanluxEngine e; 
      RandFlat rf2(e);
      w.Start();
      for (int i = 0; i < n; ++i) 
         //rf2(); 
         e.flat();
      w.Stop(); 
      time = w.CpuTime()*1.E9/n;
      cout << "Time using CLHEP::RunLux  =\t " <<   time << "\tns/call" << endl;
   }
   {
      Ranlux64Engine e; 
      RandFlat rf2(e);
      w.Start();
      for (int i = 0; i < n; ++i) 
         e.flat(); 
      w.Stop(); 
      time = w.CpuTime()*1.E9/n;
      cout << "Time using CLHEP::RunLux64  =\t " <<   time << "\tns/call" << endl;
   }
   {
      HepJamesRandom e; 
      RandFlat rf2(e);
      w.Start();
      for (int i = 0; i < n; ++i) 
         //rf2(); 
         e.flat(); 
      w.Stop(); 
      time = w.CpuTime()*1.E9/n;
      cout << "Time using CLHEP::HepJames  =\t " <<   time << "\tns/call" << endl;
   }

#endif


   // test the quality by looking at the cdf
   cout <<"\n\nTest quality of Unuran arou" << endl;
   if (! unr.Init( "normal()", "method=arou") ) {
      cout << "Error initializing unuran" << endl;
      return -1;
   }

   TH1D * h1 = new TH1D("h1","cdf on the data ",1000,0,1);
   for (int i = 0; i < 1000000; ++i) {
      double x = unr.Sample();
      h1->Fill( ROOT::Math::normal_cdf( x , 1.0) ); 
   }

   new TCanvas("c1_unuranGaus","unuran Gaus CDF");

   h1->Fit("pol0","Q");
   TF1 * f = h1->GetFunction("pol0");

   std::cout << "CDF Uniform Fit:  chi2 = " << f->GetChisquare() << " ndf = " << f->GetNDF() << std::endl;
   std::cout << "Fit Prob = " << f->GetProb() << std::endl;
   h1->Draw("E");

   if (f->GetProb() < 1.E-4) { 
      cout << "\nERROR: Test Failed !!!!"; 
      return -1; 
   }


   return 0; 

}

#ifndef __CINT__
int main(int argc, char **argv)
{
   int iret = 0;
   if (argc > 1) { 
      TApplication theApp("App",&argc,argv);
      iret = unuranSimple();
      theApp.Run();
   } 
   else 
      iret = unuranSimple();

   return iret;
}
#endif