// example of using Fumili for fitting with a user - provided chi2 function

#include "Math/FitMethodFunction.h"
#include "Math/Minimizer.h"
#include "Math/Factory.h"
#include "TRandom2.h"
#include "TF1.h"
#include "TGraphErrors.h"
#include "Fit/Chi2FCN.h"
#include "Fit/BinData.h"
#include "Math/WrappedMultiTF1.h"

// the chi2 function needed for using FUmili must inherit from the FitMethodFunction interface
class MyChi2Function : public ROOT::Math::FitMethodFunction { 

public:

   // needed for AClic
   typedef ROOT::Math::BasicFitMethodFunction<ROOT::Math::IMultiGenFunction>::Type_t  Type_t;

   MyChi2Function(int npoint, TF1 * modelFunc, const double *x, const double *y, const double * ey ) : 
      ROOT::Math::FitMethodFunction(modelFunc->GetNpar(),npoint),
      fX(x),
      fY(y),
      fEY(ey),
      fModelFunc(modelFunc)
   {
   }

   Type_t Type() const { return ROOT::Math::FitMethodFunction::kLeastSquare; }



   ROOT::Math::IMultiGenFunction * Clone() const { 
      return new MyChi2Function(NPoints(), fModelFunc, fX, fY, fEY );         
   }

   // evalution of the single data element
   double DataElement(const double *par, unsigned int i, double *g  = 0) const; 


private: 

   // evaluation of the all chi2 
   double DoEval(const double * par) const { 
      int ndata = NPoints(); 
      double chi2 = 0; 
      for (int i = 0; i <  ndata; ++i) { 
         double res = DataElement( par, i);
         chi2 += res*res; 
      }
      return chi2; 
   }

   const double * fX;   // x values 
   const double * fY;   // y values 
   const double * fEY;  // error on y values 
   mutable TF1 * fModelFunc;    // model function

};

const double sigma = 1;

void generateData(const TF1 & modelFunc, int ndata, double *x, double *y, double *ey) { 

   // generate the fitting data according to a given function modelFunc smeared with a gaussian
   // of fixed sigma 
   TRandom2 r;
   double x1,x2;
   modelFunc.GetRange(x1,x2);
   for (int i = 0; i < ndata; ++i) { 
      x[i] = r.Uniform(x1,x2);
      double y0 = modelFunc( x[i] );
      y[i] = y0 + r.Gaus(0,sigma);
      ey[i] = sigma; 
   }
   
}

// calculate the chi2 residual multiplied by its weight + the derivatives of 
//  the residual w.r.t the parameters 
double MyChi2Function::DataElement(const double *par, unsigned int i, double *g ) const { 
   double weight = 1./(fEY[i]*fEY[i]); // sigma is considered the same for all points
   double residual = weight * ( fY[i] - fModelFunc->EvalPar( &fX[i], par) ) ;
   // calculate the derivatives 
   if (g != 0)    { 
      fModelFunc->SetParameters(par);
      fModelFunc->GradientPar( &fX[i], g,1.E-9);
      // need to mutiply by -1*weight
      for (int ipar = 0; ipar < fModelFunc->GetNpar(); ++ipar) { 
         g[ipar] *= -1 * weight;
      }
   }
   return residual;
}

// routine performing the fit

int exampleFumili2(std::string fitter="Fumili2") { 

#ifdef __CINT__
   std::cout << "This macro must be run in compiled mode with ACliC" << std::endl;
   std::cout << ">>>>>  .x exampleFumili.C+ " << std::endl; 
   return -2;
#endif

   // model function to fit 
   TF1 * f1 = new TF1("f1","[0]*sin( [1] * x ) ",0,3.14);
   f1->SetParameters(5,2); // initial parameters

   const int ndata = 100;
   TGraphErrors * gr = new TGraphErrors(ndata); 
   double *x = gr->GetX(); 
   double *y = gr->GetY(); 
   double *ey = gr->GetEY(); 


   // generate the data
   generateData(*f1,ndata,x,y,ey);


   // create the chi2 function
   MyChi2Function chi2(ndata,  f1, x, y, ey);

   // create minimizer class with desired algorithm
   ROOT::Math::Minimizer * min = ROOT::Math::Factory::CreateMinimizer(fitter);
   if (min == 0) return -1; 

   // set minimizer function (chi2) 
   // for minuit it needs to be set before setting the variables
   min->SetFunction(chi2); 

   // set initial parameters (different tan true)
   f1->SetParameters(2,1);

   // set minimizer variables (parameters of model function)
   for (int i = 0; i < f1->GetNpar(); ++i) {
      min->SetVariable(i, f1->GetParName(i), f1->GetParameter(i), 0.1);
   }

   min->SetPrintLevel(1);
   min->SetMaxFunctionCalls(1000);
   
   bool ret = min->Minimize(); 
   if (!ret) { 
      std::cout << "Fitting failed !!!! " << std::endl;
   }
   // set obtained value in  TF1
   f1->SetParameters( min->X() );
   f1->SetParErrors(min->Errors() );

   gr->SetTitle("Fitted data and function");
   gr->Draw("AP");
   f1->Draw("same");

   
   // test now using standard chi2 function 

   // clone original function 
   TF1 * f2 = new TF1(*f1);
   f2->SetParameters(2,1);

   ROOT::Fit::BinData fitData(ndata, x, y, 0, ey); 
   ROOT::Math::WrappedMultiTF1 wfunc(*f2);

   ROOT::Fit::Chi2Function chi2_standard(fitData, wfunc);
   min->Clear(); 
   min->SetFunction(chi2_standard); 


   // set minimizer variables (parameters of model function)
   for (int i = 0; i < f2->GetNpar(); ++i) {
      min->SetVariable(i, f2->GetParName(i), f2->GetParameter(i), 0.1);
   }

   min->SetPrintLevel(1);
   min->SetMaxFunctionCalls(1000);
   
   ret = min->Minimize(); 
   if (!ret) { 
      std::cout << "Fitting failed !!!! " << std::endl;
   }
   // set obtained value in  TF1
   f2->SetParameters( min->X() );
   f2->SetParErrors(min->Errors() );
   f2->SetLineColor(kBlue); 

   f2->Draw("same");


}

int main() { 
   return exampleFumili();
}