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

};



//int iParGen[2][3] = {{0,1,2},{3,4,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 = 0) const { 
		// implement evaluation of single chi2 element

		int NrObs = fChi2_Vec.size();

		//Generating parameter array:
		int iParGen[NrObs][3];
		int ParIndex = 0;
		for(int i = 0; i < NrObs; i++){

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

				if(j == 2){ iParGen[i][j] = j;
				}else{
					if(ParIndex == 2) ParIndex++;
					iParGen[i][j] = ParIndex;
					ParIndex++;
				}
			}
		}

		double pars[NrObs][3]; 

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

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

				pars[i][j] = par[iParGen[i][j]];
			} 
		}

		double gPars[NrObs][3];
		double value = 0; 

		if (g != 0) { 
			for (int i = 0; i < 5; ++i) g[i] = 0; 
		}

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

		while( ipoint >= DataBlock ) {

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

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

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

		if(index == 0){

			value = (*(fChi2_Vec.at(index))).DataElement(pars[index], ipoint, gPars[index]);
			//cout << "Returning chi2(" << index << ") value of point nr." << ipoint << endl;
		}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(pars[index], jpoint, gPars[index]);
			//cout << "Returning chi2(" << index << ") value of point nr." << jpoint << endl;
		}		

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

				g[iParGen[index][i]] = gPars[index][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 {

	   int NrObs = fChi2_Vec.size();

		//Generating parameter array:
		int iParGen[NrObs][3];
		int ParIndex = 0;
		for(int i = 0; i < NrObs; i++){

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

				if(j == 2){ iParGen[i][j] = j;
				}else{
					if(ParIndex == 2) ParIndex++;
					iParGen[i][j] = ParIndex;
					ParIndex++;
				}
			}
		}

    double pars[NrObs][3]; 

      for (int i = 0; i < 3; ++i) pars[0][i] = par[iParGen[0][i] ];

      for (int i = 0; i < 3; ++i) pars[1][i] = par[iParGen[1][i] ];

      double value = 0;

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

	      value += (*fChi2_Vec[i])(pars[i]);
      }
      return value;
   } 

	 vector< const ROOT::Math::FitMethodFunction * > fChi2_Vec;
};


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;
	
	//Computing nr. of parameters
	//const int EPparams = ObsVector.size()*2+1;
	const int EPparams = 5;
	

	ROOT::Fit::Fitter EPfitter;

	double FitterParams[EPparams] = {0.35,-0.03,0.03,0.1,-0.009}; 
	//double FitterParams[EPparams] = {0.,0.,0.,0.,0.}; 
	//double* FitterParams = new double[EPparams];


	//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("Fumili2");

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