//#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 iParGen[2][3] = {{0,1,2},{3,1,4}};


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


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; 
	vector< ROOT::Math::FitMethodFunction * > 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();
	}

	cout << "Nr. of data points: " << DataSize << endl;
	
	const int EPparams = 5;
	
	GlobalChi2 EPChi2(EPparams, DataSize, chi2_vec);

	ROOT::Fit::Fitter EPfitter;

	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.Config().MinimizerOptions().SetMinimizerType("Minuit2");
	//EPfitter.Config().MinimizerOptions().SetMinimizerType("Fumili");

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

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

	result.Print(std::cout);
		
}