/*
 * EvaluationMatlab.cpp
 *
 *  Created on: Aug 5, 2015
 *      Author: jjpalacios
 */

#include "EvaluationMatlab.h"

namespace clbDEB {

//=============================================================================
//		CONSTRUCTORS / INITIALIZERS
//=============================================================================
//-----  Main constructor  ----------------------------------------------------
EvaluationMatlab::EvaluationMatlab(ParameterDB *parameters)
: Evaluation(parameters) {
	this->matlabEngine = NULL;
	this->workspace = "";
	this->model = "";
	this->simulationDays = 0;
	this->numObjectives = 1;
	this->fitnessError = MARE;
	this->matlabEvaluations = 0;
}


//-----  setup  ---------------------------------------------------------------
void EvaluationMatlab::setup(ParameterDB *parameters) {
	// Load the parameters of your parent
	Evaluation::setup(parameters);

	// Load the fitness error function (MAE by default)
	std::string value;
	value = parameters->getStringUpper(MATLAB_ERROR);
	if(value.compare(toUpper(MATLAB_ERROR_MSE))==0) this->fitnessError = MSRE;
	else this->fitnessError = MARE;
}


//-----  init  ----------------------------------------------------------------
void EvaluationMatlab::init(ParameterDB *setup,
		std::vector<ParameterDefinition *> &params) {

	// Read the path to the datasets folder
	this->workspace = setup->getString(MATLAB_DATA_FOLDER);
	if(this->workspace.length() == 0) {
		std::string err = "The data folder has not been found";
		throw ClbException(err);
	}

	// Read the name of the model to use
	this->model = setup->getString(MATLAB_MODEL);
	if(this->model.length() == 0) {
		std::string err = "The model to calibrate has not been found";
		throw ClbException(err);
	}

	// Read the number of objectives to estimate
	this->numObjectives = setup->getInteger(MATLAB_NUM_OBJECTIVES, 1);

	// Read the weights for the objectives
	std::string weigthText = setup->getString(MATLAB_WEIGHTS);
	this->setWeights(weigthText);

	// Read the number of days for simulation
	this->simulationDays = setup->getInteger(MATLAB_SIM_DAYS, -1);

	if(this->simulationDays > 0)
		this->setSimulationDays(this->simulationDays);

	// Store the information for parameter values
	this->parameters = params;

	// Read the name of the dataset to use
	std::string datasetList = setup->getString(MATLAB_DATASET);
	this->datasets = fieldsFromCSV(datasetList);
	if(this->datasets.size() == 0) {
		std::string err = "The datasets to use have not been found";
		throw ClbException(err);
	}

	// Load the history of data for each dataset
	this->loadTargetData();

	// Initialize Matlab engine and load the dataset data
	initMatlab();
}


//-----  initMatlab  ----------------------------------------------------------
void EvaluationMatlab::initMatlab() {
	// Initializes the Matlab engine
	if (!(matlabEngine = engOpen(NULL))) {
		throw ClbException("MATLAB engine cannot be started");
	}

	// Set output buffer
	matlabBuffer[MATLAB_BUFSIZE] = '\0';
	engOutputBuffer(matlabEngine, matlabBuffer, MATLAB_BUFSIZE);

	// Change working directory
	std::string command;
	command = "cd('" + this->workspace + "')";
	engEvalString(matlabEngine, command.c_str());
	if(MATLAB_DEBUG) this->printMatlab(command);

	// Load the dataset Matlab data
	for(size_t i=0; i < datasets.size(); i++) {
		std::string dataFile = this->datasets[i] + MATLAB_DATA;
		std::string initScript = this->model + MATLAB_INIT;
		initScript += "(\'" + dataFile + "\');";
		engEvalString(matlabEngine, initScript.c_str());
		if(MATLAB_DEBUG) this->printMatlab(initScript);

		// Store the global variables as local variables
		std::string command = "global ";
		command += std::string(MATLAB_GLOBAL_CHLA) + " " + MATLAB_GLOBAL_WEIGHT + " ";
		command += std::string(MATLAB_GLOBAL_TEMP) + " " + MATLAB_GLOBAL_LENGTH + " ";
		command += std::string(MATLAB_GLOBAL_DAY) + " " + MATLAB_GLOBAL_NDAYS + " ";
		command += std::string(MATLAB_GLOBAL_SPAWN) + ";";
		engEvalString(matlabEngine, command.c_str());
		if(MATLAB_DEBUG) this->printMatlab(command);

		command = this->datasets[i] + "_" + MATLAB_GLOBAL_CHLA + " = ";
		command += std::string(MATLAB_GLOBAL_CHLA) + ";";
		command += this->datasets[i] + "_" + MATLAB_GLOBAL_TEMP + " = ";
		command += std::string(MATLAB_GLOBAL_TEMP) + ";";
		command += this->datasets[i] + "_" + MATLAB_GLOBAL_LENGTH + " = ";
		command += std::string(MATLAB_GLOBAL_LENGTH) + ";";
		command += this->datasets[i] + "_" + MATLAB_GLOBAL_DAY + " = ";
		command += std::string(MATLAB_GLOBAL_DAY) + ";";
		command += this->datasets[i] + "_" + MATLAB_GLOBAL_WEIGHT + " = ";
		command += std::string(MATLAB_GLOBAL_WEIGHT) + ";";
		command += this->datasets[i] + "_" + MATLAB_GLOBAL_NDAYS + " = ";
		command += std::string(MATLAB_GLOBAL_NDAYS) + ";";
		command += this->datasets[i] + "_" + MATLAB_GLOBAL_SPAWN + " = ";
		command += std::string(MATLAB_GLOBAL_SPAWN) + ";";
		engEvalString(matlabEngine, command.c_str());
		if(MATLAB_DEBUG) this->printMatlab(command);
	}
	this->matlabEvaluations = 0;
}


//-----  loadTargetData  ------------------------------------------------------
void EvaluationMatlab::loadTargetData() {
	std::ifstream inputFile;
	std::string fileName;

	for(size_t i=0; i < datasets.size(); i++) {
		fileName = this->workspace + FSEP + this->datasets[i] + MATLAB_HISTORY;
		inputFile.open(fileName.c_str());

		// Read the data from the history file
		if(!inputFile.is_open()) {
			std::string err = "The history data file has not been found";
			throw ClbException(err);
		}

		std::string line;
		std::vector< std::string > fields;
		std::vector< std::vector< double > > histDataset;
		std::vector< double > values;
		getline(inputFile, line);	// Header line
		while(!inputFile.eof()) {
			getline(inputFile, line);
			fields = fieldsFromCSV(line);
			values.clear();

			for(unsigned int i=0; i < fields.size(); i++) {
				if(fields[i].size() < 1)
					values.push_back(-1.0);
				else
					values.push_back(atof(fields[i].c_str()));
			}
			if(fields.size() > 0)
				histDataset.push_back(values);
		}
		inputFile.close();
		history.push_back(histDataset);
	}
}





//=============================================================================
//		METHODS
//=============================================================================
//-----  setSimulationDays  ---------------------------------------------------
void EvaluationMatlab::setSimulationDays(const int simDays) {
	this->simulationDays = simDays;
	std::string command;
	command = "nDays=" + valueToString(simDays) + ";";
	engEvalString(matlabEngine, command.c_str());
	if(MATLAB_DEBUG) this->printMatlab(command);
}


//-----  evaluatePopulation method  -------------------------------------------
double EvaluationMatlab::evaluatePopulation(SharedVars *svars,
			Population *population, bool reEvaluate) {

	// If we exceed the maximum number of evaluations, restart
	// the Matlab connection to free resources
	if(this->matlabEvaluations >= MATLAB_MAX_EVALS) {
		engClose(this->matlabEngine);
		this->initMatlab();
	}

	std::string script, command;
	std::string resultsVar = "results";
	std::string paramVar = "param";
	std::string yveVar = "yve";
	std::vector< std::vector<double> > fitnessPerDataset(population->size());

	for(size_t dataset=0; dataset < this->datasets.size(); dataset++) {
		script = resultsVar + " = [ ]; \n";
		script += paramVar + " = [ ]; \n";
		script += yveVar + " = [ ]; \n";
		command = std::string(MATLAB_GLOBAL_CHLA) + "=" + this->datasets[dataset] + "_" + MATLAB_GLOBAL_CHLA + ";";
		command += std::string(MATLAB_GLOBAL_TEMP) + "=" + this->datasets[dataset] + "_" + MATLAB_GLOBAL_TEMP + ";";
		command += std::string(MATLAB_GLOBAL_LENGTH) + "=" + this->datasets[dataset] + "_" + MATLAB_GLOBAL_LENGTH + ";";
		command += std::string(MATLAB_GLOBAL_DAY) + "=" + this->datasets[dataset] + "_" + MATLAB_GLOBAL_DAY + ";";
		command += std::string(MATLAB_GLOBAL_WEIGHT) + "=" + this->datasets[dataset] + "_" + MATLAB_GLOBAL_WEIGHT + ";";
		command += std::string(MATLAB_GLOBAL_NDAYS) + "=" + this->datasets[dataset] + "_" + MATLAB_GLOBAL_NDAYS + ";";
		command += std::string(MATLAB_GLOBAL_SPAWN) + "=" + this->datasets[dataset] + "_" + MATLAB_GLOBAL_SPAWN + ";";
		engEvalString(matlabEngine, command.c_str());
		if(MATLAB_DEBUG) this->printMatlab(command);

		for(int i=0; i < population->size(); i++) {
			Individual *individual = population->individual(i);
			// Creates the command for the current individual
			command = "[var1 var2 var3] = " + this->model + MATLAB_RUN + "(";

			// Add the parameters to evaluate
			if(individual->getType() == INTEGER_INDIVIDUAL) {
				IntegerIndividual *indiv =
						dynamic_cast<IntegerIndividual *>(individual);

				for(int i=0; i < indiv->length(); i++) {
					if(i > 0)
						command += ",";
					command += valueToString(
						this->parameters[i]->getOriginalValue(indiv->gene(i)));
				}
			}
			else if(individual->getType() == DOUBLE_INDIVIDUAL) {
				DoubleIndividual *indiv =
						dynamic_cast<DoubleIndividual *>(individual);

				for(int i=0; i < indiv->length(); i++) {
					if(i > 0)
						command += ",";
					command += valueToString(
							this->parameters[i]->getOriginalValue(indiv->gene(i)));
				}
			}
			command += "); \n";
			// Appends the command to a script
			script += command;
			script += resultsVar + "=[" + resultsVar + " var1]; \n";
			script += paramVar + "=[" + paramVar + " var2]; \n";
			script += yveVar + "=[" + yveVar + " var3]; \n";
		}

		// Clear variables
		script += "clear var1 var2 var3;";
		// Run the script with the commands
		engEvalString(matlabEngine, script.c_str());
		if(MATLAB_DEBUG) this->printMatlab(script);

		// Look for feasible solutions
		std::vector<double> yveDiff = checkYve(dataset, population->size(),
				engGetVariable(matlabEngine, yveVar.c_str()));

		// Get the fitness values from the MATLAB engine
		std::vector<double> fitnessValues = computeFitness(dataset, population->size(),
				engGetVariable(matlabEngine, resultsVar.c_str()));

		for(unsigned int i=0; i < fitnessValues.size(); i++) {
			// yveDiff penalizes infeasible solutions
			fitnessPerDataset[i].push_back(fitnessValues[i] * (1 + 100.0 * yveDiff[i]));
		}
	}

	for(size_t i=0; i < fitnessPerDataset.size(); i++) {
		population->setFitness(i, aggregateFitness(fitnessPerDataset[i]));
	}

	this->matlabEvaluations += population->size();
	return 1.0;
}


//-----  EVALUATION  ----------------------------------------------------------
double EvaluationMatlab::evaluate(SharedVars *svars, Individual *individual) {
	// If we exceed the maximum number of evaluations, restart
	// the Matlab connection to free resources
	if(this->matlabEvaluations >= MATLAB_MAX_EVALS) {
		engClose(this->matlabEngine);
		this->initMatlab();
	}

	std::string command;
	std::string resultsVar = "results";
	std::vector<double> fitness;
	double error;

	for(size_t dataset=0; dataset < this->datasets.size(); dataset++) {
		command = std::string(MATLAB_GLOBAL_CHLA) + "=";
		command += this->datasets[dataset] + "_" + MATLAB_GLOBAL_CHLA + ";";
		command += std::string(MATLAB_GLOBAL_TEMP) + "=";
		command += this->datasets[dataset] + "_" + MATLAB_GLOBAL_TEMP + ";";
		command += std::string(MATLAB_GLOBAL_LENGTH) + "=";
		command += this->datasets[dataset] + "_" + MATLAB_GLOBAL_LENGTH + ";";
		command += std::string(MATLAB_GLOBAL_DAY) + "=";
		command += this->datasets[dataset] + "_" + MATLAB_GLOBAL_DAY + ";";
		command += std::string(MATLAB_GLOBAL_WEIGHT) + "=";
		command += this->datasets[dataset] + "_" + MATLAB_GLOBAL_WEIGHT + ";";
		command += std::string(MATLAB_GLOBAL_NDAYS) + "=";
		command += this->datasets[dataset] + "_" + MATLAB_GLOBAL_NDAYS + ";";
		command += std::string(MATLAB_GLOBAL_SPAWN) + "=";
		command += this->datasets[dataset] + "_" + MATLAB_GLOBAL_SPAWN + ";";
		engEvalString(matlabEngine, command.c_str());
		if(MATLAB_DEBUG) this->printMatlab(command);


		// Prepare the command to call the evaluation function in Matlab
		command = resultsVar + "=" + this->model + MATLAB_RUN + "(";

		// Add the parameters to evaluate
		if(individual->getType() == INTEGER_INDIVIDUAL) {
			IntegerIndividual *indiv =
					dynamic_cast<IntegerIndividual *>(individual);

			for(int i=0; i < indiv->length(); i++) {
				if(i > 0)
					command += ",";
				command += valueToString(
					this->parameters[i]->getOriginalValue(indiv->gene(i)));
			}
			command += ");";
		}
		else if(individual->getType() == DOUBLE_INDIVIDUAL) {
			DoubleIndividual *indiv =
					dynamic_cast<DoubleIndividual *>(individual);

			for(int i=0; i < indiv->length(); i++) {
				if(i > 0)
					command += ",";
				command += valueToString(
						this->parameters[i]->getOriginalValue(indiv->gene(i)));
			}
			command += ");";
		}

		// Run the command
		engEvalString(matlabEngine, command.c_str());
		if(MATLAB_DEBUG) this->printMatlab(command);

		error = computeFitness(dataset, 1,
				engGetVariable(matlabEngine, resultsVar.c_str()))[0];
		fitness.push_back(error);
	}
	this->matlabEvaluations++;
	return aggregateFitness(fitness);
}



//-----  computeFitness method  -----------------------------------------------
std::vector<double> EvaluationMatlab::computeFitness(int dataset,
		int numIndividuals, mxArray *matlabResult) const {

	std::vector<double> fitnessArray;
	double* results = mxGetPr(matlabResult);

	size_t nDays = mxGetM(matlabResult);		// num of rows (M)
	size_t nResults = mxGetN(matlabResult) / numIndividuals;// num of cols (N)

	size_t nDataPoints = this->history[dataset].size();

	//-------------------------------------------------------------------------
	// Compute fitness from results
	//-------------------------------------------------------------------------
	for(int ind=0; ind < numIndividuals; ind++) {
		double observed = 0.0;
		double estimated = 0.0;
		double error = 0.0;
		double fitness = 0.0;
		std::stringstream sstr;

		size_t simDay = 0;
		size_t initDay = this->history[dataset][0][0];
		std::vector<double> fitnessPerTarget;

		for (size_t target = 1; target < nResults; target++) {
			fitness = 0.0;
			int count = 0;
			for (size_t i = 1; i < nDataPoints; i++) {
				simDay = this->history[dataset][i][0] - initDay;
				if (simDay > nDays) break;

				observed = this->history[dataset][i][target];
				if(observed < 0)
					continue;

				// MATLAB arrays are column-wise
				estimated = *( results + (ind * nDays * nResults) + (target * nDays) + (simDay - 1) );
				if (isnan(estimated)) {
					throw ClbException("The MATLAB model returned a NaN value");
				}

				// Absolute error
				if(this->fitnessError == MARE) {
					error = std::abs( (observed - estimated) / observed );
					fitness += error;
				}
				// Square error
				else if(this->fitnessError == MSRE) {
					error = (observed - estimated) / observed ;
					error = error * error;
					fitness += error;
				}
				count++;
			}
			if(count > 0)
				fitnessPerTarget.push_back(fitness / count);
			else fitnessPerTarget.push_back(-Infd);
		}

		// Apply a linear combination of objectives to compute the fitness
		fitness = 0.0;
		for(size_t i=0; i < fitnessPerTarget.size(); i++)
			fitness += this->weights[i] * fitnessPerTarget[i];
		fitnessArray.push_back(fitness * 100);
	}
	mxDestroyArray(matlabResult);
	return fitnessArray;
}


//-----  checkYve method  -----------------------------------------------------
std::vector<double> EvaluationMatlab::checkYve(int dataset,
		int numIndividuals,	mxArray *matlabResult) const {

	std::vector<double> diffArray;
	double* results = mxGetPr(matlabResult);

	size_t nDays = mxGetM(matlabResult);		// num of rows (M)

	size_t nDataPoints = this->history[dataset].size();

	//-------------------------------------------------------------------------
	// Compute fitness from results
	//-------------------------------------------------------------------------
	for(int ind=0; ind < numIndividuals; ind++) {
		double yve = 0.0;
		double excess = 0.0;
		std::stringstream sstr;

		size_t simDay = 0;
		size_t initDay = this->history[dataset][0][0];

		for (size_t i = 1; i < nDataPoints; i++) {
			simDay = this->history[dataset][i][0] - initDay;
			if (simDay > nDays) break;

			// MATLAB arrays are column-wise
			yve = *( results + (ind * nDays) + (simDay - 1) );
			if (isnan(yve)) {
				throw ClbException("The MATLAB model returned a NaN value");
			}
			if(yve > YVE_LIMIT) {
				excess = std::max(excess, yve - YVE_LIMIT);
			}
		}
		diffArray.push_back(excess);
	}
	mxDestroyArray(matlabResult);
	return diffArray;
}


//-----  setWeights method  ---------------------------------------------------
void EvaluationMatlab::setWeights(std::string &textArray) {
	if(textArray.length() <= 0) {
		for(int i=0; i < this->numObjectives; i++)
			this->weights.push_back(1.0/numObjectives);
	}
	else {
		std::vector<std::string> fields = fieldsFromCSV(textArray);
		for(size_t i=0; i < fields.size(); i++)
			this->weights.push_back(atof(fields[i].c_str()));
	}
}


//-----  aggregateFitness method  ---------------------------------------------
double EvaluationMatlab::aggregateFitness(std::vector<double> fitnessValues) const {
	double sum=0.0;
	for(size_t i=0; i < fitnessValues.size(); i++)
		sum += fitnessValues[i];
	return sum / fitnessValues.size();
}


//-----  printMatlab method  --------------------------------------------------
void EvaluationMatlab::printMatlab(std::string &command) {
	std::cout << std::endl << "Introducing Matlab command:" << std::endl;
	for(size_t i=0; i < command.length(); i++) {
		std::cout << command[i];
		if(command[i] == ';')
			std::cout << std::endl;
	}
	std::cout << std::endl;
}


}
