#include "Math/IFunction.h"
#include "Math/IParamFunction.h"
#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 2;
		}

};


//Set of parameters - parameter 2 is chi square
int iparP[3] = {1,2,3};
int iparWL[3] = {4,2,5};

struct GlobalChi2 { 
	
	GlobalChi2( vector< ROOT::Fit::Chi2Function * >& chi2_vec ){
	
		for(unsigned int i = 0; i < chi2_vec.size(); i++){
					      
			fChi2_Vec.push_back( chi2_vec.at(i) );
		}
	}

	double operator() (const double *par) const {

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

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

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

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

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

		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::IMultiGenFunction * > 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 = 6.;

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

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

		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++){

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

	//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::DataRange * > DataRange_vec;

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

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

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

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

		BinData_vec.push_back(new ROOT::Fit::BinData(opt,*(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();
	}
	
	
	GlobalChi2 EPChi2(chi2_vec);

	ROOT::Fit::Fitter EPfitter;

	const int EPparams = 5;
	double FitterParams[EPparams] = {0.8,-0.03,0.03,0.8,-0.03}; 

	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().MinimizerOptions().SetPrintLevel(1);
	EPfitter.Config().MinimizerOptions().SetMinimizerType("GSLMultiFit");

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

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

	result.Print(std::cout);
		
}