@@ -273,7 +273,7 @@ void univariate_hmm(int* O, int* T, int* N, double* size, double* prob, int* max

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

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

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

-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, int* 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)

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

 {

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

@@ -6,7 +6,7 @@ static int** multiO;

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

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

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

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

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

 {

 	// Define logging level

@@ -220,6 +220,21 @@ void univariate_hmm(int* O, int* T, int* N, double* size, double* prob, int* max

 		}

 	}

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

+	}

+		

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

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

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

@@ -258,7 +273,7 @@ void univariate_hmm(int* O, int* T, int* N, double* size, double* prob, int* max

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

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

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

-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, int* states, double* A, double* proba, double* loglik, double* initial_A, double* initial_proba, bool* use_initial_params, int* num_threads, int* error, int* verbosity)

+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, int* 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)

 {

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

@@ -432,6 +447,7 @@ void multivariate_hmm(int* O, int* T, int* N, int *Nmod, double* comb_states, do

 			}

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

 		}

 // 	}

 // 	else

@@ -506,46 +522,22 @@ void array3D_which_max(double* array3D, int* dim, int* ind_max)

 }

-

-// ====================================================================================

-// C version of apply(array2D, MARGIN = 1, FUN = mean)

-// ====================================================================================

-void array2D_mean(double* array2D, int* dim, double* mean)

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

+// 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 matrix in R

-  double sum=0;

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

   {

-    sum = 0;

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

     {

-      sum += array2D[i1*dim[0] + i0];

-    }

-    mean[i0] = sum / dim[1];

-  }

-	

-}

-

-

-// ====================================================================================

-// C version of apply(array3D, MARGIN = c(1,2), FUN = mean)

-// ====================================================================================

-void array3D_mean(double* array3D, int* dim, double* mean)

-{

-  // array3D is actually a vector, but is intended to originate from a 3D array in R

-  double sum=0;

-  for (int i0=0; i0<dim[0]; i0++)

-  {

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

-    {

-      sum = 0;

-      for (int i2=0; i2<dim[2]; i2++)

-      {

-        sum += array3D[(i2*dim[1] + i1)*dim[0] + i0];

-        // Rprintf("i0=%d, i1=%d, i2=%d, array3D[(i2*dim[1] + i1)*dim[0] + i0 = %d] = %g\n", i0, i1, i2, (i2*dim[1] + i1)*dim[0] + i0, array3D[(i2*dim[1] + i1)*dim[0] + i0]);

-      }

-      mean[i1*dim[0] + i0] = sum / dim[2];

+			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

@@ -483,3 +483,69 @@ void multivariate_cleanup(int* Nmod)

 	FreeIntMatrix(multiO, *Nmod);

 }

+

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

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

+  }

+	

+}

+

+

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

+// C version of apply(array2D, MARGIN = 1, FUN = mean)

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

+void array2D_mean(double* array2D, int* dim, double* mean)

+{

+  // array2D is actually a vector, but is intended to originate from a matrix in R

+  double sum=0;

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

+  {

+    sum = 0;

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

+    {

+      sum += array2D[i1*dim[0] + i0];

+    }

+    mean[i0] = sum / dim[1];

+  }

+	

+}

+

+

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

+// C version of apply(array3D, MARGIN = c(1,2), FUN = mean)

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

+void array3D_mean(double* array3D, int* dim, double* mean)

+{

+  // array3D is actually a vector, but is intended to originate from a 3D array in R

+  double sum=0;

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

+  {

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

+    {

+      sum = 0;

+      for (int i2=0; i2<dim[2]; i2++)

+      {

+        sum += array3D[(i2*dim[1] + i1)*dim[0] + i0];

+        // Rprintf("i0=%d, i1=%d, i2=%d, array3D[(i2*dim[1] + i1)*dim[0] + i0 = %d] = %g\n", i0, i1, i2, (i2*dim[1] + i1)*dim[0] + i0, array3D[(i2*dim[1] + i1)*dim[0] + i0]);

+      }

+      mean[i1*dim[0] + i0] = sum / dim[2];

+    }

+  }

+	

+}

