src/R_interface.cpp
06e7e30e
 #include "R_interface.h"
6d1b8b48
 
e2cf7b55
 static ScaleHMM* hmm; // declare as static outside the function because we only need one and this enables memory-cleanup on R_CheckUserInterrupt()
 static int** multiO;
 
c1d2ec41
 // ===================================================================================================================================================
557a10bd
 // This function takes parameters from R, creates a univariate HMM object, creates the distributions, runs the Baum-Welch and returns the result to R.
c1d2ec41
 // ===================================================================================================================================================
b166b29f
 void univariate_hmm(int* O, int* T, int* N, double* size, double* prob, int* maxiter, int* maxtime, double* eps, double* posteriors, double* densities, bool* keep_densities, int* states, double* maxPosterior, double* A, double* proba, double* loglik, double* weights, int* iniproc, 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* verbosity)
36d28a3c
 {
557a10bd
 
 	// Define logging level
a7bd4f47
 	//FILE* pFile = fopen("chromStar.log", "w");
  	//Output2FILE::Stream() = pFile;
  	//FILELog::ReportingLevel() = FILELog::FromString("ERROR");
  	//FILELog::ReportingLevel() = FILELog::FromString("DEBUG2");
557a10bd
 
632fb5c9
 	//FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;
557a10bd
 	// Parallelization settings
88815e4c
 	#ifdef _OPENMP
9c0d5651
 	omp_set_num_threads(*num_threads);
88815e4c
 	#endif
557a10bd
 
 	// Print some information
632fb5c9
 	//FILE_LOG(logINFO) << "number of states = " << *N;
22b31710
 	if (*verbosity>=1) Rprintf("HMM: number of states = %d\n", *N);
632fb5c9
 	//FILE_LOG(logINFO) << "number of bins = " << *T;
22b31710
 	if (*verbosity>=1) Rprintf("HMM: number of bins = %d\n", *T);
557a10bd
 	if (*maxiter < 0)
 	{
632fb5c9
 		//FILE_LOG(logINFO) << "maximum number of iterations = none";
22b31710
 		if (*verbosity>=1) Rprintf("HMM: maximum number of iterations = none\n");
557a10bd
 	} else {
632fb5c9
 		//FILE_LOG(logINFO) << "maximum number of iterations = " << *maxiter;
22b31710
 		if (*verbosity>=1) Rprintf("HMM: maximum number of iterations = %d\n", *maxiter);
557a10bd
 	}
 	if (*maxtime < 0)
 	{
632fb5c9
 		//FILE_LOG(logINFO) << "maximum running time = none";
22b31710
 		if (*verbosity>=1) Rprintf("HMM: maximum running time = none\n");
557a10bd
 	} else {
632fb5c9
 		//FILE_LOG(logINFO) << "maximum running time = " << *maxtime << " sec";
22b31710
 		if (*verbosity>=1) Rprintf("HMM: maximum running time = %d sec\n", *maxtime);
557a10bd
 	}
632fb5c9
 	//FILE_LOG(logINFO) << "epsilon = " << *eps;
22b31710
 	if (*verbosity>=1) Rprintf("HMM: epsilon = %g\n", *eps);
6cade34f
 
632fb5c9
 	//FILE_LOG(logDEBUG3) << "observation vector";
557a10bd
 	for (int t=0; t<50; t++) {
632fb5c9
 		//FILE_LOG(logDEBUG3) << "O["<<t<<"] = " << O[t];
557a10bd
 	}
 
3c231119
 	// Flush Rprintf statements to console
 	R_FlushConsole();
 
557a10bd
 	// Create the HMM
632fb5c9
 	//FILE_LOG(logDEBUG1) << "Creating a univariate HMM";
22b31710
 	hmm = new ScaleHMM(*T, *N, *verbosity);
62f53d65
 	hmm->set_cutoff(*read_cutoff);
557a10bd
 	// Initialize the transition probabilities and proba
62f53d65
 	hmm->initialize_transition_probs(initial_A, *use_initial_params);
 	hmm->initialize_proba(initial_proba, *use_initial_params);
557a10bd
     
 	// Calculate mean and variance of data
690bc637
 	double Tadjust = 0, mean = 0, variance = 0;
557a10bd
 	for(int t=0; t<*T; t++)
 	{
690bc637
 		if (O[t]>0)
 		{
 			mean += O[t];
 			Tadjust += 1;
 		}
557a10bd
 	}
690bc637
 	mean = mean / Tadjust;
557a10bd
 	for(int t=0; t<*T; t++)
 	{
690bc637
 		if (O[t]>0)
 		{
 			variance += pow(O[t] - mean, 2);
 		}
557a10bd
 	}
690bc637
 	variance = variance / Tadjust;
632fb5c9
 	//FILE_LOG(logINFO) << "data mean = " << mean << ", data variance = " << variance;		
22b31710
 	if (*verbosity>=1) Rprintf("HMM: data mean = %g, data variance = %g\n", mean, variance);		
557a10bd
 	
62f53d65
 	// Go through all states of the hmm and assign the density functions
53067d60
 	double imean=0, ivariance=0;
62f53d65
 	for (int i_state=0; i_state<*N; i_state++)
557a10bd
 	{
 		if (*use_initial_params) {
632fb5c9
 			//FILE_LOG(logINFO) << "Using given parameters for size and prob";
22b31710
 			if (*verbosity>=1) Rprintf("HMM: Using given parameters for size and prob\n");
62f53d65
 			imean = (1-initial_prob[i_state])*initial_size[i_state] / initial_prob[i_state];
 			ivariance = imean / initial_prob[i_state];
632fb5c9
 			//FILE_LOG(logDEBUG2) << "imean = " << imean;
 			//FILE_LOG(logDEBUG2) << "ivariance = " << ivariance;
557a10bd
 		} else {
 
6d87d317
 			if (*iniproc == 1)
a46a7bcc
 			{
6d87d317
 				// Simple initialization, seems to give the fastest convergence
 				if (i_state == 1)
 				{
632fb5c9
 					//FILE_LOG(logDEBUG) << "Initializing size and prob for state 1";
6d87d317
 					imean = mean;
 					ivariance = variance;
 				}
 				else if (i_state == 2)
 				{
632fb5c9
 					//FILE_LOG(logDEBUG) << "Initializing size and prob for state 2";
1e9da246
 					imean = mean*1.5;
6a814487
 					ivariance = variance*2;
6d87d317
 				} 
 				// Make sure variance is greater than mean
 				if (imean >= ivariance)
 				{
 					ivariance = imean + 1;
 				}
 			}
 			else if (*iniproc == 2)
a46a7bcc
 			{
6d87d317
 				// Disturb mean and variance for use as randomized initial parameters
632fb5c9
 				//FILE_LOG(logINFO) << "Using random initialization for size and prob";
22b31710
 				if (*verbosity>=1) Rprintf("HMM: Using random initialization for size and prob\n");
06e7e30e
 					imean = runif(0, 10*mean);
 					ivariance = imean + runif(0, 20*imean); // variance has to be greater than mean, otherwise r will be negative
632fb5c9
 				//FILE_LOG(logDEBUG2) << "imean = " << imean;
 				//FILE_LOG(logDEBUG2) << "ivariance = " << ivariance;
6d87d317
 			}
 			else if (*iniproc == 3)
a46a7bcc
 			{
6d87d317
 				// Empirical initialization
 				if (i_state == 1)
 				{
632fb5c9
 					//FILE_LOG(logINFO) << "Initializing r and p empirically for state 1";
22b31710
 					if (*verbosity>=1) Rprintf("HMM: Initializing r and p empirically for state 1\n");
6d87d317
 					imean = mean/2;
 					ivariance = imean*2;
 				}
 				else if (i_state == 2)
 				{
632fb5c9
 					//FILE_LOG(logINFO) << "Initializing r and p empirically for state 2";
22b31710
 					if (*verbosity>=1) Rprintf("HMM: Initializing r and p empirically for state 2\n");
6d87d317
 					imean = mean*2;
 					ivariance = imean*2;
 				} 
a46a7bcc
 			}
6d87d317
 
557a10bd
 			// Calculate r and p from mean and variance
62f53d65
 			initial_size[i_state] = pow(imean,2)/(ivariance-imean);
 			initial_prob[i_state] = imean/ivariance;
6d87d317
 
557a10bd
 		}
 
8e8433f5
 		if (i_state >= 1)
557a10bd
 		{
632fb5c9
 			//FILE_LOG(logDEBUG1) << "Using negative binomial for state " << i_state;
62f53d65
 			NegativeBinomial *d = new NegativeBinomial(O, *T, initial_size[i_state], initial_prob[i_state]); // delete is done inside ~ScaleHMM()
 			hmm->densityFunctions.push_back(d);
557a10bd
 		}
8e8433f5
 		else if (i_state == 0)
557a10bd
 		{
632fb5c9
 			//FILE_LOG(logDEBUG1) << "Using only zeros for state " << i_state;
62f53d65
 			ZeroInflation *d = new ZeroInflation(O, *T); // delete is done inside ~ScaleHMM()
 			hmm->densityFunctions.push_back(d);
557a10bd
 		}
 		else
 		{
632fb5c9
 			//FILE_LOG(logWARNING) << "Density not specified, using default negative binomial for state " << i_state;
62f53d65
 			NegativeBinomial *d = new NegativeBinomial(O, *T, initial_size[i_state], initial_prob[i_state]);
 			hmm->densityFunctions.push_back(d);
557a10bd
 		}
 	}
 
3c231119
 	// Flush Rprintf statements to console
 	R_FlushConsole();
 
557a10bd
 	// Do the Baum-Welch to estimate the parameters
632fb5c9
 	//FILE_LOG(logDEBUG1) << "Starting Baum-Welch estimation";
36d28a3c
 	try
 	{
62f53d65
 		hmm->baumWelch(maxiter, maxtime, eps);
36d28a3c
 	}
 	catch (std::exception& e)
 	{
632fb5c9
 		//FILE_LOG(logERROR) << "Error in Baum-Welch: " << e.what();
22b31710
 		if (*verbosity>=1) Rprintf("HMM: Error in Baum-Welch: %s\n", e.what());
53067d60
 		if (strcmp(e.what(),"nan detected")==0) { *error = 1; }
a46a7bcc
 		else { *error = 2; }
36d28a3c
 	}
a46a7bcc
 		
632fb5c9
 	//FILE_LOG(logDEBUG1) << "Finished with Baum-Welch estimation";
22b31710
 
 	// Get the posteriors and save results directly to the R pointer
632fb5c9
 	//FILE_LOG(logDEBUG1) << "Recode posteriors into column representation";
22b31710
 	if (*verbosity>=1) Rprintf("HMM: Recoding posteriors ...\n");
987f993a
 	R_FlushConsole();
557a10bd
 	#pragma omp parallel for
62f53d65
 	for (int iN=0; iN<*N; iN++)
557a10bd
 	{
62f53d65
 		for (int t=0; t<*T; t++)
557a10bd
 		{
62f53d65
 			posteriors[t + iN * (*T)] = hmm->get_posterior(iN, t);
557a10bd
 		}
 	}
 
22b31710
 	// Get the densities and save results directly to the R pointer
 	if (*keep_densities == true)
 	{
 		//FILE_LOG(logDEBUG1) << "Recode posteriors into column representation";
 		if (*verbosity>=1) Rprintf("HMM: Recoding densities ...\n");
 		R_FlushConsole();
 		#pragma omp parallel for
 		for (int iN=0; iN<*N; iN++)
 		{
 			for (int t=0; t<*T; t++)
 			{
 				densities[t + iN * (*T)] = hmm->get_density(iN, t);
 			}
 		}
 	}
 
b166b29f
 	// 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] = ind_max + 1;
 		maxPosterior[t] = posterior_per_t[ind_max];
 	}
 		
