#include "R_interface.h"

 
 static ScaleHMM* hmm; // declare as static outside the function because we only need one and this enables memory-cleanup on R_CheckUserInterrupt()
 static double** multiD;

 // ===================================================================================================================================================

 // This function takes parameters from R, creates a univariate HMM object, creates the distributions, runs the EM and returns the result to R.

 // ===================================================================================================================================================

 void univariate_hmm(int* O, int* T, int* N, int* state_labels, double* size, double* prob, int* maxiter, int* maxtime, double* eps, double* maxPosterior, int* states, double* A, double* proba, double* loglik, double* weights, int* distr_type, double* initial_size, double* initial_prob, double* initial_A, double* initial_proba, bool* use_initial_params, int* num_threads, int* error, int* read_cutoff, int* algorithm, int* verbosity)

 {

 
 	// Define logging level

 // 	FILE* pFile = fopen("chromStar.log", "w");

 // 	Output2FILE::Stream() = pFile;

 //  	FILELog::ReportingLevel() = FILELog::FromString("NONE");

 //  	FILELog::ReportingLevel() = FILELog::FromString("DEBUG2");

 // 	// Parallelization settings

 // 	omp_set_num_threads(*num_threads);

 
 	// Print some information

 	//FILE_LOG(logINFO) << "number of states = " << *N;

 	if (*verbosity>=1) Rprintf("number of states = %d\n", *N);

 	//FILE_LOG(logINFO) << "number of bins = " << *T;

 	if (*verbosity>=1) Rprintf("number of bins = %d\n", *T);

 	if (*maxiter < 0)
 	{

 		//FILE_LOG(logINFO) << "maximum number of iterations = none";

 		if (*verbosity>=1) Rprintf("maximum number of iterations = none\n");

 	} else {

 		//FILE_LOG(logINFO) << "maximum number of iterations = " << *maxiter;

 		if (*verbosity>=1) Rprintf("maximum number of iterations = %d\n", *maxiter);

 	}
 	if (*maxtime < 0)
 	{

 		//FILE_LOG(logINFO) << "maximum running time = none";

 		if (*verbosity>=1) Rprintf("maximum running time = none\n");

 	} else {

 		//FILE_LOG(logINFO) << "maximum running time = " << *maxtime << " sec";

 		if (*verbosity>=1) Rprintf("maximum running time = %d sec\n", *maxtime);

 	}

 	//FILE_LOG(logINFO) << "epsilon = " << *eps;

 	if (*verbosity>=1) Rprintf("epsilon = %g\n", *eps);

 	//FILE_LOG(logDEBUG3) << "observation vector";

 	for (int t=0; t<50; t++) {

 		//FILE_LOG(logDEBUG3) << "O["<<t<<"] = " << O[t];

 	}
 

 	// Flush if (*verbosity>=1) Rprintf statements to console

 	R_FlushConsole();
 

 	// Create the HMM

 	//FILE_LOG(logDEBUG1) << "Creating a univariate HMM";
 	hmm = new ScaleHMM(*T, *N);

 // 	LogHMM* hmm = new LogHMM(*T, *N);

 	hmm->set_cutoff(*read_cutoff);

 	// Initialize the transition probabilities and proba

 	hmm->initialize_transition_probs(initial_A, *use_initial_params);
 	hmm->initialize_proba(initial_proba, *use_initial_params);

     
 	// Calculate mean and variance of data
 	double mean = 0, variance = 0;
 	for(int t=0; t<*T; t++)
 	{
 		mean+= O[t];
 	}
 	mean = mean / *T;
 	for(int t=0; t<*T; t++)
 	{
 		variance+= pow(O[t] - mean, 2);
 	}
 	variance = variance / *T;

 	//FILE_LOG(logINFO) << "data mean = " << mean << ", data variance = " << variance;		

 	if (*verbosity>=1) Rprintf("data mean = %g, data variance = %g\n", mean, variance);		

 	// Create the emission densities and initialize

 	for (int i_state=0; i_state<*N; i_state++)
 	{
 		if (distr_type[i_state] == 1)

 		{

 			//FILE_LOG(logDEBUG1) << "Using delta distribution for state " << i_state;

 			ZeroInflation *d = new ZeroInflation(O, *T); // delete is done inside ~ScaleHMM()

 			hmm->densityFunctions.push_back(d);
 		}

 		else if (distr_type[i_state] == 2)

 		{

 			//FILE_LOG(logDEBUG1) << "Using geometric distribution for state " << i_state;

 			Geometric *d = new Geometric(O, *T, initial_prob[i_state]); // delete is done inside ~ScaleHMM()

 			hmm->densityFunctions.push_back(d);
 		}

 		else if (distr_type[i_state] == 3)

 		{

 			//FILE_LOG(logDEBUG1) << "Using negative binomial for state " << i_state;

 			NegativeBinomial *d = new NegativeBinomial(O, *T, initial_size[i_state], initial_prob[i_state]); // delete is done inside ~ScaleHMM()

 			hmm->densityFunctions.push_back(d);
 		}

 		else if (distr_type[i_state] == 4)
 		{

 			//FILE_LOG(logDEBUG1) << "Using binomial for state " << i_state;

 			NegativeBinomial *d = new NegativeBinomial(O, *T, initial_size[i_state], initial_prob[i_state]); // delete is done inside ~ScaleHMM()
 			hmm->densityFunctions.push_back(d);
 		}

 		else
 		{

 			//FILE_LOG(logWARNING) << "Density not specified, using default negative binomial for state " << i_state;

 			NegativeBinomial *d = new NegativeBinomial(O, *T, initial_size[i_state], initial_prob[i_state]);

 			hmm->densityFunctions.push_back(d);
 		}

 	}
 

 	// Flush if (*verbosity>=1) Rprintf statements to console

 	R_FlushConsole();
 

 	// Do the EM to estimate the parameters

 	try
 	{

 		if (*algorithm == 1)
 		{
 			hmm->baumWelch();
 		}
 		else if (*algorithm == 3)
 		{
 			//FILE_LOG(logDEBUG1) << "Starting EM estimation";
 			hmm->EM(maxiter, maxtime, eps);
 			//FILE_LOG(logDEBUG1) << "Finished with EM estimation";
 		}

 	}
 	catch (std::exception& e)
 	{

 		//FILE_LOG(logERROR) << "Error in EM/baumWelch: " << e.what();

 		if (*verbosity>=1) Rprintf("Error in EM/baumWelch: %s\n", e.what());

 		if (strcmp(e.what(),"nan detected")==0) { *error = 1; }

 		else { *error = 2; }
 	}

 	// // Compute the posteriors and save results directly to the R pointer
 	// //FILE_LOG(logDEBUG1) << "Recode posteriors into column representation";
 	// #pragma omp parallel for
 	// for (int iN=0; iN<*N; iN++)
 	// {
 	// 	for (int t=0; t<*T; t++)
 	// 	{
 	// 		posteriors[t + iN * (*T)] = hmm->get_posterior(iN, t);
 	// 	}
 	// }

 
 	// Compute the states from posteriors

 	//FILE_LOG(logDEBUG1) << "Computing states from posteriors";
 	int ind_max;
 	std::vector<double> posterior_per_t(*N);

 	for (int t=0; t<*T; t++)

 	{

 		for (int iN=0; iN<*N; iN++)

 		{

 			posterior_per_t[iN] = hmm->get_posterior(iN, t);

 		}

 		ind_max = std::distance(posterior_per_t.begin(), std::max_element(posterior_per_t.begin(), posterior_per_t.end()));
 		states[t] = state_labels[ind_max];

 		maxPosterior[t] = posterior_per_t[ind_max];

 	}
 

 	//FILE_LOG(logDEBUG1) << "Return parameters";

 	// also return the estimated transition matrix and the initial probs

 	for (int i=0; i<*N; i++)

 	{

 		proba[i] = hmm->get_proba(i);
 		for (int j=0; j<*N; j++)

 		{

 			A[i * (*N) + j] = hmm->get_A(j,i);

 		}
 	}
 
 	// copy the estimated distribution params

 	for (int i=0; i<*N; i++)

 	{

 		if (hmm->densityFunctions[i]->get_name() == NEGATIVE_BINOMIAL) 
 		{
 			NegativeBinomial* d = (NegativeBinomial*)(hmm->densityFunctions[i]);
 			size[i] = d->get_size();
 			prob[i] = d->get_prob();
 		}

 		else if (hmm->densityFunctions[i]->get_name() == GEOMETRIC)
 		{
 			Geometric* d = (Geometric*)(hmm->densityFunctions[i]);
 			size[i] = 0;
 			prob[i] = d->get_prob();
 		}

 		else if (hmm->densityFunctions[i]->get_name() == ZERO_INFLATION)
 		{
 			// These values for a Negative Binomial define a zero-inflation (delta distribution)
 			size[i] = 0;
 			prob[i] = 1;
 		}

 		else if (hmm->densityFunctions[i]->get_name() == BINOMIAL) 
 		{
 			Binomial* d = (Binomial*)(hmm->densityFunctions[i]);
 			size[i] = d->get_size();
 			prob[i] = d->get_prob();
 		}

 	}

 	*loglik = hmm->get_logP();
 	hmm->calc_weights(weights);

 	//FILE_LOG(logDEBUG1) << "Deleting the hmm";

 	delete hmm;

 	hmm = NULL; // assign NULL to defuse the additional delete in on.exit() call

 }
 

 
 // =====================================================================================================================================================

 // This function takes parameters from R, creates a multivariate HMM object, runs the EM and returns the result to R.

 // =====================================================================================================================================================

 void multivariate_hmm(double* D, int* T, int* N, int *Nmod, int* comb_states, int* maxiter, int* maxtime, double* eps, double* maxPosterior, int* states, double* A, double* proba, double* loglik, double* initial_A, double* initial_proba, bool* use_initial_params, int* num_threads, int* error, int* algorithm, int* verbosity)

 {
 
 	// Define logging level {"ERROR", "WARNING", "INFO", "ITERATION", "DEBUG", "DEBUG1", "DEBUG2", "DEBUG3", "DEBUG4"}
 // 	FILE* pFile = fopen("chromStar.log", "w");
 // 	Output2FILE::Stream() = pFile;
 //  	FILELog::ReportingLevel() = FILELog::FromString("ITERATION");
 //  	FILELog::ReportingLevel() = FILELog::FromString("DEBUG2");

 //  	FILELog::ReportingLevel() = FILELog::FromString("ERROR");

 
 	// Parallelization settings

 // 	omp_set_num_threads(*num_threads);

 
 	// Print some information

 	//FILE_LOG(logINFO) << "number of states = " << *N;

 	if (*verbosity>=1) Rprintf("number of states = %d\n", *N);

 	//FILE_LOG(logINFO) << "number of bins = " << *T;

 	if (*verbosity>=1) Rprintf("number of bins = %d\n", *T);

 	if (*maxiter < 0)
 	{

 		//FILE_LOG(logINFO) << "maximum number of iterations = none";

 		if (*verbosity>=1) Rprintf("maximum number of iterations = none\n");

 	} else {

 		//FILE_LOG(logINFO) << "maximum number of iterations = " << *maxiter;

 		if (*verbosity>=1) Rprintf("maximum number of iterations = %d\n", *maxiter);

 	}
 	if (*maxtime < 0)
 	{

 		//FILE_LOG(logINFO) << "maximum running time = none";

 		if (*verbosity>=1) Rprintf("maximum running time = none\n");

 	} else {

 		//FILE_LOG(logINFO) << "maximum running time = " << *maxtime << " sec";

 		if (*verbosity>=1) Rprintf("maximum running time = %d sec\n", *maxtime);

 	}

 	//FILE_LOG(logINFO) << "epsilon = " << *eps;

 	if (*verbosity>=1) Rprintf("epsilon = %g\n", *eps);

 	//FILE_LOG(logINFO) << "number of modifications = " << *Nmod;

 	if (*verbosity>=1) Rprintf("number of modifications = %d\n", *Nmod);

 	// Flush if (*verbosity>=1) Rprintf statements to console

 	R_FlushConsole();
 
 	// Recode the densities vector to matrix representation
 // 	clock_t clocktime = clock(), dtime;

 	multiD = CallocDoubleMatrix(*N, *T);

 	for (int iN=0; iN<*N; iN++)
 	{
 		for (int t=0; t<*T; t++)
 		{
 			multiD[iN][t] = D[iN*(*T)+t];
 		}
 	}
 // 	dtime = clock() - clocktime;

 // 	//FILE_LOG(logDEBUG1) << "recoding densities vector to matrix representation: " << dtime << " clicks";

 
 	// Create the HMM

 	//FILE_LOG(logDEBUG1) << "Creating the multivariate HMM";
 	hmm = new ScaleHMM(*T, *N, *Nmod, multiD);

 	// Initialize the transition probabilities and proba
 	hmm->initialize_transition_probs(initial_A, *use_initial_params);
 	hmm->initialize_proba(initial_proba, *use_initial_params);
 	
 	// Print logproba and A
 // 	for (int iN=0; iN<*N; iN++)
 // 	{

 // 		//FILE_LOG(logDEBUG) << "proba["<<iN<<"] = " <<exp(hmm->logproba[iN]);

 // 		for (int jN=0; jN<*N; jN++)
 // 		{

 // 			//FILE_LOG(logDEBUG) << "A["<<iN<<"]["<<jN<<"] = " << hmm->A[iN][jN];

 // 		}
 // 	}
 

 	// Do the EM to estimate the parameters

 	try
 	{

 		if (*algorithm == 1)
 		{
 			hmm->baumWelch();
 		}
 		else if (*algorithm == 3)
 		{
 			//FILE_LOG(logDEBUG1) << "Starting EM estimation";
 			hmm->EM(maxiter, maxtime, eps);
 			//FILE_LOG(logDEBUG1) << "Finished with EM estimation";
 		}

 	}
 	catch (std::exception& e)
 	{

 		//FILE_LOG(logERROR) << "Error in EM/baumWelch: " << e.what();

 		if (*verbosity>=1) Rprintf("Error in EM/baumWelch: %s\n", e.what());

 		if (strcmp(e.what(),"nan detected")==0) { *error = 1; }

 		else { *error = 2; }
 	}
 	
 // 	// Compute the posteriors and save results directly to the R pointer

 // 	//FILE_LOG(logDEBUG1) << "Recode posteriors into column representation";

 // 	for (int iN=0; iN<*N; iN++)
 // 	{
 // 		for (int t=0; t<*T; t++)
 // 		{
 // 			posteriors[t + iN * (*T)] = hmm->get_posterior(iN, t);
 // 		}
 // 	}
 
 	// Compute the states from posteriors

 	//FILE_LOG(logDEBUG1) << "Computing states from posteriors";
 	int ind_max;
 	std::vector<double> posterior_per_t(*N);

 	for (int t=0; t<*T; t++)
 	{
 		for (int iN=0; iN<*N; iN++)
 		{
 			posterior_per_t[iN] = hmm->get_posterior(iN, t);
 		}

 		ind_max = std::distance(posterior_per_t.begin(), std::max_element(posterior_per_t.begin(), posterior_per_t.end()));
 		states[t] = comb_states[ind_max];

 		maxPosterior[t] = posterior_per_t[ind_max];

 	}
 	

 	//FILE_LOG(logDEBUG1) << "Return parameters";

 	// also return the estimated transition matrix and the initial probs
 	for (int i=0; i<*N; i++)
 	{
 		proba[i] = hmm->get_proba(i);
 		for (int j=0; j<*N; j++)
 		{
 				A[i * (*N) + j] = hmm->get_A(j,i);
 		}
 	}
 	*loglik = hmm->get_logP();
 

 	//FILE_LOG(logDEBUG1) << "Deleting the hmm";

 	delete hmm;

 	hmm = NULL; // assign NULL to defuse the additional delete in on.exit() call
 // 	FreeDoubleMatrix(multiD, *N);

 }
 

 
 // =======================================================
 // This function make a cleanup if anything was left over
 // =======================================================

 void univariate_cleanup()

 {
 // 	//FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__; // This message will be shown if interrupt happens before start of C-code
 	delete hmm;
 }
 

 void multivariate_cleanup(int* N)

 {
 	delete hmm;
 	FreeDoubleMatrix(multiD, *N);
 }
 

 
 // ====================================================================
 // C version of apply(array2D, 1, which.max) and apply(array2D, 1, max)
 // ====================================================================
 void array2D_which_max(double* array2D, int* dim, int* ind_max, double* value_max)
 {
   // array2D is actually a vector, but is intended to originate from a 2D array in R
 	std::vector<double> value_per_i0(dim[1]);
   for (int i0=0; i0<dim[0]; i0++)
   {
     for (int i1=0; i1<dim[1]; i1++)
     {
 			value_per_i0[i1] = array2D[i1 * dim[0] + i0];

       // if (*verbosity>=1) Rprintf("i0=%d, i1=%d, value_per_i0[%d] = %g\n", i0, i1, i1, value_per_i0[i1]);

     }
 		ind_max[i0] = 1 + std::distance(value_per_i0.begin(), std::max_element(value_per_i0.begin(), value_per_i0.end()));
     value_max[i0] = *std::max_element(value_per_i0.begin(), value_per_i0.end());
   }
 	

 }