@@ -214,7 +214,7 @@ void univariate_hmm(int* O, int* T, int* N, int* state_labels, double* size, dou

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

 // 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, 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)

+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"}

@@ -330,6 +330,7 @@ void multivariate_hmm(double* D, int* T, int* N, int *Nmod, int* comb_states, in

 		}

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

@@ -9,7 +9,7 @@ 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)

+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

@@ -23,34 +23,34 @@ void univariate_hmm(int* O, int* T, int* N, int* state_labels, double* size, dou

 	// Print some information

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

-	Rprintf("number of states = %d\n", *N);

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

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

-	Rprintf("number of bins = %d\n", *T);

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

 	if (*maxiter < 0)

 	{

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

-		Rprintf("maximum number of iterations = none\n");

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

 	} else {

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

-		Rprintf("maximum number of iterations = %d\n", *maxiter);

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

 	}

 	if (*maxtime < 0)

 	{

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

-		Rprintf("maximum running time = none\n");

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

 	} else {

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

-		Rprintf("maximum running time = %d sec\n", *maxtime);

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

 	}

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

-	Rprintf("epsilon = %g\n", *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 Rprintf statements to console

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

 	R_FlushConsole();

 	// Create the HMM

@@ -75,7 +75,7 @@ void univariate_hmm(int* O, int* T, int* N, int* state_labels, double* size, dou

 	}

 	variance = variance / *T;

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

-	Rprintf("data mean = %g, data variance = %g\n", mean, 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++)

@@ -112,7 +112,7 @@ void univariate_hmm(int* O, int* T, int* N, int* state_labels, double* size, dou

 		}

 	}

-	// Flush Rprintf statements to console

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

 	R_FlushConsole();

 	// Do the EM to estimate the parameters

@@ -132,7 +132,7 @@ void univariate_hmm(int* O, int* T, int* N, int* state_labels, double* size, dou

 	catch (std::exception& e)

 	{

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

-		Rprintf("Error in EM/baumWelch: %s\n", 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; }

 	}

@@ -214,7 +214,7 @@ void univariate_hmm(int* O, int* T, int* N, int* state_labels, double* size, dou

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

 // 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, 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)

+void multivariate_hmm(double* D, int* T, int* N, int *Nmod, int* comb_states, int* maxiter, int* maxtime, double* eps, 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"}

@@ -229,31 +229,31 @@ void multivariate_hmm(double* D, int* T, int* N, int *Nmod, int* comb_states, in

 	// Print some information

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

-	Rprintf("number of states = %d\n", *N);

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

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

-	Rprintf("number of bins = %d\n", *T);

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

 	if (*maxiter < 0)

 	{

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

-		Rprintf("maximum number of iterations = none\n");

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

 	} else {

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

-		Rprintf("maximum number of iterations = %d\n", *maxiter);

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

 	}

 	if (*maxtime < 0)

 	{

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

-		Rprintf("maximum running time = none\n");

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

 	} else {

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

-		Rprintf("maximum running time = %d sec\n", *maxtime);

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

 	}

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

-	Rprintf("epsilon = %g\n", *eps);

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

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

-	Rprintf("number of modifications = %d\n", *Nmod);

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

-	// Flush Rprintf statements to console

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

 	R_FlushConsole();

 	// Recode the densities vector to matrix representation

@@ -303,7 +303,7 @@ void multivariate_hmm(double* D, int* T, int* N, int *Nmod, int* comb_states, in

 	catch (std::exception& e)

 	{

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

-		Rprintf("Error in EM/baumWelch: %s\n", 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; }

 	}

@@ -379,10 +379,10 @@ void array2D_which_max(double* array2D, int* dim, int* ind_max, double* value_ma

     for (int i1=0; i1<dim[1]; i1++)

     {

 			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]);

+      // 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());

   }

-}

\ No newline at end of file

+}

@@ -9,7 +9,7 @@ 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, 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)

+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)

 {

 	// Define logging level

@@ -137,16 +137,16 @@ void univariate_hmm(int* O, int* T, int* N, int* state_labels, double* size, dou

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

@@ -160,6 +160,7 @@ void univariate_hmm(int* O, int* T, int* N, int* state_labels, double* size, dou

 		}

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

@@ -365,3 +366,23 @@ void multivariate_cleanup(int* N)

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

+      // 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());

+  }

+	

+}

\ No newline at end of file

@@ -9,8 +9,7 @@ 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.

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

-extern "C" {

-void R_univariate_hmm(int* O, int* T, int* N, int* state_labels, double* size, double* prob, int* maxiter, int* maxtime, double* eps, 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)

+void univariate_hmm(int* O, int* T, int* N, int* state_labels, double* size, double* prob, int* maxiter, int* maxtime, double* eps, 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)

 {

 	// Define logging level

@@ -209,14 +208,12 @@ void R_univariate_hmm(int* O, int* T, int* N, int* state_labels, double* size, d

 	delete hmm;

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

 }

-} //extern

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

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

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

-extern "C" {

-void R_multivariate_hmm(double* D, int* T, int* N, int *Nmod, int* comb_states, int* maxiter, int* maxtime, double* eps, 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)

+void multivariate_hmm(double* D, int* T, int* N, int *Nmod, int* comb_states, int* maxiter, int* maxtime, double* eps, 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)

 {

 	// Define logging level {"ERROR", "WARNING", "INFO", "ITERATION", "DEBUG", "DEBUG1", "DEBUG2", "DEBUG3", "DEBUG4"}

@@ -351,25 +348,20 @@ void R_multivariate_hmm(double* D, int* T, int* N, int *Nmod, int* comb_states,

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

 // 	FreeDoubleMatrix(multiD, *N);

 }

-} //extern

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

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

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

-extern "C" {

-void R_univariate_cleanup()

+void univariate_cleanup()

 {

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

 	delete hmm;

 }

-} //extern

-extern "C" {

-void R_multivariate_cleanup(int* N)

+void multivariate_cleanup(int* N)

 {

 	delete hmm;

 	FreeDoubleMatrix(multiD, *N);

 }

-} //extern

@@ -20,7 +20,7 @@ void R_univariate_hmm(int* O, int* T, int* N, int* state_labels, double* size, d

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

 // 	// Parallelization settings

-	omp_set_num_threads(*num_threads);

+// 	omp_set_num_threads(*num_threads);

 	// Print some information

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

@@ -227,7 +227,7 @@ void R_multivariate_hmm(double* D, int* T, int* N, int *Nmod, int* comb_states,

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

 	// Parallelization settings

-	omp_set_num_threads(*num_threads);

+// 	omp_set_num_threads(*num_threads);

 	// Print some information

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

@@ -1,11 +1,7 @@

-#include "utility.h"

-#include "scalehmm.h"

-#include "loghmm.h"

-#include <omp.h> // parallelization options

-#include <string> // strcmp

+#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;

@@ -213,7 +209,7 @@ void R_univariate_hmm(int* O, int* T, int* N, int* state_labels, double* size, d

 	delete hmm;

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

 }

-} // extern C

+} //extern

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

@@ -355,7 +351,7 @@ void R_multivariate_hmm(double* D, int* T, int* N, int *Nmod, int* comb_states,

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

 // 	FreeDoubleMatrix(multiD, *N);

 }

-} // extern C

+} //extern

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

@@ -367,7 +363,7 @@ void R_univariate_cleanup()

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

 	delete hmm;

 }