632fb5c9
 	//FILE_LOG(logDEBUG1) << "Return parameters";
557a10bd
 	// also return the estimated transition matrix and the initial probs
62f53d65
 	for (int i=0; i<*N; i++)
557a10bd
 	{
62f53d65
 		proba[i] = hmm->get_proba(i);
 		for (int j=0; j<*N; j++)
557a10bd
 		{
ddbcdfac
 			A[i * (*N) + j] = hmm->get_A(j,i);
557a10bd
 		}
 	}
 
 	// copy the estimated distribution params
62f53d65
 	for (int i=0; i<*N; i++)
557a10bd
 	{
62f53d65
 		if (hmm->densityFunctions[i]->get_name() == NEGATIVE_BINOMIAL) 
a46a7bcc
 		{
62f53d65
 			NegativeBinomial* d = (NegativeBinomial*)(hmm->densityFunctions[i]);
 			size[i] = d->get_size();
 			prob[i] = d->get_prob();
a46a7bcc
 		}
62f53d65
 		else if (hmm->densityFunctions[i]->get_name() == ZERO_INFLATION)
a46a7bcc
 		{
 			// These values for a Negative Binomial define a zero-inflation (delta distribution)
62f53d65
 			size[i] = 0;
 			prob[i] = 1;
a46a7bcc
 		}
557a10bd
 	}