\ No newline at end of file

@@ -258,7 +258,7 @@ void univariate_hmm(int* O, int* T, int* N, double* size, double* prob, int* max

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

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

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

-void multivariate_hmm(int* O, int* T, int* N, int *Nmod, int* 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, int* states, double* A, double* proba, double* loglik, double* initial_A, double* initial_proba, bool* use_initial_params, int* num_threads, int* error, int* verbosity)

+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, int* states, double* A, double* proba, double* loglik, double* initial_A, double* initial_proba, bool* use_initial_params, int* num_threads, int* error, int* verbosity)

 {

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

@@ -334,14 +334,15 @@ void multivariate_hmm(int* O, int* T, int* N, int *Nmod, int* comb_states, doubl

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

 	//FILE_LOG(logDEBUG1) << "Preparing the binary_states vector";

+	double res;

 	bool **binary_states = CallocBoolMatrix(*N, *Nmod);

-	for(int iN=0; iN < *N; iN++) //for each comb state considered

+	for (int iN=0; iN < *N; iN++) //for each comb state considered

 	{

-		for(int imod=0; imod < *Nmod; imod++) //for each modification of this comb state

+		res = comb_states[iN];

+		for (int imod=(*Nmod-1); imod >= 0; imod--) //for each modification of this comb state

 		{

-			binary_states[iN][imod] = comb_states[iN]&(int)pow(2,*Nmod-imod-1);//if =0->hidden state comb_states[iN] has modification imod non enriched; if !=0->enriched

-			if (binary_states[iN][imod] != 0 )

-				binary_states[iN][imod] = 1;

+			binary_states[iN][imod] = (bool)fmod(res,2);

+			res = (res - (double)binary_states[iN][imod]) / 2.0;

 		}

 	}

old mode 100755

new mode 100644

old mode 100644

new mode 100755

@@ -17,7 +17,9 @@ void univariate_hmm(int* O, int* T, int* N, double* size, double* prob, int* max

 	//FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;

 	// Parallelization settings

+	#ifdef _OPENMP

 	omp_set_num_threads(*num_threads);

+	#endif

 	// Print some information

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

@@ -267,7 +269,9 @@ void multivariate_hmm(int* O, int* T, int* N, int *Nmod, int* comb_states, doubl

 	//FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;

 	// Parallelization settings

+	#ifdef _OPENMP

 	omp_set_num_threads(*num_threads);

+	#endif

 	// Print some information

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

@@ -1,13 +1,4 @@

-#include "utility.h"

-#include "scalehmm.h"

-#include <vector> // storing density functions in multivariate

-#include <string> // strcmp

-

-#if defined TARGET_OS_MAC || defined __APPLE__

-#include <libiomp/omp.h> // parallelization options on mac

-#elif defined __linux__ || defined _WIN32 || defined _WIN64

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

-#endif

+#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 int** multiO;

@@ -15,8 +6,7 @@ static int** multiO;

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

 // 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, double* size, double* prob, int* maxiter, int* maxtime, double* eps, double* posteriors, double* densities, bool* keep_densities, 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)

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

 {

 	// Define logging level

@@ -130,15 +120,8 @@ void R_univariate_hmm(int* O, int* T, int* N, double* size, double* prob, int* m

 				// Disturb mean and variance for use as randomized initial parameters

 				//FILE_LOG(logINFO) << "Using random initialization for size and prob";

 				if (*verbosity>=1) Rprintf("HMM: Using random initialization for size and prob\n");

-				srand (clock());

-				int rand1, rand2;

-				rand1 = rand();

-				rand2 = rand();

-				imean = (double)rand1/(double)RAND_MAX * 10*mean;

-				ivariance = imean + (double)rand2/(double)RAND_MAX * 20*imean; // variance has to be greater than mean, otherwise r will be negative

-				//FILE_LOG(logDEBUG2) << "RAND_MAX = " << RAND_MAX;

-				//FILE_LOG(logDEBUG2) << "rand1 = " << rand1;

-				//FILE_LOG(logDEBUG2) << "rand2 = " << rand2;

+					imean = runif(0, 10*mean);

+					ivariance = imean + runif(0, 20*imean); // variance has to be greater than mean, otherwise r will be negative

 				//FILE_LOG(logDEBUG2) << "imean = " << imean;

 				//FILE_LOG(logDEBUG2) << "ivariance = " << ivariance;

 			}

@@ -269,20 +252,18 @@ void R_univariate_hmm(int* O, int* T, int* N, double* size, double* prob, int* m

 	delete hmm;

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

 }

-} // extern C

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

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

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

-extern "C" {

-void R_multivariate_hmm(int* O, int* T, int* N, int *Nmod, int* 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, int* states, double* A, double* proba, double* loglik, double* initial_A, double* initial_proba, bool* use_initial_params, int* num_threads, int* error, int* verbosity)

+void multivariate_hmm(int* O, int* T, int* N, int *Nmod, int* 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, int* states, double* A, double* proba, double* loglik, double* initial_A, double* initial_proba, bool* use_initial_params, int* num_threads, int* error, 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("ERROR");

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

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

 	//FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;

 	// Parallelization settings

@@ -480,25 +461,20 @@ void R_multivariate_hmm(int* O, int* T, int* N, int *Nmod, int* comb_states, dou

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

 // 	FreeIntMatrix(multiO, *Nmod); // free on.exit() in R code

 }

-} // extern C

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

 // 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 "C" {

-void R_multivariate_cleanup(int* Nmod)

+void multivariate_cleanup(int* Nmod)

 {

 	delete hmm;

 	FreeIntMatrix(multiO, *Nmod);

 }

-}

@@ -116,7 +116,7 @@ void R_univariate_hmm(int* O, int* T, int* N, double* size, double* prob, int* m

 				else if (i_state == 2)

 				{

 					//FILE_LOG(logDEBUG) << "Initializing size and prob for state 2";

-					imean = mean*2;

+					imean = mean*1.5;

 					ivariance = variance*2;

 				} 

 				// Make sure variance is greater than mean

@@ -116,7 +116,7 @@ void R_univariate_hmm(int* O, int* T, int* N, double* size, double* prob, int* m

 				else if (i_state == 2)

 				{

 					//FILE_LOG(logDEBUG) << "Initializing size and prob for state 2";

-					imean = mean+1;

+					imean = mean*2;

 					ivariance = variance*2;

 				} 

 				// Make sure variance is greater than mean

@@ -70,17 +70,24 @@ void R_univariate_hmm(int* O, int* T, int* N, double* size, double* prob, int* m

 	hmm->initialize_proba(initial_proba, *use_initial_params);

 	// Calculate mean and variance of data

-	double mean = 0, variance = 0;

+	double Tadjust = 0, mean = 0, variance = 0;

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

 	{

-		mean+= O[t];

+		if (O[t]>0)

+		{

+			mean += O[t];

+			Tadjust += 1;

+		}

 	}

-	mean = mean / *T;

+	mean = mean / Tadjust;

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

 	{

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

+		if (O[t]>0)

+		{

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

+		}

 	}

-	variance = variance / *T;

+	variance = variance / Tadjust;

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

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

@@ -110,7 +110,7 @@ void R_univariate_hmm(int* O, int* T, int* N, double* size, double* prob, int* m

 				{

 					//FILE_LOG(logDEBUG) << "Initializing size and prob for state 2";

 					imean = mean+1;

-					ivariance = variance;

+					ivariance = variance*2;

 				} 

 				// Make sure variance is greater than mean

 				if (imean >= ivariance)

@@ -20,10 +20,10 @@ void R_univariate_hmm(int* O, int* T, int* N, double* size, double* prob, int* m

 {

 	// Define logging level

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

-// 	Output2//FILE::Stream() = pFile;

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

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

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

+ 	//Output2FILE::Stream() = pFile;

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

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

 	//FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;

 	// Parallelization settings

@@ -272,10 +272,10 @@ void R_multivariate_hmm(int* O, int* T, int* N, int *Nmod, int* comb_states, dou

 {

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

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

-// 	Output2//FILE::Stream() = pFile;

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

-//  	//FILELog::ReportingLevel() = //FILELog::FromString("DEBUG3");

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

+	//Output2FILE::Stream() = pFile;

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

+ 	FILELog::ReportingLevel() = FILELog::FromString("DEBUG3");

 	//FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;

 	// Parallelization settings

@@ -304,8 +304,8 @@ void R_multivariate_hmm(int* O, int* T, int* N, int *Nmod, int* comb_states, dou

 	}

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

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

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

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

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

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

 	// Flush Rprintf statements to console

 	R_FlushConsole();

@@ -1,9 +1,14 @@

 #include "utility.h"

 #include "scalehmm.h"

 #include <vector> // storing density functions in multivariate

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

 #include <string> // strcmp

+#if defined TARGET_OS_MAC || defined __APPLE__

+#include <libiomp/omp.h> // parallelization options on mac

+#elif defined __linux__ || defined _WIN32 || defined _WIN64

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

+#endif

+

 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;

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

 #include "utility.h"

 #include "scalehmm.h"

 #include <vector> // storing density functions in multivariate

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

@@ -22,7 +22,7 @@ void R_univariate_hmm(int* O, int* T, int* N, double* size, double* prob, int* m

 	//FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;

 	// Parallelization settings

-// 	omp_set_num_threads(*num_threads);

+	omp_set_num_threads(*num_threads);

 	// Print some information

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

@@ -274,7 +274,7 @@ void R_multivariate_hmm(int* O, int* T, int* N, int *Nmod, int* comb_states, dou

 	//FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;

 	// 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,7 +1,7 @@

 #include "utility.h"

 #include "scalehmm.h"

 #include <vector> // storing density functions in multivariate

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

@@ -22,7 +22,7 @@ void R_univariate_hmm(int* O, int* T, int* N, double* size, double* prob, int* m

 	//FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;

 	// Parallelization settings

-	omp_set_num_threads(*num_threads);

+// 	omp_set_num_threads(*num_threads);

 	// Print some information

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

@@ -274,7 +274,7 @@ void R_multivariate_hmm(int* O, int* T, int* N, int *Nmod, int* comb_states, dou

 	//FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;

 	// 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,7 +1,7 @@

 #include "utility.h"

 #include "scalehmm.h"

 #include <vector> // storing density functions in multivariate

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

@@ -22,7 +22,7 @@ void R_univariate_hmm(int* O, int* T, int* N, double* size, double* prob, int* m

 	//FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;

 	// Parallelization settings

-// 	omp_set_num_threads(*num_threads);

+	omp_set_num_threads(*num_threads);

 	// Print some information

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

@@ -274,7 +274,7 @@ void R_multivariate_hmm(int* O, int* T, int* N, int *Nmod, int* comb_states, dou

 	//FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;

 	// 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,7 +1,7 @@

 #include "utility.h"

 #include "scalehmm.h"

 #include <vector> // storing density functions in multivariate

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

@@ -22,7 +22,7 @@ void R_univariate_hmm(int* O, int* T, int* N, double* size, double* prob, int* m

 	//FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;

 	// Parallelization settings

-	omp_set_num_threads(*num_threads);

+// 	omp_set_num_threads(*num_threads);

 	// Print some information

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

@@ -274,7 +274,7 @@ void R_multivariate_hmm(int* O, int* T, int* N, int *Nmod, int* comb_states, dou

 	//FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;

 	// Parallelization settings

-	omp_set_num_threads(*num_threads);

+// 	omp_set_num_threads(*num_threads);

 	// Print some information

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

@@ -434,7 +434,6 @@ void R_multivariate_hmm(int* O, int* T, int* N, int *Nmod, int* comb_states, dou

 			}

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

-// 			if (*verbosity>=1) Rprintf("HMM: comb_states[t=%d, ind_max=%d] = %d\n", t, ind_max, states[t]);

 		}

 // 	}

 // 	else

@@ -11,7 +11,7 @@ static int** multiO;

 // 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, double* size, double* prob, int* maxiter, int* maxtime, double* eps, double* posteriors, 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)

+void R_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, 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)

 {

 	// Define logging level

@@ -26,27 +26,27 @@ void R_univariate_hmm(int* O, int* T, int* N, double* size, double* prob, int* m

 	// Print some information

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

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

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

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

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

+	if (*verbosity>=1) Rprintf("HMM: 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("HMM: 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("HMM: 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("HMM: 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("HMM: maximum running time = %d sec\n", *maxtime);

 	}

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

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

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

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

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

@@ -58,7 +58,7 @@ void R_univariate_hmm(int* O, int* T, int* N, double* size, double* prob, int* m

 	// Create the HMM

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

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

+	hmm = new ScaleHMM(*T, *N, *verbosity);

 	hmm->set_cutoff(*read_cutoff);

 	// Initialize the transition probabilities and proba

 	hmm->initialize_transition_probs(initial_A, *use_initial_params);

@@ -77,7 +77,7 @@ void R_univariate_hmm(int* O, int* T, int* N, double* size, double* prob, int* m

 	}

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

 	// Go through all states of the hmm and assign the density functions

 	double imean=0, ivariance=0;

@@ -85,7 +85,7 @@ void R_univariate_hmm(int* O, int* T, int* N, double* size, double* prob, int* m

 	{

 		if (*use_initial_params) {

 			//FILE_LOG(logINFO) << "Using given parameters for size and prob";

-			Rprintf("Using given parameters for size and prob\n");

+			if (*verbosity>=1) Rprintf("HMM: Using given parameters for size and prob\n");

 			imean = (1-initial_prob[i_state])*initial_size[i_state] / initial_prob[i_state];

 			ivariance = imean / initial_prob[i_state];

 			//FILE_LOG(logDEBUG2) << "imean = " << imean;

@@ -117,7 +117,7 @@ void R_univariate_hmm(int* O, int* T, int* N, double* size, double* prob, int* m

 			{

 				// Disturb mean and variance for use as randomized initial parameters

 				//FILE_LOG(logINFO) << "Using random initialization for size and prob";

-				Rprintf("Using random initialization for size and prob\n");

+				if (*verbosity>=1) Rprintf("HMM: Using random initialization for size and prob\n");

 				srand (clock());

 				int rand1, rand2;

 				rand1 = rand();

@@ -136,14 +136,14 @@ void R_univariate_hmm(int* O, int* T, int* N, double* size, double* prob, int* m

 				if (i_state == 1)

 				{

 					//FILE_LOG(logINFO) << "Initializing r and p empirically for state 1";

-					Rprintf("Initializing r and p empirically for state 1\n");

+					if (*verbosity>=1) Rprintf("HMM: Initializing r and p empirically for state 1\n");

 					imean = mean/2;

 					ivariance = imean*2;

 				}

 				else if (i_state == 2)

 				{

 					//FILE_LOG(logINFO) << "Initializing r and p empirically for state 2";

-					Rprintf("Initializing r and p empirically for state 2\n");

+					if (*verbosity>=1) Rprintf("HMM: Initializing r and p empirically for state 2\n");

 					imean = mean*2;

 					ivariance = imean*2;

 				} 

@@ -187,15 +187,16 @@ void R_univariate_hmm(int* O, int* T, int* N, double* size, double* prob, int* m

 	catch (std::exception& e)

 	{

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

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

+		if (*verbosity>=1) Rprintf("HMM: Error in Baum-Welch: %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

+

+	// Get the posteriors and save results directly to the R pointer

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

-	Rprintf("Recoding posteriors ...\n");

+	if (*verbosity>=1) Rprintf("HMM: Recoding posteriors ...\n");

 	R_FlushConsole();

 	#pragma omp parallel for

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

@@ -206,6 +207,22 @@ void R_univariate_hmm(int* O, int* T, int* N, double* size, double* prob, int* m

 		}

 	}

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

+			}

+		}

+	}

+

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

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

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

@@ -246,7 +263,7 @@ void R_univariate_hmm(int* O, int* T, int* N, double* size, double* prob, int* m

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

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

 extern "C" {

-void R_multivariate_hmm(int* O, int* T, int* N, int *Nmod, int* 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, 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(int* O, int* T, int* N, int *Nmod, int* 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, int* states, double* A, double* proba, double* loglik, double* initial_A, double* initial_proba, bool* use_initial_params, int* num_threads, int* error, int* verbosity)

 {

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

@@ -261,29 +278,29 @@ void R_multivariate_hmm(int* O, int* T, int* N, int *Nmod, int* comb_states, dou

 	// Print some information

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

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

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

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

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

+	if (*verbosity>=1) Rprintf("HMM: 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("HMM: 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("HMM: 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("HMM: 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("HMM: maximum running time = %d sec\n", *maxtime);

 	}

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

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

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

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

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

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

 	// Flush Rprintf statements to console

 	R_FlushConsole();

@@ -303,7 +320,7 @@ void R_multivariate_hmm(int* O, int* T, int* N, int *Nmod, int* comb_states, dou

 	// Create the HMM

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

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

+	hmm = new ScaleHMM(*T, *N, *Nmod, *verbosity);

 	// Initialize the transition probabilities and proba

 	hmm->initialize_transition_probs(initial_A, *use_initial_params);

 	hmm->initialize_proba(initial_proba, *use_initial_params);

@@ -365,17 +382,17 @@ void R_multivariate_hmm(int* O, int* T, int* N, int *Nmod, int* comb_states, dou

 	catch (std::exception& e)

 	{

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

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

+		if (*verbosity>=1) Rprintf("HMM: Error in Baum-Welch: %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

+	// Get the posteriors and save results directly to the R pointer

 	if (*keep_posteriors == true)

 	{

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

-		Rprintf("Recoding posteriors ...\n");

+		if (*verbosity>=1) Rprintf("HMM: Recoding posteriors ...\n");

 		R_FlushConsole();

 		#pragma omp parallel for

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

@@ -387,6 +404,22 @@ void R_multivariate_hmm(int* O, int* T, int* N, int *Nmod, int* comb_states, dou

 		}

 	}

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

+			}

+		}

+	}

+

 	// Compute the states from posteriors

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

 // 	if (*fdr == -1)

@@ -401,7 +434,7 @@ void R_multivariate_hmm(int* O, int* T, int* N, int *Nmod, int* comb_states, dou

 			}

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

-// 			Rprintf("comb_states[t=%d, ind_max=%d] = %d\n", t, ind_max, states[t]);

+// 			if (*verbosity>=1) Rprintf("HMM: comb_states[t=%d, ind_max=%d] = %d\n", t, ind_max, states[t]);

 		}

 // 	}

 // 	else

@@ -391,6 +391,7 @@ void R_multivariate_hmm(int* O, int* T, int* N, int *Nmod, int* comb_states, dou

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

 // 	if (*fdr == -1)

 // 	{

+		int ind_max;

 		std::vector<double> posterior_per_t(*N);

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

 		{

@@ -398,7 +399,9 @@ void R_multivariate_hmm(int* O, int* T, int* N, int *Nmod, int* comb_states, dou

 			{

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

 			}

-			states[t] = comb_states[(int)*std::max_element(posterior_per_t.begin(), posterior_per_t.end())];

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

+// 			Rprintf("comb_states[t=%d, ind_max=%d] = %d\n", t, ind_max, states[t]);

 		}

 // 	}

 // 	else

@@ -15,49 +15,49 @@ void R_univariate_hmm(int* O, int* T, int* N, double* size, double* prob, int* m

 {

 	// Define logging level

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

-// 	Output2FILE::Stream() = pFile;

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

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

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

+// 	Output2//FILE::Stream() = pFile;

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

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

-	FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;