-}

+} //extern

 extern "C" {

 void R_multivariate_cleanup(int* N)

@@ -375,5 +371,5 @@ void R_multivariate_cleanup(int* N)

 	delete hmm;

 	FreeDoubleMatrix(multiD, *N);

 }

-}

+} //extern

@@ -11,10 +11,10 @@ static ScaleHMM* hmm; // declare as static outside the function because we only

 static double** multiD;

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

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

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

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

 extern "C" {

-void R_univariate_hmm(int* O, int* T, int* N, int* state_labels, double* size, double* prob, int* maxiter, int* maxtime, double* eps, 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)

+void R_univariate_hmm(int* O, int* T, int* N, int* state_labels, double* size, double* prob, int* maxiter, int* maxtime, double* eps, 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)

 {

 	// Define logging level

@@ -120,21 +120,27 @@ void R_univariate_hmm(int* O, int* T, int* N, int* state_labels, double* size, d

 	// Flush Rprintf statements to console

 	R_FlushConsole();

-	// Do the Baum-Welch to estimate the parameters

-	//FILE_LOG(logDEBUG1) << "Starting Baum-Welch estimation";

+	// Do the EM to estimate the parameters

 	try

 	{

-		hmm->baumWelch(maxiter, maxtime, eps);

+		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 Baum-Welch: " << e.what();

-		Rprintf("Error in Baum-Welch: %s\n", e.what());

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

+		Rprintf("Error in EM/baumWelch: %s\n", e.what());

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

 		else { *error = 2; }

 	}

-		

-	//FILE_LOG(logDEBUG1) << "Finished with Baum-Welch estimation";

 // 	// Compute the posteriors and save results directly to the R pointer

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

@@ -211,10 +217,10 @@ void R_univariate_hmm(int* O, int* T, int* N, int* state_labels, double* size, d

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

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

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

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

 extern "C" {

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

+void R_multivariate_hmm(double* D, int* T, int* N, int *Nmod, int* comb_states, int* maxiter, int* maxtime, double* eps, 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)

 {

 	// Define logging level {"ERROR", "WARNING", "INFO", "ITERATION", "DEBUG", "DEBUG1", "DEBUG2", "DEBUG3", "DEBUG4"}

@@ -286,20 +292,27 @@ void R_multivariate_hmm(double* D, int* T, int* N, int *Nmod, int* comb_states,

 // 		}

 // 	}

-	// Estimate the parameters

-	//FILE_LOG(logDEBUG1) << "Starting Baum-Welch estimation";

+	// Do the EM to estimate the parameters

 	try

 	{

-		hmm->baumWelch(maxiter, maxtime, eps);

+		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 Baum-Welch: " << e.what();

-		Rprintf("Error in Baum-Welch: %s\n", e.what());

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

+		Rprintf("Error in EM/baumWelch: %s\n", e.what());

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

 		else { *error = 2; }

 	}

-	//FILE_LOG(logDEBUG1) << "Finished with Baum-Welch estimation";

 // 	// Compute the posteriors and save results directly to the R pointer

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

@@ -1,18 +1,3 @@

-//aneufinder - An R-package for CNV detection in whole-genome single cell sequencing data

-//Copyright (C) 2015  Aaron Taudt

-//

-//This program is free software: you can redistribute it and/or modify

-//it under the terms of the GNU General Public License as published by

-//the Free Software Foundation, either version 3 of the License, or

-//(at your option) any later version.

-//

-//This program is distributed in the hope that it will be useful,

-//but WITHOUT ANY WARRANTY; without even the implied warranty of

-//MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the

-//GNU General Public License for more details.

-//

-//You should have received a copy of the GNU General Public License

-//along with this program.  If not, see <http://www.gnu.org/licenses/>.

@@ -1,3 +1,21 @@

+//aneufinder - An R-package for CNV detection in whole-genome single cell sequencing data

+//Copyright (C) 2015  Aaron Taudt

+//

+//This program is free software: you can redistribute it and/or modify

+//it under the terms of the GNU General Public License as published by

+//the Free Software Foundation, either version 3 of the License, or

+//(at your option) any later version.

+//

+//This program is distributed in the hope that it will be useful,

+//but WITHOUT ANY WARRANTY; without even the implied warranty of

+//MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the

+//GNU General Public License for more details.

+//

+//You should have received a copy of the GNU General Public License

+//along with this program.  If not, see <http://www.gnu.org/licenses/>.

+

+

+

 #include "utility.h"

 #include "scalehmm.h"

 #include "loghmm.h"

@@ -11,7 +11,7 @@ static double** multiD;

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

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

 extern "C" {

-void R_univariate_hmm(int* O, int* T, int* N, int* state_labels, double* size, double* prob, double* lambda, int* maxiter, int* maxtime, double* eps, int* states, double* A, double* proba, double* loglik, double* weights, int* distr_type, double* initial_size, double* initial_prob, double* initial_lambda, double* initial_A, double* initial_proba, bool* use_initial_params, int* num_threads, int* error, int* read_cutoff)

+void R_univariate_hmm(int* O, int* T, int* N, int* state_labels, double* size, double* prob, int* maxiter, int* maxtime, double* eps, 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)

 {

 	// Define logging level

@@ -101,12 +101,6 @@ void R_univariate_hmm(int* O, int* T, int* N, int* state_labels, double* size, d

 			hmm->densityFunctions.push_back(d);

 		}

 		else if (distr_type[i_state] == 4)

-		{

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

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

-			hmm->densityFunctions.push_back(d);

-		}

-		else if (distr_type[i_state] == 5)

 		{

 			//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()

@@ -196,11 +190,6 @@ void R_univariate_hmm(int* O, int* T, int* N, int* state_labels, double* size, d

 			size[i] = 0;

 			prob[i] = 1;

 		}

-		else if (hmm->densityFunctions[i]->get_name() == POISSON)

-		{

-			Poisson* d = (Poisson*)(hmm->densityFunctions[i]);

-			lambda[i] = d->get_lambda();

-		}

 		else if (hmm->densityFunctions[i]->get_name() == BINOMIAL) 

 		{

 			Binomial* d = (Binomial*)(hmm->densityFunctions[i]);

@@ -1,59 +1,63 @@

 #include "utility.h"

 #include "scalehmm.h"

 #include "loghmm.h"

-// #include <omp.h> // parallelization options

+#include <omp.h> // parallelization options

+#include <string> // strcmp

+

+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 Baum-Welch and returns the result to R.

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

 extern "C" {

-void R_univariate_hmm(int* O, int* T, int* N, int* state_labels, double* size, double* prob, double* lambda, int* maxiter, int* maxtime, double* eps, int* states, double* A, double* proba, double* loglik, double* weights, int* distr_type, int* iniproc, double* initial_size, double* initial_prob, double* initial_lambda, double* initial_A, double* initial_proba, bool* use_initial_params, int* num_threads, int* error, int* read_cutoff)

+void R_univariate_hmm(int* O, int* T, int* N, int* state_labels, double* size, double* prob, double* lambda, int* maxiter, int* maxtime, double* eps, int* states, double* A, double* proba, double* loglik, double* weights, int* distr_type, double* initial_size, double* initial_prob, double* initial_lambda, double* initial_A, double* initial_proba, bool* use_initial_params, int* num_threads, int* error, int* read_cutoff)

 {

 	// Define logging level

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

 // 	Output2FILE::Stream() = pFile;

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

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

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

 // 	// Parallelization settings

-// 	omp_set_num_threads(*num_threads);

+	omp_set_num_threads(*num_threads);

 	// Print some information

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

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

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

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

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

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

 	if (*maxiter < 0)

 	{

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

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

 		Rprintf("maximum number of iterations = none\n");

 	} else {

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

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

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

 	}

 	if (*maxtime < 0)

 	{

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

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

 		Rprintf("maximum running time = none\n");

 	} else {

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

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

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

 	}

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

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

 	Rprintf("epsilon = %g\n", *eps);

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

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

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

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

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

 	}

 	// Flush Rprintf statements to console

 	R_FlushConsole();

 	// Create the HMM

-	FILE_LOG(logDEBUG1) << "Creating a univariate HMM";

-	ScaleHMM* hmm = new ScaleHMM(*T, *N);

+	//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

@@ -72,7 +76,7 @@ void R_univariate_hmm(int* O, int* T, int* N, int* state_labels, double* size, d

 		variance+= pow(O[t] - mean, 2);

 	}

 	variance = variance / *T;

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

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

 	Rprintf("data mean = %g, data variance = %g\n", mean, variance);		

 	// Create the emission densities and initialize

@@ -80,37 +84,37 @@ void R_univariate_hmm(int* O, int* T, int* N, int* state_labels, double* size, d

 	{

 		if (distr_type[i_state] == 1)

 		{

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

+			//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;

+			//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;

+			//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 poisson for state " << i_state;

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

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

 			hmm->densityFunctions.push_back(d);

 		}

 		else if (distr_type[i_state] == 5)

 		{

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

+			//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;

+			//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);

 		}

@@ -120,23 +124,23 @@ void R_univariate_hmm(int* O, int* T, int* N, int* state_labels, double* size, d

 	R_FlushConsole();

 	// Do the Baum-Welch to estimate the parameters

-	FILE_LOG(logDEBUG1) << "Starting Baum-Welch estimation";

+	//FILE_LOG(logDEBUG1) << "Starting Baum-Welch estimation";

 	try

 	{

 		hmm->baumWelch(maxiter, maxtime, eps);

 	}

 	catch (std::exception& e)

 	{

-		FILE_LOG(logERROR) << "Error in Baum-Welch: " << e.what();

+		//FILE_LOG(logERROR) << "Error in Baum-Welch: " << e.what();

 		Rprintf("Error in Baum-Welch: %s\n", e.what());

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

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

 		else { *error = 2; }

 	}

-	FILE_LOG(logDEBUG1) << "Finished with Baum-Welch estimation";

+	//FILE_LOG(logDEBUG1) << "Finished with Baum-Welch estimation";

 // 	// Compute the posteriors and save results directly to the R pointer

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

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

 // 	#pragma omp parallel for

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

 // 	{

@@ -147,18 +151,20 @@ void R_univariate_hmm(int* O, int* T, int* N, int* state_labels, double* size, d

 // 	}

 	// Compute the states from posteriors

-	FILE_LOG(logDEBUG1) << "Computing states from posteriors";

-	double posterior_per_t [*N];

+	//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);

 		}

-		states[t] = state_labels[argMax(posterior_per_t, *N)];

+		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];

 	}

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

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

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

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

 	{

@@ -205,8 +211,9 @@ void R_univariate_hmm(int* O, int* T, int* N, int* state_labels, double* size, d

 	*loglik = hmm->get_logP();

 	hmm->calc_weights(weights);

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

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

 	delete hmm;

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

 }

 } // extern C

@@ -223,35 +230,35 @@ void R_multivariate_hmm(double* D, int* T, int* N, int *Nmod, int* comb_states,

 // 	Output2FILE::Stream() = pFile;

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

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

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

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

 	// Parallelization settings

-// 	omp_set_num_threads(*num_threads);

+	omp_set_num_threads(*num_threads);

 	// Print some information

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

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

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

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

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

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

 	if (*maxiter < 0)

 	{

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

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

 		Rprintf("maximum number of iterations = none\n");

 	} else {

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

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

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

 	}

 	if (*maxtime < 0)

 	{

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

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

 		Rprintf("maximum running time = none\n");

 	} else {

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

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

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

 	}

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

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

 	Rprintf("epsilon = %g\n", *eps);

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

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

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

 	// Flush Rprintf statements to console

@@ -259,7 +266,7 @@ void R_multivariate_hmm(double* D, int* T, int* N, int *Nmod, int* comb_states,

 	// Recode the densities vector to matrix representation

 // 	clock_t clocktime = clock(), dtime;

-	double** multiD = allocDoubleMatrix(*N, *T);

+	multiD = CallocDoubleMatrix(*N, *T);

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

 	{

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

@@ -268,11 +275,11 @@ void R_multivariate_hmm(double* D, int* T, int* N, int *Nmod, int* comb_states,

 		}

 	}

 // 	dtime = clock() - clocktime;

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

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

 	// Create the HMM

-	FILE_LOG(logDEBUG1) << "Creating the multivariate HMM";

-	ScaleHMM* hmm = new ScaleHMM(*T, *N, *Nmod, multiD);

+	//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);

@@ -280,30 +287,30 @@ void R_multivariate_hmm(double* D, int* T, int* N, int *Nmod, int* comb_states,

 	// Print logproba and A

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

 // 	{

-// 		FILE_LOG(logDEBUG) << "proba["<<iN<<"] = " <<exp(hmm->logproba[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];

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

 // 		}

 // 	}

 	// Estimate the parameters

-	FILE_LOG(logDEBUG1) << "Starting Baum-Welch estimation";

+	//FILE_LOG(logDEBUG1) << "Starting Baum-Welch estimation";

 	try

 	{

 		hmm->baumWelch(maxiter, maxtime, eps);

 	}

 	catch (std::exception& e)

 	{

-		FILE_LOG(logERROR) << "Error in Baum-Welch: " << e.what();

+		//FILE_LOG(logERROR) << "Error in Baum-Welch: " << e.what();

 		Rprintf("Error in Baum-Welch: %s\n", e.what());

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

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

 		else { *error = 2; }

 	}

-	FILE_LOG(logDEBUG1) << "Finished with Baum-Welch estimation";

+	//FILE_LOG(logDEBUG1) << "Finished with Baum-Welch estimation";

 // 	// Compute the posteriors and save results directly to the R pointer

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

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

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

 // 	{

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

@@ -313,18 +320,20 @@ void R_multivariate_hmm(double* D, int* T, int* N, int *Nmod, int* comb_states,

 // 	}

 	// Compute the states from posteriors

-	FILE_LOG(logDEBUG1) << "Computing states from posteriors";

-	double posterior_per_t [*N];

+	//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);

 		}

-		states[t] = comb_states[argMax(posterior_per_t, *N)];

+		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];

 	}

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

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

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

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

 	{

@@ -336,9 +345,30 @@ void R_multivariate_hmm(double* D, int* T, int* N, int *Nmod, int* comb_states,

 	}

 	*loglik = hmm->get_logP();

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

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

 	delete hmm;

-	freeDoubleMatrix(multiD, *N);

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

+// 	FreeDoubleMatrix(multiD, *N);

 }

 } // extern C

+

+// =======================================================

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

+// =======================================================

+extern "C" {

+void R_univariate_cleanup()