//C++ includes
#include <string>
#include <vector>
#include <valarray>
#include <iostream>
#include <errno.h> 
#include <math.h>

//Root includes
#include "TBrowser.h"
#include "TCanvas.h"
#include "TFile.h"
#include "TH1D.h"
#include "TH1F.h"
#include "TLeaf.h"
#include "TTree.h"
#include "TMatrixD.h"
#include "TGraph.h"
#include "TPad.h"
#include "TLatex.h"
#include "TRandom.h"
#include "TMatrixD.h"
#include "TCanvas.h"
#include "TH2D.h"
#include "TMatrixDEigen.h"
#include "TRandom.h"
#include "TSystem.h"
#include "TStopwatch.h"

//KDE Includes
//#include "../include/KDE.h"
//#include "../include/FileUtils.h"
//#include "../include/StringUtils.h"

using namespace std;
//using namespace pde_tools;

int split(const string& original, const string& delimiter, vector<string>& results)
{
	int numFound = 0;
	int currentPos=0, findPos=0;
	results.clear();
	findPos = original.find( delimiter, 0);
	while( 0 <= findPos){
		assert( currentPos < original.size() );
		if( findPos - currentPos > 1){ //make sure the token value is not empty
			numFound++;
			results.push_back( original.substr(currentPos, findPos - currentPos) );
		}
		//cout<<"token("<<currentPos<<", "<<findPos<<"): "<< original.substr(currentPos, findPos - currentPos) <<endl;
		//cout<<"original.size()="<<original.size()<<endl;
		currentPos = findPos + delimiter.size();
		findPos = original.find( delimiter, currentPos);
	}
	
	if(original.size() - currentPos >= 1){ //add the remainder of the original to the output
		//cout<<"token("<<currentPos<<", "<<findPos<<"): "<< original.substr(currentPos) <<endl;
		results.push_back( original.substr(currentPos) );
		numFound++;
	}

	return numFound;
}//end of split


/////////////////////////////////////////////////////////////////////
void ReadRootTreeIntoTMatrixD(TTree *dataTree, TMatrixD *dataMatrix, const char *variableList){
/////////////////////////////////////////////////////////////////////
/// Extracts data from tree, applies cuts, and copies the values into
/// the TMatrixD object. TMatrixD objects are the data structure used
/// by PDE so this function is useful in converting data from the 
/// the commonly used storage format, TTree, to format used for PDE 
/// processing
/////////////////////////////////////////////////////////////////////
///
/// @param dataTree Tree containing the data to be extracted
/// @param mat TMatrixD object to hold the data from the tree
/// @returns TCut object with all of the cuts for this extractor
/// @bug The cuts should not be hard coded
/////////////////////////////////////////////////////////////////////	
	cout<<"Copying tree "<< dataTree->GetName()<<" to TMatrixD . . ."<<endl;
	vector<string> variablesVec;
	string variableString(variableList);
	Int_t numberVars = split(variableString, ":", variablesVec);
	assert( numberVars == (Int_t)variablesVec.size() );
	printf("%d variables in variableVec\n", variablesVec.size() );
									
	int rows = dataTree->GetEntries();
	int columns = variablesVec.size();
	dataMatrix->Clear();
	dataMatrix->ResizeTo( rows, columns );
	
	for(int i=0; i < rows; ++i){
		Long64_t bytesRead = dataTree->GetEntry(i);
		for(int j=0; j< columns; j++){
			//cout<<"getting "<<variablesVec[j].c_str()<<" . . ."<<endl;
			(*dataMatrix)[i][j] = dataTree->GetLeaf( variablesVec[j].c_str() )->GetValue( );
			//cout<< dataTree->GetLeaf( variablesVec[j].c_str() )->GetValue(0) << endl;
		}
	}
	
	cout << "Built TMatrixD: (" << dataMatrix->GetNrows() << "," << dataMatrix->GetNcols() << ") matrix"<<endl;
} //end of ReadRootTreeIntoTMatrixD	
//////////////////////////////////////////////////////////////////////////////////////////////

