#include "Fit/Fitter.h"
#include "Fit/BinData.h"
#include "Fit/Chi2FCN.h"
#include "TH1.h"
#include "TList.h"
#include "Math/WrappedMultiTF1.h"
#include "Math/FitMethodFunction.h"
#include "HFitInterface.h"
#include "TCanvas.h"
#include "TStyle.h"
#include "TGraph2DErrors.h"
#include <fstream>
#include <iostream>
#include <iomanip>
#include <cstdlib>
#include <vector>

using namespace std;

int det_entries(ostringstream& oss){
	
	char string[256];
	char GetLine[1000];

	ifstream file;
	file.open(oss.str().c_str());

	if(file.good() == 0){
		cout << "File does not exist, terminating ..." << endl;
		exit(0);
	}

	int nCounter = -1;
	int nNoLines;

	for(int i = 0; i < 3; i++){
		file.getline(GetLine,1000);
	}

	while(file.good() && !file.eof()){
		for(int i = 0; i < 9; i++){
		file >> string;
		cout << string << " ";
		}
		cout << endl; 
		nCounter++;
	}

	nNoLines = nCounter;
	
	cout << "file has " << nNoLines << " lines"  << endl;

	file.close();

	return nNoLines;
}

void read_fit_file(ostringstream& ossFit, int nEntriesFit, vector< vector<double> >& dMassFit,vector< vector<double> >& dVolumeFit,vector< vector<double> >& dChargeFit,vector< vector<double> >& dErrorFit){

	string Value; 
	char GetLine[1000];

	ifstream file;
	file.open(ossFit.str().c_str());

	if(file.good() == 0){
		cout << "File does not exist, terminating ..." << endl;
		exit(0);
	}

	int nCounter = 0;

	for(int i = 0; i < 3; i++){
		file.getline(GetLine,1000);
	}

	while((nCounter < nEntriesFit)){
		file >> Value;
		dMassFit.back().push_back(atof(Value.c_str()));	
		file >> Value;
		dErrorFit.back().push_back(atof(Value.c_str()));	
		file >> Value;
		dVolumeFit.back().push_back(atof(Value.c_str()));	
		file >> Value;
		dChargeFit.back().push_back(atoi(Value.c_str()));					

		for(int i = 0; i < 5; i++){
		file >> Value;
		}


		nCounter++;
	}

	file.close();

}

void DisplayFitData(vector<string>& ObsVector,vector< vector<double> >& dMassFit,vector< vector<double> >& dVolumeFit,vector< vector<double> >& dChargeFit,vector< vector<double> >& dErrorFit, int ObsNr){

	cout << "Observables Nr.1: " << ObsVector.at(ObsNr) << endl;

	for(int i = 0; i < 12; i++){

		cout << dMassFit.at(ObsNr).at(i) << " " << dVolumeFit.at(ObsNr).at(i) << " " << dChargeFit.at(ObsNr).at(i) << " " << dErrorFit.at(ObsNr).at(i) << endl;
	}
}


class MyParametricFunction: public ROOT::Math::IParametricFunctionMultiDim
{
	private:
		const double* pars;

	public:
		double DoEvalPar(const double* x, const double* p) const
		{
			return (p[0] + p[1]/(p[2]*x[0]) - p[1]/(p[2]*p[2]*x[0]*x[0])*x[1]*x[1]);		
		}

		unsigned int NDim() const
		{
			return 2;
		}

		ROOT::Math::IParametricFunctionMultiDim* Clone() const
		{
			return new MyParametricFunction();
		}

		const double* Parameters() const 
		{
			return pars;
		}

		void SetParameters(const double* p)
		{
			pars = p;
		}

		unsigned int NPar() const
		{
			return 3;
		}

};


//Set of parameters - parameter 2 is chi_t 
//int iParGen[2][3] = {{0,1,2},{3,4,2}};
//int iParGen[1][3] = {{0,1,2}};