62f53d65
 	*loglik = hmm->get_logP();
 	hmm->calc_weights(weights);
557a10bd
 	
632fb5c9
 	//FILE_LOG(logDEBUG1) << "Deleting the hmm";
62f53d65
 	delete hmm;
6290af09
 	hmm = NULL; // assign NULL to defuse the additional delete in on.exit() call
557a10bd
 }
 
c1d2ec41
 // =====================================================================================================================================================
557a10bd
 // This function takes parameters from R, creates a multivariate HMM object, creates the distributions, runs the Baum-Welch and returns the result to R.
c1d2ec41
 // =====================================================================================================================================================
e1079a11
 void multivariate_hmm(int* O, int* T, int* N, int *Nmod, double* comb_states, double* size, double* prob, double* w, double* cor_matrix_inv, double* det, int* maxiter, int* maxtime, double* eps, double* posteriors, bool* keep_posteriors, double* densities, bool* keep_densities, double* states, double* maxPosterior, double* A, double* proba, double* loglik, double* initial_A, double* initial_proba, bool* use_initial_params, int* num_threads, int* error, int* verbosity)
36d28a3c
 {
557a10bd
 
 	// Define logging level {"ERROR", "WARNING", "INFO", "ITERATION", "DEBUG", "DEBUG1", "DEBUG2", "DEBUG3", "DEBUG4"}
a7bd4f47
  	//FILE* pFile = fopen("chromStar.log", "w");
 	//Output2FILE::Stream() = pFile;
  	//FILELog::ReportingLevel() = FILELog::FromString("ERROR");
06e7e30e
  	//FILELog::ReportingLevel() = FILELog::FromString("DEBUG3");
557a10bd
 
632fb5c9
 	//FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;
557a10bd
 	// Parallelization settings
88815e4c
 	#ifdef _OPENMP
9c0d5651
 	omp_set_num_threads(*num_threads);
88815e4c
 	#endif
557a10bd
 
 	// Print some information
632fb5c9
 	//FILE_LOG(logINFO) << "number of states = " << *N;
22b31710
 	if (*verbosity>=1) Rprintf("HMM: number of states = %d\n", *N);
632fb5c9
 	//FILE_LOG(logINFO) << "number of bins = " << *T;
22b31710
 	if (*verbosity>=1) Rprintf("HMM: number of bins = %d\n", *T);
557a10bd
 	if (*maxiter < 0)
 	{
632fb5c9
 		//FILE_LOG(logINFO) << "maximum number of iterations = none";
22b31710
 		if (*verbosity>=1) Rprintf("HMM: maximum number of iterations = none\n");
557a10bd
 	} else {
632fb5c9
 		//FILE_LOG(logINFO) << "maximum number of iterations = " << *maxiter;
22b31710
 		if (*verbosity>=1) Rprintf("HMM: maximum number of iterations = %d\n", *maxiter);
557a10bd
 	}
 	if (*maxtime < 0)
 	{
632fb5c9
 		//FILE_LOG(logINFO) << "maximum running time = none";
22b31710
 		if (*verbosity>=1) Rprintf("HMM: maximum running time = none\n");
557a10bd
 	} else {
632fb5c9
 		//FILE_LOG(logINFO) << "maximum running time = " << *maxtime << " sec";
22b31710
 		if (*verbosity>=1) Rprintf("HMM: maximum running time = %d sec\n", *maxtime);
557a10bd
 	}
632fb5c9
 	//FILE_LOG(logINFO) << "epsilon = " << *eps;
22b31710
 	if (*verbosity>=1) Rprintf("HMM: epsilon = %g\n", *eps);
35dbc0ae
 	//FILE_LOG(logINFO) << "number of experiments = " << *Nmod;
 	if (*verbosity>=1) Rprintf("HMM: number of experiments = %d\n", *Nmod);
557a10bd
 
3c231119
 	// Flush Rprintf statements to console
 	R_FlushConsole();
 
557a10bd
 	// Recode the observation vector to matrix representation
36d28a3c
 // 	clock_t clocktime = clock(), dtime;
e2cf7b55
 	multiO = CallocIntMatrix(*Nmod, *T);
557a10bd
 	for (int imod=0; imod<*Nmod; imod++)
 	{
 		for (int t=0; t<*T; t++)
 		{
 			multiO[imod][t] = O[imod*(*T)+t];
 		}
 	}
36d28a3c
 // 	dtime = clock() - clocktime;
632fb5c9
 // 	//FILE_LOG(logDEBUG1) << "recoding observation vector to matrix representation: " << dtime << " clicks";
557a10bd
 
 	// Create the HMM
632fb5c9
 	//FILE_LOG(logDEBUG1) << "Creating the multivariate HMM";
22b31710
 	hmm = new ScaleHMM(*T, *N, *Nmod, *verbosity);
557a10bd
 	// Initialize the transition probabilities and proba
62f53d65
 	hmm->initialize_transition_probs(initial_A, *use_initial_params);
 	hmm->initialize_proba(initial_proba, *use_initial_params);
557a10bd
 	
 	// Print logproba and A
 // 	for (int iN=0; iN<*N; iN++)
 // 	{
632fb5c9
 // 		//FILE_LOG(logDEBUG) << "proba["<<iN<<"] = " <<exp(hmm->logproba[iN]);
557a10bd
 // 		for (int jN=0; jN<*N; jN++)
 // 		{
632fb5c9
 // 			//FILE_LOG(logDEBUG) << "A["<<iN<<"]["<<jN<<"] = " << hmm->A[iN][jN];
557a10bd
 // 		}
 // 	}
 
83cae045
 	// Prepare the binary_states (univariate) vector: binary_states[N][Nmod], e.g., binary_states[iN][imod] tells me at state comb_states[iN], modification imod is non-enriched (0) or enriched (1)
632fb5c9
 	//FILE_LOG(logDEBUG1) << "Preparing the binary_states vector";
b6f70d74
 	double res;
33410caa
 	bool **binary_states = CallocBoolMatrix(*N, *Nmod);
b6f70d74
 	for (int iN=0; iN < *N; iN++) //for each comb state considered
557a10bd
 	{
b6f70d74
 		res = comb_states[iN];
 		for (int imod=(*Nmod-1); imod >= 0; imod--) //for each modification of this comb state
557a10bd
 		{
b6f70d74
 			binary_states[iN][imod] = (bool)fmod(res,2);
 			res = (res - (double)binary_states[iN][imod]) / 2.0;
557a10bd
 		}
 	}
 
 	/* initialize the distributions */
632fb5c9
 	//FILE_LOG(logDEBUG1) << "Initializing the distributions";
62f53d65
 	for (int iN=0; iN<*N; iN++) //for each combinatorial state
557a10bd
 	{
36d28a3c
 		std::vector <Density*> tempMarginals;            
62f53d65
 		for (int imod=0; imod < *Nmod; imod++) //for each modification
557a10bd
 		{
 			Density *d;
 			if (binary_states[iN][imod]) //construct the marginal density function for modification imod being enriched
 			{
62f53d65
 				d = new NegativeBinomial(multiO[imod], *T, size[2*imod+1], prob[2*imod+1]); // delete is done inside ~MVCopulaApproximation()
557a10bd
 			}
 			else //construct the density function for modification imod being non-enriched
 			{
62f53d65
 				d = new ZiNB(multiO[imod], *T, size[2*imod], prob[2*imod], w[imod]); // delete is done inside ~MVCopulaApproximation()
557a10bd
 			}
 			tempMarginals.push_back(d);
 		}
 		//MVCopulaApproximation *tempMVdens = new MVCopulaApproximation(O, tempMarginals, &(cor_matrix_inv[iN*Nmod*Nmod]), det[iN]);
632fb5c9
 		//FILE_LOG(logDEBUG1) << "Calling MVCopulaApproximation for state " << iN;
62f53d65
 		MVCopulaApproximation *tempMVdens = new MVCopulaApproximation(multiO, *T, tempMarginals, &(cor_matrix_inv[iN**Nmod**Nmod]), det[iN]); // delete is done inside ~ScaleHMM()
 		hmm->densityFunctions.push_back(tempMVdens);
557a10bd
 	}
33410caa
 	FreeBoolMatrix(binary_states, *N);
557a10bd
 	
 	// Estimate the parameters
632fb5c9
 	//FILE_LOG(logDEBUG1) << "Starting Baum-Welch estimation";
36d28a3c
 	try
 	{
62f53d65
 		hmm->baumWelch(maxiter, maxtime, eps);
36d28a3c
 	}
 	catch (std::exception& e)
 	{
632fb5c9
 		//FILE_LOG(logERROR) << "Error in Baum-Welch: " << e.what();
22b31710
 		if (*verbosity>=1) Rprintf("HMM: Error in Baum-Welch: %s\n", e.what());
53067d60
 		if (strcmp(e.what(),"nan detected")==0) { *error = 1; }
a46a7bcc
 		else { *error = 2; }
36d28a3c
 	}
632fb5c9
 	//FILE_LOG(logDEBUG1) << "Finished with Baum-Welch estimation";
557a10bd
 	
22b31710
 	// Get the posteriors and save results directly to the R pointer
ee7a8379
 	if (*keep_posteriors == true)
 	{
632fb5c9
 		//FILE_LOG(logDEBUG1) << "Recode posteriors into column representation";
22b31710
 		if (*verbosity>=1) Rprintf("HMM: Recoding posteriors ...\n");
ee7a8379
 		R_FlushConsole();
 		#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);
 			}
 		}
 	}