///////////////////////////////////////////////////////////////////////////////
/// GetListOfVariables
///////////////////////////////////////////////////////////////////////////////
/// returns a : delimitted list of variables x0:x1:. . .:x[numberOfVariables]
///////////////////////////////////////////////////////////////////////////////
void GetListOfVariables(int numberOfVariables, string &variables){

	char variablesArray[1000], buffer[1000];
	variablesArray[0] = '\0';
	string delimiter("");
	for(int j=0; j<numberOfVariables; ++j){
		sprintf(&buffer[0], "%s%sy%d", variablesArray, delimiter.c_str(), j);
		sprintf(&variablesArray[0], "%s", buffer);
		delimiter = ":";
	}

	variables = variablesArray;
}


///////////////////////////////////////////////////////////////////////////////
/// GetLeavesString
///////////////////////////////////////////////////////////////////////////////
/// returns string formatted to create ROOT branch numberOfVariables leaves.
///////////////////////////////////////////////////////////////////////////////
void GetLeavesString(int numberOfVariables, string &leaves){

	char leavesString[1000], buffer[1000];
	leavesString[0] = '\0';
	string delimiter("");
	for(int j=0; j<numberOfVariables; ++j){
		sprintf(&buffer[0], "%s%sy%d/D", leavesString, delimiter.c_str(), j);
		sprintf(&leavesString[0], "%s", buffer);
		delimiter = ":";
	}
	leaves = leavesString;
}


///////////////////////////////////////////////////////////////////////////////
/// GetLeavesString
///////////////////////////////////////////////////////////////////////////////
/// populates tree with numberOfEvents events,  numberOfVariables, 
/// randomly distributed as Gaussian(mean, std). The random seed can be specified
/// by variable randomSeed.
///////////////////////////////////////////////////////////////////////////////
void GetGaussianTree(TTree &tree, int numberOfVariables, double numberOfEvents, double mean, double std, int randomSeed=12 ){

	printf("Filling tree %s with Gaussian data.\n", tree.GetName() );
	TRandom rand(randomSeed);
	Double_t *leaves = new Double_t[numberOfVariables];
	string leavesString;
	GetLeavesString(numberOfVariables, leavesString);
	printf("%s\n", leavesString.c_str());

	TBranch *branch = tree.Branch("TestPDE", leaves, leavesString.c_str() );

	for(int i=0; i < numberOfEvents; ++i){
		for(int j=0; j<numberOfVariables; ++j){
			leaves[j] = rand.Gaus(mean, std);
		}
		branch->Fill( );
		tree.Fill( );
	}
	
	delete[] leaves;
}

int memError(){

	static const Int_t numberOfTrainingEvents = 3000;
	static const Int_t numberOfTestEvents = 30;
	static const Int_t numberOfVariables = 100;

	cout<<"Testing KDE . . ."<<endl;

	TTree *signalTree = new TTree("signalTree","signal");
	TTree *bgTree = new TTree("bgTree", "background ");
	TTree *testTree = new TTree("dataTree", "test events");
	TMatrixD *signalData = new TMatrixD(numberOfTrainingEvents, numberOfVariables);
	TMatrixD *bgData = new TMatrixD(numberOfTrainingEvents, numberOfVariables);
	TMatrixD *testData = new TMatrixD(numberOfTestEvents , numberOfVariables);

	GetGaussianTree(*signalTree, numberOfVariables, numberOfTrainingEvents, 0.2, 1);
	GetGaussianTree(*bgTree, numberOfVariables, numberOfTrainingEvents, 0, 1);
	GetGaussianTree(*testTree, numberOfVariables, numberOfTestEvents, -0.1, 1);

	//Copy trees into TMatrixD's for input to KDE
	string variableList;
	GetListOfVariables(numberOfVariables, variableList);
	cout<<"vars: "<<variableList<<endl;
	ReadRootTreeIntoTMatrixD( signalTree, signalData, variableList.c_str() );
	cout << "Copied " << signalTree->GetName()<<endl;
	ReadRootTreeIntoTMatrixD( bgTree, bgData, variableList.c_str() );
	cout << "Copied " << bgTree->GetName()<<endl;
	ReadRootTreeIntoTMatrixD( testTree, testData, variableList.c_str() );
	cout << "Copied " << testTree->GetName()<<endl;

	signalData->Delete();
	bgData->Delete();
	testData->Delete();

	testTree->Delete();
	signalTree->Delete();
	bgTree->Delete();
	 
	cout<<"\ntest complete."<<endl;
	return 0;
}


int main(int argc, char **argv){
	gSystem->Load("libTree");
	memError();
}