/*
class GlobalChi2 : public ROOT::Math::FitMethodFunction { 
	
	public:
		
		GlobalChi2( int dim, int npoints, vector< ROOT::Math::FitMethodFunction * >& chi2_vec ): ROOT::Math::FitMethodFunction(dim,npoints){
			 
			for(unsigned int i = 0; i < chi2_vec.size(); i++){
						      
				cout << "GlogablChi2 constr.: Address of chi2_vec.at(" << i << ") to fChi2_Vec ..." << endl;
				fChi2_Vec.push_back( chi2_vec.at(i) );
			}
			cout << "Addresses assigned ..." << endl;
		}

		ROOT::Math::IMultiGenFunction * Clone() const { 
		// copy using default copy-ctor
		// i.e. function pointer will be copied (and not the functions)
			cout << "IMultiGenfunction default copy-ctor ..." << endl;
			return new GlobalChi2(*this);   
		}

		double DataElement(const double *par, unsigned int ipoint, double* g ) const {

			vector< vector<double> > pVec;
			vector< double > dummyVec;

			//cout << "DataElement(): assigning parameters ..." << endl;

			for (int i = 0; i < 3; ++i){

				dummyVec.push_back( par[iParGen[0][i] ]);
			}
			pVec.push_back(dummyVec);

			
			dummyVec.clear();
			
			for (int i = 0; i < 3; ++i){ 

				dummyVec.push_back( par[iParGen[1][i] ]);
			}
			pVec.push_back(dummyVec);
			

			//cout << "DataElement(): Determining residual element ..." << endl;
			//cout << "DataElement(): Point nr. " << ipoint << endl;

			double g_gen[3];
			double value = 0;

			int index = 0;	
			int DataBlock = (*(fChi2_Vec.at(index))).NPoints();

			while( ipoint > DataBlock ) {

				index++;
				DataBlock += (*(fChi2_Vec.at(index))).NPoints(); 
			}

			//cout << "Size of data block to be substracted from ipoint nr." << ipoint << " : " << DataBlock << endl;
			DataBlock = DataBlock - (*(fChi2_Vec.at(index))).NPoints();
			//cout << "Size of data block to be substracted from ipoint nr." << ipoint << " : " << DataBlock << endl;

			//cout << "ipoint is in data block nr." << index << endl;

			if(index == 0){

				value = (*(fChi2_Vec.at(index))).DataElement(&(pVec[index][0]), ipoint, g_gen);
			}else{

				//unsigned int jpoint = ipoint - (*(fChi2_Vec.at(index))).NPoints();
				unsigned int jpoint = ipoint - DataBlock;

				//cout << "ipoint: " << ipoint << " --> " << "jpoint: " << jpoint << endl;

				assert (jpoint < (*(fChi2_Vec.at(index))).NPoints() );

				value = (*(fChi2_Vec.at(index))).DataElement(&(pVec[index][0]), jpoint, g_gen);
			}		

			for(int i = 0; i < 3; i++){

				//cout << "loop nr." << i << ":"<< endl;
				//cout << "iParGen[" << index << "][" << i << "]: " << iParGen[index][i] << endl; 
				//cout << "g_gen[" << i << "]: " << g_gen[i] << endl;
                                if (g) { 
                                   //cout << "g[" <<iParGen[index][i]<<"]: "; cout << g[iParGen[index][i]] << endl;
                                   g[iParGen[index][i]] = g_gen[i];	
                                }
			}

			return value; 

			//cout << "DataElement(): succeed ..." << endl;
		}

		// needed if want to use Fumili or Fumili2
		virtual Type_t Type() const { return ROOT::Math::FitMethodFunction::kLeastSquare; }

	private:
	
		virtual double DoEval (const double *par) const {

			vector< vector<double> > pVec;
			vector< double > dummyVec;

			//cout << "DoEval(): Assigning parameters ..." << endl;

			for (int i = 0; i < 3; ++i){

				dummyVec.push_back(par[iParGen[0][i] ]);
			}
			pVec.push_back(dummyVec);

			dummyVec.clear();
			
			for (int i = 0; i < 3; ++i){ 

				dummyVec.push_back(par[iParGen[1][i] ]);
			}
			pVec.push_back(dummyVec);
			

			//cout << "DoEval(): evaluating ..." << endl;

			double val = 0;

			for( size_t i = 0; i < fChi2_Vec.size(); i++ ){

				val += (*fChi2_Vec[i])(&(pVec[i][0]));
			}

			return val;
		}
		
	 vector< const ROOT::Math::FitMethodFunction * > fChi2_Vec;
};
*/
// definition of shared parameter
int iparObsA[3] = {0,1,2};

// signal + background function 
int iparObsB[3] = {3,4,2};

class GlobalChi2 : public ROOT::Math::FitMethodFunction { 

public:
   GlobalChi2( int dim, int npoints, ROOT::Math::FitMethodFunction & f1,  ROOT::Math::FitMethodFunction & f2) :
      ROOT::Math::FitMethodFunction(dim,npoints),
      fChi2_1(&f1), fChi2_2(&f2) {}