83cae045
 
22b31710
 	// Get the densities and save results directly to the R pointer
 	if (*keep_densities == true)
 	{
 		//FILE_LOG(logDEBUG1) << "Recode posteriors into column representation";
 		if (*verbosity>=1) Rprintf("HMM: Recoding densities ...\n");
 		R_FlushConsole();
 		#pragma omp parallel for
 		for (int iN=0; iN<*N; iN++)
 		{
 			for (int t=0; t<*T; t++)
 			{
 				densities[t + iN * (*T)] = hmm->get_density(iN, t);
 			}
 		}
 	}
 
83cae045
 	// Compute the states from posteriors
632fb5c9
 	//FILE_LOG(logDEBUG1) << "Computing states from posteriors";
eb486655
 // 	if (*fdr == -1)
 // 	{
860bc753
 		int ind_max;
53067d60
 		std::vector<double> posterior_per_t(*N);
eb486655
 		for (int t=0; t<*T; t++)
557a10bd
 		{
eb486655
 			for (int iN=0; iN<*N; iN++)
 			{
 				posterior_per_t[iN] = hmm->get_posterior(iN, t);
 			}
860bc753
 			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];
b166b29f
 			maxPosterior[t] = posterior_per_t[ind_max];
557a10bd
 		}
eb486655
 // 	}
 // 	else
 // 	{
 // 		double** transformed_posteriors = CallocDoubleMatrix(*T, *Nmod);
 // 		for (int t=0; t<*T; t++)
 // 		{
 // 			for (int iN=0; iN<*N; iN++)
 // 			{
 // 				for (int iNmod=0; iNmod<*Nmod; iNmod++)
 // 				{
 // 					transformed_posteriors[t][iNmod] += (double)binary_states[iN][iNmod] * hmm->get_posterior(iN, t);
 // 				}
 // 			}
 // 		}
 // 	}