   ROOT::Math::IMultiGenFunction * Clone() const { 
      // copy using default copy-ctor
      // i.e. function pointer will be copied (and not the functions)
      return new GlobalChi2(*this);   
   }

   double  DataElement(const double *par, unsigned int ipoint, double *g = 0) const { 
      // implement evaluation of single chi2 element
      double p1[3];
      for (int i = 0; i < 3; ++i) p1[i] = par[iparObsA[i] ];

      double p2[3]; 
      for (int i = 0; i < 3; ++i) p2[i] = par[iparObsB[i] ];

      double g1[3]; 
      double g2[3];
      double value = 0; 
      if (g != 0) { 
         for (int i = 0; i < 6; ++i) g[i] = 0; 
      }

      if (ipoint < fChi2_1->NPoints() ) {
         value = fChi2_1->DataElement(p1, ipoint, g1);
         // update gradient values
         if (g != 0) { 
            for (int i= 0; i < 3; ++i) g[iparObsA[i]] = g1[i];         
         }
      }
      else { 
         unsigned int jpoint = ipoint- fChi2_1->NPoints();
         assert (jpoint < fChi2_2->NPoints() );
         value =  fChi2_2->DataElement(p2, jpoint, g2);
         // update gradient values
         if ( g != 0) { 
            for (int i= 0; i < 3; ++i) g[iparObsB[i]] = g2[i];
         }
      }
        
      return value; 
      
   }


   // needed if want to use Fumili or Fumili2
   virtual Type_t Type() const { return ROOT::Math::FitMethodFunction::kLeastSquare; }

private:
   // parameter vector is first background (in common 1 and 2) and then is signal (only in 2)
   virtual double DoEval (const double *par) const {
      double p1[3];
      for (int i = 0; i < 3; ++i) p1[i] = par[iparObsA[i] ];

      double p2[3]; 
      for (int i = 0; i < 3; ++i) p2[i] = par[iparObsB[i] ];

      return (*fChi2_1)(p1) + (*fChi2_2)(p2);
   } 

   const  ROOT::Math::FitMethodFunction * fChi2_1;
   const  ROOT::Math::FitMethodFunction * fChi2_2;
};