557a10bd
 	
632fb5c9
 	//FILE_LOG(logDEBUG1) << "Return parameters";
557a10bd
 	// also return the estimated transition matrix and the initial probs
62f53d65
 	for (int i=0; i<*N; i++)
557a10bd
 	{
62f53d65
 		proba[i] = hmm->get_proba(i);
 		for (int j=0; j<*N; j++)
557a10bd
 		{
62f53d65
 				A[i * (*N) + j] = hmm->get_A(i,j);
557a10bd
 		}
 	}
62f53d65
 	*loglik = hmm->get_logP();
557a10bd
 
632fb5c9
 	//FILE_LOG(logDEBUG1) << "Deleting the hmm";
62f53d65
 	delete hmm;
6290af09
 	hmm = NULL; // assign NULL to defuse the additional delete in on.exit() call
 // 	FreeIntMatrix(multiO, *Nmod); // free on.exit() in R code
557a10bd
 }
 
 
e2cf7b55
 // =======================================================
 // This function make a cleanup if anything was left over
 // =======================================================
06e7e30e
 void univariate_cleanup()
e2cf7b55
 {
632fb5c9
 // 	//FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__; // This message will be shown if interrupt happens before start of C-code
e2cf7b55
 	delete hmm;
 }
557a10bd
 