void NonLinLSFit() { 

	vector< vector<double> > dMassFit;
	vector< vector<double> > dVolumeFit;
	vector< vector<double> > dChargeFit;
	vector< vector<double> > dErrorFit;

	double qmass = 0.03;
	double beta = 5.;

	ostringstream ossFit;
	ossFit << "fit_mass_table_beta_" << (int)beta << "_qmass_" << fixed << setprecision(2) << qmass << "_tc_";

	ostringstream ossObs;

	//Vector which contains the observable functions to be fitted
	vector<string> ObsVector;
	ObsVector.push_back("Pion");
	//ObsVector.push_back("Pion");
	ObsVector.push_back("WL_exp_1");

	TGraph2DErrors** GR = new TGraph2DErrors*[ObsVector.size()]();

	cout << "Nr of observables: " << ObsVector.size() << endl;

	//Filling the vector with data from the data files and generating the TGraph2DErrors objects
	for(unsigned int i = 0; i < ObsVector.size(); i++){

		dMassFit.push_back(vector<double>());
		dVolumeFit.push_back(vector<double>());		
		dChargeFit.push_back(vector<double>());
		dErrorFit.push_back(vector<double>()); 		

		ossObs << ossFit.str() << ObsVector.at(i) << ".dat"; 
		cout << ossObs.str() << endl;

		int nEntriesFit = det_entries(ossObs);
		cout << "Entries in file: " << nEntriesFit << endl;

		read_fit_file(ossObs, nEntriesFit, dMassFit, dVolumeFit, dChargeFit, dErrorFit);

		GR[i] = new TGraph2DErrors();

		//DisplayFitData(ObsVector, dMassFit, dVolumeFit, dChargeFit, dErrorFit, i);

		ossObs.str("");
	}


	//Setting the points for the graph
	for(unsigned int i = 0; i < ObsVector.size(); i++){

		cout << "Adding data points to GR[" << i << "]" << endl;
		for(unsigned int j = 0; j < dVolumeFit.at(i).size(); j++){

			GR[i]->SetPoint(j,dVolumeFit.at(i).at(j),dChargeFit.at(i).at(j),dMassFit.at(i).at(j));
			GR[i]->SetPointError(j,0,0,dErrorFit.at(i).at(j));

			cout << "GR[" << i << "]->SetPoint(" << j << "," << dVolumeFit.at(i).at(j)  << "," << dChargeFit.at(i).at(j)<< "," << dMassFit.at(i).at(j) << "),  GR[" << i << "]->SetPointError(" << j << ",0,0," << dErrorFit.at(i).at(j) << ")" << endl;
		}
		cout << endl;
	}
	

	//Nr of parameters
	int npar = 3;

	MyParametricFunction* fobj = new MyParametricFunction[ObsVector.size()];

	TF1** MyFunc = new TF1*[ObsVector.size()]();

	vector< ROOT::Math::WrappedMultiTF1 * > WrappedMultiTF1_vec;
	//ROOT::Math::WrappedMultiTF1 EP(*(MyFunc[i]),2);	

	for(unsigned int i = 0; i < ObsVector.size(); i++){

		MyFunc[i] = new TF1("MyFunc",fobj[i],0,1,npar,"MyParametricFunction");

		WrappedMultiTF1_vec.push_back( new ROOT::Math::WrappedMultiTF1(*(MyFunc[i]),2) );
	}

	
	//ROOT::Fit::DataOptions opt;
	vector< ROOT::Fit::DataOptions * > DataOptions_vec;

	vector< ROOT::Fit::DataRange * > DataRange_vec;

	vector< ROOT::Fit::BinData * > BinData_vec;

	vector< ROOT::Fit::Chi2Function * > chi2_vec; 
	//vector< ROOT::Math::FitMethodFunction * > chi2_vec; 
	
	int DataSize = 0;

	for(unsigned int i = 0; i < ObsVector.size(); i++){

		DataOptions_vec.push_back(new ROOT::Fit::DataOptions);

		DataRange_vec.push_back(new ROOT::Fit::DataRange);

		BinData_vec.push_back(new ROOT::Fit::BinData(*(DataOptions_vec.at(i)),*(DataRange_vec.at(i))));

		ROOT::Fit::FillData(*(BinData_vec.at(i)), GR[i]);
		
		chi2_vec.push_back( new ROOT::Fit::Chi2Function(*(BinData_vec.at(i)),*(WrappedMultiTF1_vec.at(i))) );

		DataSize += (*(BinData_vec.at(i))).Size();
	}

	cout << "Nr. of data points: " << DataSize << endl;
	
	const int EPparams = 5;
	//const int EPparams = 3;
	

	ROOT::Fit::Fitter EPfitter;

	double FitterParams[EPparams] = {0.35,-0.03,0.03,0.1,-0.009}; 
	//double FitterParams[EPparams] = {0.4,-0.03,0.003}; 

	//EPfitter.Config().SetParamsSettings(5,FitterParams);
	EPfitter.Config().SetParamsSettings(5,FitterParams);

	//In the following loop i corresponds to the number of the parameter
	
	for(int i = 0; i < EPparams; i++){

		//EPfitter.Config().ParSettings(i).SetLimits(-2.0,2.0);
		EPfitter.Config().ParSettings(i).SetStepSize(0.0001);
	}
	
	EPfitter.Config().ParSettings(0).SetLimits(0.,2.0);
	EPfitter.Config().ParSettings(1).SetLimits(-0.5,0.5);
	EPfitter.Config().ParSettings(2).SetLimits(0.,0.5);
	EPfitter.Config().ParSettings(3).SetLimits(0,2.0);
	EPfitter.Config().ParSettings(4).SetLimits(-0.5,0.5);
	
	EPfitter.Config().MinimizerOptions().SetPrintLevel(1);
	EPfitter.Config().MinimizerOptions().SetMinimizerType("GSLMultiFit");
	//EPfitter.Config().MinimizerOptions().SetMinimizerType("Minuit2");
	//EPfitter.Config().MinimizerOptions().SetMinimizerType("Fumili");

	//GlobalChi2 EPChi2(EPparams, DataSize, chi2_vec);
	GlobalChi2 EPChi2(EPparams,DataSize,*(chi2_vec.at(0)),*(chi2_vec.at(1)));

	EPfitter.FitFCN(EPChi2, 0, DataSize, true);

	ROOT::Fit::FitResult result = EPfitter.Result();

	cout << "Size of fitted data: " << DataSize << endl;

	result.Print(std::cout);
		
}