06e7e30e
 void multivariate_cleanup(int* Nmod)
e2cf7b55
 {
62f53d65
 	delete hmm;
e2cf7b55
 	FreeIntMatrix(multiO, *Nmod);
557a10bd
 }
e2cf7b55
 
1f2bb2a0
 
 // =======================================================
 // C version of apply(array3D, 1, which.max)
 // =======================================================
 void array3D_which_max(double* array3D, int* dim, int* ind_max)
 {
   // array3D is actually a vector, but is intended to originate from a 3D array in R
 	std::vector<double> value_per_i0(dim[1] * dim[2]);
   for (int i0=0; i0<dim[0]; i0++)
   {
     for (int i1=0; i1<dim[1]; i1++)
     {
       for (int i2=0; i2<dim[2]; i2++)
       {
   			value_per_i0[i1*dim[2]+i2] = array3D[(i1*dim[2]+i2) * dim[0] + i0];
         // Rprintf("i0=%d, i1=%d, i2=%d, value_per_i0[%d] = %g\n", i0, i1, i2, i1*dim[2]+i2, value_per_i0[i1*dim[2]+i2]);
       }
     }
 		ind_max[i0] = 1 + std::distance(value_per_i0.begin(), std::max_element(value_per_i0.begin(), value_per_i0.end()));
   }
 	
 }
 
b166b29f
 // ====================================================================
 // 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)
1f2bb2a0
 {
b166b29f
   // array2D is actually a vector, but is intended to originate from a 2D array in R
 	std::vector<double> value_per_i0(dim[1]);
1f2bb2a0
   for (int i0=0; i0<dim[0]; i0++)
   {
     for (int i1=0; i1<dim[1]; i1++)
     {
b166b29f
 			value_per_i0[i1] = array2D[i1 * dim[0] + i0];
       // Rprintf("i0=%d, i1=%d, value_per_i0[%d] = %g\n", i0, i1, i1, value_per_i0[i1]);
1f2bb2a0
     }
b166b29f
 		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());
1f2bb2a0
   }
 	
 }