Browse code

experiment.table can contain multiple of the same files in column "file", save "state" column as double instead of integer

chakalakka authored on 18/05/2017 11:49:42
Showing1 changed files
... ...
@@ -273,7 +273,7 @@ void univariate_hmm(int* O, int* T, int* N, double* size, double* prob, int* max
273 273
 // =====================================================================================================================================================
274 274
 // This function takes parameters from R, creates a multivariate HMM object, creates the distributions, runs the Baum-Welch and returns the result to R.
275 275
 // =====================================================================================================================================================
276
-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)
276
+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)
277 277
 {
278 278
 
279 279
 	// Define logging level {"ERROR", "WARNING", "INFO", "ITERATION", "DEBUG", "DEBUG1", "DEBUG2", "DEBUG3", "DEBUG4"}
Browse code

more memory efficient stepsize calculation

chakalakka authored on 03/05/2017 15:39:22
Showing1 changed files
... ...
@@ -6,7 +6,7 @@ static int** multiO;
6 6
 // ===================================================================================================================================================
7 7
 // This function takes parameters from R, creates a univariate HMM object, creates the distributions, runs the Baum-Welch and returns the result to R.
8 8
 // ===================================================================================================================================================
9
-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)
9
+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)
10 10
 {
11 11
 
12 12
 	// Define logging level
... ...
@@ -220,6 +220,21 @@ void univariate_hmm(int* O, int* T, int* N, double* size, double* prob, int* max
220 220
 		}
221 221
 	}
222 222
 
223
+	// Compute the states from posteriors
224
+	//FILE_LOG(logDEBUG1) << "Computing states from posteriors";
225
+	int ind_max;
226
+	std::vector<double> posterior_per_t(*N);
227
+	for (int t=0; t<*T; t++)
228
+	{
229
+		for (int iN=0; iN<*N; iN++)
230
+		{
231
+			posterior_per_t[iN] = hmm->get_posterior(iN, t);
232
+		}
233
+		ind_max = std::distance(posterior_per_t.begin(), std::max_element(posterior_per_t.begin(), posterior_per_t.end()));
234
+		states[t] = ind_max + 1;
235
+		maxPosterior[t] = posterior_per_t[ind_max];
236
+	}
237
+		
223 238
 	//FILE_LOG(logDEBUG1) << "Return parameters";
224 239
 	// also return the estimated transition matrix and the initial probs
225 240
 	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
258 273
 // =====================================================================================================================================================
259 274
 // This function takes parameters from R, creates a multivariate HMM object, creates the distributions, runs the Baum-Welch and returns the result to R.
260 275
 // =====================================================================================================================================================
261
-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)
276
+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)
262 277
 {
263 278
 
264 279
 	// 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
432 447
 			}
433 448
 			ind_max = std::distance(posterior_per_t.begin(), std::max_element(posterior_per_t.begin(), posterior_per_t.end()));
434 449
 			states[t] = comb_states[ind_max];
450
+			maxPosterior[t] = posterior_per_t[ind_max];
435 451
 		}
436 452
 // 	}
437 453
 // 	else
... ...
@@ -506,46 +522,22 @@ void array3D_which_max(double* array3D, int* dim, int* ind_max)
506 522
 	
507 523
 }
508 524
 
509
-
510
-// ====================================================================================
511
-// C version of apply(array2D, MARGIN = 1, FUN = mean)
512
-// ====================================================================================
513
-void array2D_mean(double* array2D, int* dim, double* mean)
525
+// ====================================================================
526
+// C version of apply(array2D, 1, which.max) and apply(array2D, 1, max)
527
+// ====================================================================
528
+void array2D_which_max(double* array2D, int* dim, int* ind_max, double* value_max)
514 529
 {
515
-  // array2D is actually a vector, but is intended to originate from a matrix in R
516
-  double sum=0;
530
+  // array2D is actually a vector, but is intended to originate from a 2D array in R
531
+	std::vector<double> value_per_i0(dim[1]);
517 532
   for (int i0=0; i0<dim[0]; i0++)
518 533
   {
519
-    sum = 0;
520 534
     for (int i1=0; i1<dim[1]; i1++)
521 535
     {
522
-      sum += array2D[i1*dim[0] + i0];
523
-    }
524
-    mean[i0] = sum / dim[1];
525
-  }
526
-	
527
-}
528
-
529
-
530
-// ====================================================================================
531
-// C version of apply(array3D, MARGIN = c(1,2), FUN = mean)
532
-// ====================================================================================
533
-void array3D_mean(double* array3D, int* dim, double* mean)
534
-{
535
-  // array3D is actually a vector, but is intended to originate from a 3D array in R
536
-  double sum=0;
537
-  for (int i0=0; i0<dim[0]; i0++)
538
-  {
539
-    for (int i1=0; i1<dim[1]; i1++)
540
-    {
541
-      sum = 0;
542
-      for (int i2=0; i2<dim[2]; i2++)
543
-      {
544
-        sum += array3D[(i2*dim[1] + i1)*dim[0] + i0];
545
-        // 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]);
546
-      }
547
-      mean[i1*dim[0] + i0] = sum / dim[2];
536
+			value_per_i0[i1] = array2D[i1 * dim[0] + i0];
537
+      // Rprintf("i0=%d, i1=%d, value_per_i0[%d] = %g\n", i0, i1, i1, value_per_i0[i1]);
548 538
     }
539
+		ind_max[i0] = 1 + std::distance(value_per_i0.begin(), std::max_element(value_per_i0.begin(), value_per_i0.end()));
540
+    value_max[i0] = *std::max_element(value_per_i0.begin(), value_per_i0.end());
549 541
   }
550 542
 	
551 543
 }
552 544
\ No newline at end of file
Browse code

C helper functions for time-intensive calculations in stepsize posterior selection

chakalakka authored on 09/04/2017 12:54:46
Showing1 changed files
... ...
@@ -483,3 +483,69 @@ void multivariate_cleanup(int* Nmod)
483 483
 	FreeIntMatrix(multiO, *Nmod);
484 484
 }
485 485
 
486
+
487
+// =======================================================
488
+// C version of apply(array3D, 1, which.max)
489
+// =======================================================
490
+void array3D_which_max(double* array3D, int* dim, int* ind_max)
491
+{
492
+  // array3D is actually a vector, but is intended to originate from a 3D array in R
493
+	std::vector<double> value_per_i0(dim[1] * dim[2]);
494
+  for (int i0=0; i0<dim[0]; i0++)
495
+  {
496
+    for (int i1=0; i1<dim[1]; i1++)
497
+    {
498
+      for (int i2=0; i2<dim[2]; i2++)
499
+      {
500
+  			value_per_i0[i1*dim[2]+i2] = array3D[(i1*dim[2]+i2) * dim[0] + i0];
501
+        // 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]);
502
+      }
503
+    }
504
+		ind_max[i0] = 1 + std::distance(value_per_i0.begin(), std::max_element(value_per_i0.begin(), value_per_i0.end()));
505
+  }
506
+	
507
+}
508
+
509
+
510
+// ====================================================================================
511
+// C version of apply(array2D, MARGIN = 1, FUN = mean)
512
+// ====================================================================================
513
+void array2D_mean(double* array2D, int* dim, double* mean)
514
+{
515
+  // array2D is actually a vector, but is intended to originate from a matrix in R
516
+  double sum=0;
517
+  for (int i0=0; i0<dim[0]; i0++)
518
+  {
519
+    sum = 0;
520
+    for (int i1=0; i1<dim[1]; i1++)
521
+    {
522
+      sum += array2D[i1*dim[0] + i0];
523
+    }
524
+    mean[i0] = sum / dim[1];
525
+  }
526
+	
527
+}
528
+
529
+
530
+// ====================================================================================
531
+// C version of apply(array3D, MARGIN = c(1,2), FUN = mean)
532
+// ====================================================================================
533
+void array3D_mean(double* array3D, int* dim, double* mean)
534
+{
535
+  // array3D is actually a vector, but is intended to originate from a 3D array in R
536
+  double sum=0;
537
+  for (int i0=0; i0<dim[0]; i0++)
538
+  {
539
+    for (int i1=0; i1<dim[1]; i1++)
540
+    {
541
+      sum = 0;
542
+      for (int i2=0; i2<dim[2]; i2++)
543
+      {
544
+        sum += array3D[(i2*dim[1] + i1)*dim[0] + i0];
545
+        // 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]);
546
+      }
547
+      mean[i1*dim[0] + i0] = sum / dim[2];
548
+    }
549
+  }
550
+	
551
+}
486 552
\ No newline at end of file
Browse code

added support for up to 53 experiments in multivariate

chakalakka authored on 09/06/2016 13:55:17
Showing1 changed files
... ...
@@ -258,7 +258,7 @@ void univariate_hmm(int* O, int* T, int* N, double* size, double* prob, int* max
258 258
 // =====================================================================================================================================================
259 259
 // This function takes parameters from R, creates a multivariate HMM object, creates the distributions, runs the Baum-Welch and returns the result to R.
260 260
 // =====================================================================================================================================================
261
-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)
261
+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)
262 262
 {
263 263
 
264 264
 	// 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
334 334
 
335 335
 	// 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)
336 336
 	//FILE_LOG(logDEBUG1) << "Preparing the binary_states vector";
337
+	double res;
337 338
 	bool **binary_states = CallocBoolMatrix(*N, *Nmod);
338
-	for(int iN=0; iN < *N; iN++) //for each comb state considered
339
+	for (int iN=0; iN < *N; iN++) //for each comb state considered
339 340
 	{
340
-		for(int imod=0; imod < *Nmod; imod++) //for each modification of this comb state
341
+		res = comb_states[iN];
342
+		for (int imod=(*Nmod-1); imod >= 0; imod--) //for each modification of this comb state
341 343
 		{
342
-			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
343
-			if (binary_states[iN][imod] != 0 )
344
-				binary_states[iN][imod] = 1;
344
+			binary_states[iN][imod] = (bool)fmod(res,2);
345
+			res = (res - (double)binary_states[iN][imod]) / 2.0;
345 346
 		}
346 347
 	}
347 348
 
Browse code

corrected file permissions

chakalakka authored on 08/06/2016 09:00:02
Showing1 changed files
1 1
old mode 100755
2 2
new mode 100644
Browse code

fixed seqlengths assign bug, ubuntu reinstall

chakalakka authored on 07/06/2016 09:48:47
Showing1 changed files
1 1
old mode 100644
2 2
new mode 100755
Browse code

compiles on OSX

chakalakka authored on 16/04/2016 08:02:18
Showing1 changed files
... ...
@@ -17,7 +17,9 @@ void univariate_hmm(int* O, int* T, int* N, double* size, double* prob, int* max
17 17
 
18 18
 	//FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;
19 19
 	// Parallelization settings
20
+	#ifdef _OPENMP
20 21
 	omp_set_num_threads(*num_threads);
22
+	#endif
21 23
 
22 24
 	// Print some information
23 25
 	//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
267 269
 
268 270
 	//FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;
269 271
 	// Parallelization settings
272
+	#ifdef _OPENMP
270 273
 	omp_set_num_threads(*num_threads);
274
+	#endif
271 275
 
272 276
 	// Print some information
273 277
 	//FILE_LOG(logINFO) << "number of states = " << *N;
Browse code

on the way to BiocCheck

chakalakka authored on 30/03/2016 12:09:31
Showing1 changed files
... ...
@@ -1,13 +1,4 @@
1
-#include "utility.h"
2
-#include "scalehmm.h"
3
-#include <vector> // storing density functions in multivariate
4
-#include <string> // strcmp
5
-
6
-#if defined TARGET_OS_MAC || defined __APPLE__
7
-#include <libiomp/omp.h> // parallelization options on mac
8
-#elif defined __linux__ || defined _WIN32 || defined _WIN64
9
-#include <omp.h> // parallelization options
10
-#endif
1
+#include "R_interface.h"
11 2
 
12 3
 static ScaleHMM* hmm; // declare as static outside the function because we only need one and this enables memory-cleanup on R_CheckUserInterrupt()
13 4
 static int** multiO;
... ...
@@ -15,8 +6,7 @@ static int** multiO;
15 6
 // ===================================================================================================================================================
16 7
 // This function takes parameters from R, creates a univariate HMM object, creates the distributions, runs the Baum-Welch and returns the result to R.
17 8
 // ===================================================================================================================================================
18
-extern "C" {
19
-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)
9
+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)
20 10
 {
21 11
 
22 12
 	// Define logging level
... ...
@@ -130,15 +120,8 @@ void R_univariate_hmm(int* O, int* T, int* N, double* size, double* prob, int* m
130 120
 				// Disturb mean and variance for use as randomized initial parameters
131 121
 				//FILE_LOG(logINFO) << "Using random initialization for size and prob";
132 122
 				if (*verbosity>=1) Rprintf("HMM: Using random initialization for size and prob\n");
133
-				srand (clock());
134
-				int rand1, rand2;
135
-				rand1 = rand();
136
-				rand2 = rand();
137
-				imean = (double)rand1/(double)RAND_MAX * 10*mean;
138
-				ivariance = imean + (double)rand2/(double)RAND_MAX * 20*imean; // variance has to be greater than mean, otherwise r will be negative
139
-				//FILE_LOG(logDEBUG2) << "RAND_MAX = " << RAND_MAX;
140
-				//FILE_LOG(logDEBUG2) << "rand1 = " << rand1;
141
-				//FILE_LOG(logDEBUG2) << "rand2 = " << rand2;
123
+					imean = runif(0, 10*mean);
124
+					ivariance = imean + runif(0, 20*imean); // variance has to be greater than mean, otherwise r will be negative
142 125
 				//FILE_LOG(logDEBUG2) << "imean = " << imean;
143 126
 				//FILE_LOG(logDEBUG2) << "ivariance = " << ivariance;
144 127
 			}
... ...
@@ -269,20 +252,18 @@ void R_univariate_hmm(int* O, int* T, int* N, double* size, double* prob, int* m
269 252
 	delete hmm;
270 253
 	hmm = NULL; // assign NULL to defuse the additional delete in on.exit() call
271 254
 }
272
-} // extern C
273 255
 
274 256
 // =====================================================================================================================================================
275 257
 // This function takes parameters from R, creates a multivariate HMM object, creates the distributions, runs the Baum-Welch and returns the result to R.
276 258
 // =====================================================================================================================================================
277
-extern "C" {
278
-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)
259
+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)
279 260
 {
280 261
 
281 262
 	// Define logging level {"ERROR", "WARNING", "INFO", "ITERATION", "DEBUG", "DEBUG1", "DEBUG2", "DEBUG3", "DEBUG4"}
282 263
  	//FILE* pFile = fopen("chromStar.log", "w");
283 264
 	//Output2FILE::Stream() = pFile;
284 265
  	//FILELog::ReportingLevel() = FILELog::FromString("ERROR");
285
- 	FILELog::ReportingLevel() = FILELog::FromString("DEBUG3");
266
+ 	//FILELog::ReportingLevel() = FILELog::FromString("DEBUG3");
286 267
 
287 268
 	//FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;
288 269
 	// Parallelization settings
... ...
@@ -480,25 +461,20 @@ void R_multivariate_hmm(int* O, int* T, int* N, int *Nmod, int* comb_states, dou
480 461
 	hmm = NULL; // assign NULL to defuse the additional delete in on.exit() call
481 462
 // 	FreeIntMatrix(multiO, *Nmod); // free on.exit() in R code
482 463
 }
483
-} // extern C
484 464
 
485 465
 
486 466
 // =======================================================
487 467
 // This function make a cleanup if anything was left over
488 468
 // =======================================================
489
-extern "C" {
490
-void R_univariate_cleanup()
469
+void univariate_cleanup()
491 470
 {
492 471
 // 	//FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__; // This message will be shown if interrupt happens before start of C-code
493 472
 	delete hmm;
494 473
 }
495
-}
496 474
 
497
-extern "C" {
498
-void R_multivariate_cleanup(int* Nmod)
475
+void multivariate_cleanup(int* Nmod)
499 476
 {
500 477
 	delete hmm;
501 478
 	FreeIntMatrix(multiO, *Nmod);
502 479
 }
503
-}
504 480
 
Browse code

test change standard init

chakalakka authored on 11/02/2016 11:15:29
Showing1 changed files
... ...
@@ -116,7 +116,7 @@ void R_univariate_hmm(int* O, int* T, int* N, double* size, double* prob, int* m
116 116
 				else if (i_state == 2)
117 117
 				{
118 118
 					//FILE_LOG(logDEBUG) << "Initializing size and prob for state 2";
119
-					imean = mean*2;
119
+					imean = mean*1.5;
120 120
 					ivariance = variance*2;
121 121
 				} 
122 122
 				// Make sure variance is greater than mean
Browse code

changed standard init again

chakalakka authored on 11/02/2016 09:28:05
Showing1 changed files
... ...
@@ -116,7 +116,7 @@ void R_univariate_hmm(int* O, int* T, int* N, double* size, double* prob, int* m
116 116
 				else if (i_state == 2)
117 117
 				{
118 118
 					//FILE_LOG(logDEBUG) << "Initializing size and prob for state 2";
119
-					imean = mean+1;
119
+					imean = mean*2;
120 120
 					ivariance = variance*2;
121 121
 				} 
122 122
 				// Make sure variance is greater than mean
Browse code

added pre-fitting capabilities to univariate

chakalakka authored on 10/02/2016 18:32:59
Showing1 changed files
... ...
@@ -70,17 +70,24 @@ void R_univariate_hmm(int* O, int* T, int* N, double* size, double* prob, int* m
70 70
 	hmm->initialize_proba(initial_proba, *use_initial_params);
71 71
     
72 72
 	// Calculate mean and variance of data
73
-	double mean = 0, variance = 0;
73
+	double Tadjust = 0, mean = 0, variance = 0;
74 74
 	for(int t=0; t<*T; t++)
75 75
 	{
76
-		mean+= O[t];
76
+		if (O[t]>0)
77
+		{
78
+			mean += O[t];
79
+			Tadjust += 1;
80
+		}
77 81
 	}
78
-	mean = mean / *T;
82
+	mean = mean / Tadjust;
79 83
 	for(int t=0; t<*T; t++)
80 84
 	{
81
-		variance+= pow(O[t] - mean, 2);
85
+		if (O[t]>0)
86
+		{
87
+			variance += pow(O[t] - mean, 2);
88
+		}
82 89
 	}
83
-	variance = variance / *T;
90
+	variance = variance / Tadjust;
84 91
 	//FILE_LOG(logINFO) << "data mean = " << mean << ", data variance = " << variance;		
85 92
 	if (*verbosity>=1) Rprintf("HMM: data mean = %g, data variance = %g\n", mean, variance);		
86 93
 	
Browse code

changed standard initialization procedure

chakalakka authored on 10/02/2016 17:29:44
Showing1 changed files
... ...
@@ -110,7 +110,7 @@ void R_univariate_hmm(int* O, int* T, int* N, double* size, double* prob, int* m
110 110
 				{
111 111
 					//FILE_LOG(logDEBUG) << "Initializing size and prob for state 2";
112 112
 					imean = mean+1;
113
-					ivariance = variance;
113
+					ivariance = variance*2;
114 114
 				} 
115 115
 				// Make sure variance is greater than mean
116 116
 				if (imean >= ivariance)
Browse code

changed all densities=0 behaviour

chakalakka authored on 06/10/2015 14:19:40
Showing1 changed files
... ...
@@ -20,10 +20,10 @@ void R_univariate_hmm(int* O, int* T, int* N, double* size, double* prob, int* m
20 20
 {
21 21
 
22 22
 	// Define logging level
23
-// 	//FILE* pFile = fopen("chromStar.log", "w");
24
-// 	Output2//FILE::Stream() = pFile;
25
- 	//FILELog::ReportingLevel() = //FILELog::FromString("ERROR");
26
-//  	//FILELog::ReportingLevel() = //FILELog::FromString("DEBUG2");
23
+	//FILE* pFile = fopen("chromStar.log", "w");
24
+ 	//Output2FILE::Stream() = pFile;
25
+ 	//FILELog::ReportingLevel() = FILELog::FromString("ERROR");
26
+ 	//FILELog::ReportingLevel() = FILELog::FromString("DEBUG2");
27 27
 
28 28
 	//FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;
29 29
 	// Parallelization settings
... ...
@@ -272,10 +272,10 @@ void R_multivariate_hmm(int* O, int* T, int* N, int *Nmod, int* comb_states, dou
272 272
 {
273 273
 
274 274
 	// Define logging level {"ERROR", "WARNING", "INFO", "ITERATION", "DEBUG", "DEBUG1", "DEBUG2", "DEBUG3", "DEBUG4"}
275
-// 	//FILE* pFile = fopen("chromStar.log", "w");
276
-// 	Output2//FILE::Stream() = pFile;
277
- 	//FILELog::ReportingLevel() = //FILELog::FromString("ERROR");
278
-//  	//FILELog::ReportingLevel() = //FILELog::FromString("DEBUG3");
275
+ 	//FILE* pFile = fopen("chromStar.log", "w");
276
+	//Output2FILE::Stream() = pFile;
277
+ 	//FILELog::ReportingLevel() = FILELog::FromString("ERROR");
278
+ 	FILELog::ReportingLevel() = FILELog::FromString("DEBUG3");
279 279
 
280 280
 	//FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;
281 281
 	// Parallelization settings
Browse code

changed print text

chakalakka authored on 24/09/2015 12:37:23
Showing1 changed files
... ...
@@ -304,8 +304,8 @@ void R_multivariate_hmm(int* O, int* T, int* N, int *Nmod, int* comb_states, dou
304 304
 	}
305 305
 	//FILE_LOG(logINFO) << "epsilon = " << *eps;
306 306
 	if (*verbosity>=1) Rprintf("HMM: epsilon = %g\n", *eps);
307
-	//FILE_LOG(logINFO) << "number of modifications = " << *Nmod;
308
-	if (*verbosity>=1) Rprintf("HMM: number of modifications = %d\n", *Nmod);
307
+	//FILE_LOG(logINFO) << "number of experiments = " << *Nmod;
308
+	if (*verbosity>=1) Rprintf("HMM: number of experiments = %d\n", *Nmod);
309 309
 
310 310
 	// Flush Rprintf statements to console
311 311
 	R_FlushConsole();
Browse code

added parallelization for mac

chakalakka authored on 09/06/2015 08:30:26
Showing1 changed files
... ...
@@ -1,9 +1,14 @@
1 1
 #include "utility.h"
2 2
 #include "scalehmm.h"
3 3
 #include <vector> // storing density functions in multivariate
4
-#include <omp.h> // parallelization options
5 4
 #include <string> // strcmp
6 5
 
6
+#if defined TARGET_OS_MAC || defined __APPLE__
7
+#include <libiomp/omp.h> // parallelization options on mac
8
+#elif defined __linux__ || defined _WIN32 || defined _WIN64
9
+#include <omp.h> // parallelization options
10
+#endif
11
+
7 12
 static ScaleHMM* hmm; // declare as static outside the function because we only need one and this enables memory-cleanup on R_CheckUserInterrupt()
8 13
 static int** multiO;
9 14
 
Browse code

readded parallelization

chakalakka authored on 04/06/2015 17:12:09
Showing1 changed files
... ...
@@ -1,7 +1,7 @@
1 1
 #include "utility.h"
2 2
 #include "scalehmm.h"
3 3
 #include <vector> // storing density functions in multivariate
4
-// #include <omp.h> // parallelization options
4
+#include <omp.h> // parallelization options
5 5
 #include <string> // strcmp
6 6
 
7 7
 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
22 22
 
23 23
 	//FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;
24 24
 	// Parallelization settings
25
-// 	omp_set_num_threads(*num_threads);
25
+	omp_set_num_threads(*num_threads);
26 26
 
27 27
 	// Print some information
28 28
 	//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
274 274
 
275 275
 	//FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;
276 276
 	// Parallelization settings
277
-// 	omp_set_num_threads(*num_threads);
277
+	omp_set_num_threads(*num_threads);
278 278
 
279 279
 	// Print some information
280 280
 	//FILE_LOG(logINFO) << "number of states = " << *N;
Browse code

no parallelization

chakalakka authored on 04/06/2015 16:04:39
Showing1 changed files
... ...
@@ -1,7 +1,7 @@
1 1
 #include "utility.h"
2 2
 #include "scalehmm.h"
3 3
 #include <vector> // storing density functions in multivariate
4
-#include <omp.h> // parallelization options
4
+// #include <omp.h> // parallelization options
5 5
 #include <string> // strcmp
6 6
 
7 7
 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
22 22
 
23 23
 	//FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;
24 24
 	// Parallelization settings
25
-	omp_set_num_threads(*num_threads);
25
+// 	omp_set_num_threads(*num_threads);
26 26
 
27 27
 	// Print some information
28 28
 	//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
274 274
 
275 275
 	//FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;
276 276
 	// Parallelization settings
277
-	omp_set_num_threads(*num_threads);
277
+// 	omp_set_num_threads(*num_threads);
278 278
 
279 279
 	// Print some information
280 280
 	//FILE_LOG(logINFO) << "number of states = " << *N;
Browse code

reinclude parallelization

chakalakka authored on 01/06/2015 12:41:50
Showing1 changed files
... ...
@@ -1,7 +1,7 @@
1 1
 #include "utility.h"
2 2
 #include "scalehmm.h"
3 3
 #include <vector> // storing density functions in multivariate
4
-// #include <omp.h> // parallelization options
4
+#include <omp.h> // parallelization options
5 5
 #include <string> // strcmp
6 6
 
7 7
 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
22 22
 
23 23
 	//FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;
24 24
 	// Parallelization settings
25
-// 	omp_set_num_threads(*num_threads);
25
+	omp_set_num_threads(*num_threads);
26 26
 
27 27
 	// Print some information
28 28
 	//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
274 274
 
275 275
 	//FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;
276 276
 	// Parallelization settings
277
-// 	omp_set_num_threads(*num_threads);
277
+	omp_set_num_threads(*num_threads);
278 278
 
279 279
 	// Print some information
280 280
 	//FILE_LOG(logINFO) << "number of states = " << *N;
Browse code

no parallelization

chakalakka authored on 01/06/2015 11:26:56
Showing1 changed files
... ...
@@ -1,7 +1,7 @@
1 1
 #include "utility.h"
2 2
 #include "scalehmm.h"
3 3
 #include <vector> // storing density functions in multivariate
4
-#include <omp.h> // parallelization options
4
+// #include <omp.h> // parallelization options
5 5
 #include <string> // strcmp
6 6
 
7 7
 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
22 22
 
23 23
 	//FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;
24 24
 	// Parallelization settings
25
-	omp_set_num_threads(*num_threads);
25
+// 	omp_set_num_threads(*num_threads);
26 26
 
27 27
 	// Print some information
28 28
 	//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
274 274
 
275 275
 	//FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;
276 276
 	// Parallelization settings
277
-	omp_set_num_threads(*num_threads);
277
+// 	omp_set_num_threads(*num_threads);
278 278
 
279 279
 	// Print some information
280 280
 	//FILE_LOG(logINFO) << "number of states = " << *N;
Browse code

improved binning

chakalakka authored on 17/02/2015 11:54:42
Showing1 changed files
... ...
@@ -434,7 +434,6 @@ void R_multivariate_hmm(int* O, int* T, int* N, int *Nmod, int* comb_states, dou
434 434
 			}
435 435
 			ind_max = std::distance(posterior_per_t.begin(), std::max_element(posterior_per_t.begin(), posterior_per_t.end()));
436 436
 			states[t] = comb_states[ind_max];
437
-// 			if (*verbosity>=1) Rprintf("HMM: comb_states[t=%d, ind_max=%d] = %d\n", t, ind_max, states[t]);
438 437
 		}
439 438
 // 	}
440 439
 // 	else
Browse code

added verbosity option

chakalakka authored on 13/02/2015 14:47:43
Showing1 changed files
... ...
@@ -11,7 +11,7 @@ static int** multiO;
11 11
 // This function takes parameters from R, creates a univariate HMM object, creates the distributions, runs the Baum-Welch and returns the result to R.
12 12
 // ===================================================================================================================================================
13 13
 extern "C" {
14
-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)
14
+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)
15 15
 {
16 16
 
17 17
 	// Define logging level
... ...
@@ -26,27 +26,27 @@ void R_univariate_hmm(int* O, int* T, int* N, double* size, double* prob, int* m
26 26
 
27 27
 	// Print some information
28 28
 	//FILE_LOG(logINFO) << "number of states = " << *N;
29
-	Rprintf("number of states = %d\n", *N);
29
+	if (*verbosity>=1) Rprintf("HMM: number of states = %d\n", *N);
30 30
 	//FILE_LOG(logINFO) << "number of bins = " << *T;
31
-	Rprintf("number of bins = %d\n", *T);
31
+	if (*verbosity>=1) Rprintf("HMM: number of bins = %d\n", *T);
32 32
 	if (*maxiter < 0)
33 33
 	{
34 34
 		//FILE_LOG(logINFO) << "maximum number of iterations = none";
35
-		Rprintf("maximum number of iterations = none\n");
35
+		if (*verbosity>=1) Rprintf("HMM: maximum number of iterations = none\n");
36 36
 	} else {
37 37
 		//FILE_LOG(logINFO) << "maximum number of iterations = " << *maxiter;
38
-		Rprintf("maximum number of iterations = %d\n", *maxiter);
38
+		if (*verbosity>=1) Rprintf("HMM: maximum number of iterations = %d\n", *maxiter);
39 39
 	}
40 40
 	if (*maxtime < 0)
41 41
 	{
42 42
 		//FILE_LOG(logINFO) << "maximum running time = none";
43
-		Rprintf("maximum running time = none\n");
43
+		if (*verbosity>=1) Rprintf("HMM: maximum running time = none\n");
44 44
 	} else {
45 45
 		//FILE_LOG(logINFO) << "maximum running time = " << *maxtime << " sec";
46
-		Rprintf("maximum running time = %d sec\n", *maxtime);
46
+		if (*verbosity>=1) Rprintf("HMM: maximum running time = %d sec\n", *maxtime);
47 47
 	}
48 48
 	//FILE_LOG(logINFO) << "epsilon = " << *eps;
49
-	Rprintf("epsilon = %g\n", *eps);
49
+	if (*verbosity>=1) Rprintf("HMM: epsilon = %g\n", *eps);
50 50
 
51 51
 	//FILE_LOG(logDEBUG3) << "observation vector";
52 52
 	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
58 58
 
59 59
 	// Create the HMM
60 60
 	//FILE_LOG(logDEBUG1) << "Creating a univariate HMM";
61
-	hmm = new ScaleHMM(*T, *N);
61
+	hmm = new ScaleHMM(*T, *N, *verbosity);
62 62
 	hmm->set_cutoff(*read_cutoff);
63 63
 	// Initialize the transition probabilities and proba
64 64
 	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
77 77
 	}
78 78
 	variance = variance / *T;
79 79
 	//FILE_LOG(logINFO) << "data mean = " << mean << ", data variance = " << variance;		
80
-	Rprintf("data mean = %g, data variance = %g\n", mean, variance);		
80
+	if (*verbosity>=1) Rprintf("HMM: data mean = %g, data variance = %g\n", mean, variance);		
81 81
 	
82 82
 	// Go through all states of the hmm and assign the density functions
83 83
 	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
85 85
 	{
86 86
 		if (*use_initial_params) {
87 87
 			//FILE_LOG(logINFO) << "Using given parameters for size and prob";
88
-			Rprintf("Using given parameters for size and prob\n");
88
+			if (*verbosity>=1) Rprintf("HMM: Using given parameters for size and prob\n");
89 89
 			imean = (1-initial_prob[i_state])*initial_size[i_state] / initial_prob[i_state];
90 90
 			ivariance = imean / initial_prob[i_state];
91 91
 			//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
117 117
 			{
118 118
 				// Disturb mean and variance for use as randomized initial parameters
119 119
 				//FILE_LOG(logINFO) << "Using random initialization for size and prob";
120
-				Rprintf("Using random initialization for size and prob\n");
120
+				if (*verbosity>=1) Rprintf("HMM: Using random initialization for size and prob\n");
121 121
 				srand (clock());
122 122
 				int rand1, rand2;
123 123
 				rand1 = rand();
... ...
@@ -136,14 +136,14 @@ void R_univariate_hmm(int* O, int* T, int* N, double* size, double* prob, int* m
136 136
 				if (i_state == 1)
137 137
 				{
138 138
 					//FILE_LOG(logINFO) << "Initializing r and p empirically for state 1";
139
-					Rprintf("Initializing r and p empirically for state 1\n");
139
+					if (*verbosity>=1) Rprintf("HMM: Initializing r and p empirically for state 1\n");
140 140
 					imean = mean/2;
141 141
 					ivariance = imean*2;
142 142
 				}
143 143
 				else if (i_state == 2)
144 144
 				{
145 145
 					//FILE_LOG(logINFO) << "Initializing r and p empirically for state 2";
146
-					Rprintf("Initializing r and p empirically for state 2\n");
146
+					if (*verbosity>=1) Rprintf("HMM: Initializing r and p empirically for state 2\n");
147 147
 					imean = mean*2;
148 148
 					ivariance = imean*2;
149 149
 				} 
... ...
@@ -187,15 +187,16 @@ void R_univariate_hmm(int* O, int* T, int* N, double* size, double* prob, int* m
187 187
 	catch (std::exception& e)
188 188
 	{
189 189
 		//FILE_LOG(logERROR) << "Error in Baum-Welch: " << e.what();
190
-		Rprintf("Error in Baum-Welch: %s\n", e.what());
190
+		if (*verbosity>=1) Rprintf("HMM: Error in Baum-Welch: %s\n", e.what());
191 191
 		if (strcmp(e.what(),"nan detected")==0) { *error = 1; }
192 192
 		else { *error = 2; }
193 193
 	}
194 194
 		
195 195
 	//FILE_LOG(logDEBUG1) << "Finished with Baum-Welch estimation";
196
-	// Compute the posteriors and save results directly to the R pointer
196
+
197
+	// Get the posteriors and save results directly to the R pointer
197 198
 	//FILE_LOG(logDEBUG1) << "Recode posteriors into column representation";
198
-	Rprintf("Recoding posteriors ...\n");
199
+	if (*verbosity>=1) Rprintf("HMM: Recoding posteriors ...\n");
199 200
 	R_FlushConsole();
200 201
 	#pragma omp parallel for
201 202
 	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
206 207
 		}
207 208
 	}
208 209
 
210
+	// Get the densities and save results directly to the R pointer
211
+	if (*keep_densities == true)
212
+	{
213
+		//FILE_LOG(logDEBUG1) << "Recode posteriors into column representation";
214
+		if (*verbosity>=1) Rprintf("HMM: Recoding densities ...\n");
215
+		R_FlushConsole();
216
+		#pragma omp parallel for
217
+		for (int iN=0; iN<*N; iN++)
218
+		{
219
+			for (int t=0; t<*T; t++)
220
+			{
221
+				densities[t + iN * (*T)] = hmm->get_density(iN, t);
222
+			}
223
+		}
224
+	}
225
+
209 226
 	//FILE_LOG(logDEBUG1) << "Return parameters";
210 227
 	// also return the estimated transition matrix and the initial probs
211 228
 	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
246 263
 // This function takes parameters from R, creates a multivariate HMM object, creates the distributions, runs the Baum-Welch and returns the result to R.
247 264
 // =====================================================================================================================================================
248 265
 extern "C" {
249
-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)
266
+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)
250 267
 {
251 268
 
252 269
 	// 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
261 278
 
262 279
 	// Print some information
263 280
 	//FILE_LOG(logINFO) << "number of states = " << *N;
264
-	Rprintf("number of states = %d\n", *N);
281
+	if (*verbosity>=1) Rprintf("HMM: number of states = %d\n", *N);
265 282
 	//FILE_LOG(logINFO) << "number of bins = " << *T;
266
-	Rprintf("number of bins = %d\n", *T);
283
+	if (*verbosity>=1) Rprintf("HMM: number of bins = %d\n", *T);
267 284
 	if (*maxiter < 0)
268 285
 	{
269 286
 		//FILE_LOG(logINFO) << "maximum number of iterations = none";
270
-		Rprintf("maximum number of iterations = none\n");
287
+		if (*verbosity>=1) Rprintf("HMM: maximum number of iterations = none\n");
271 288
 	} else {
272 289
 		//FILE_LOG(logINFO) << "maximum number of iterations = " << *maxiter;
273
-		Rprintf("maximum number of iterations = %d\n", *maxiter);
290
+		if (*verbosity>=1) Rprintf("HMM: maximum number of iterations = %d\n", *maxiter);
274 291
 	}
275 292
 	if (*maxtime < 0)
276 293
 	{
277 294
 		//FILE_LOG(logINFO) << "maximum running time = none";
278
-		Rprintf("maximum running time = none\n");
295
+		if (*verbosity>=1) Rprintf("HMM: maximum running time = none\n");
279 296
 	} else {
280 297
 		//FILE_LOG(logINFO) << "maximum running time = " << *maxtime << " sec";
281
-		Rprintf("maximum running time = %d sec\n", *maxtime);
298
+		if (*verbosity>=1) Rprintf("HMM: maximum running time = %d sec\n", *maxtime);
282 299
 	}
283 300
 	//FILE_LOG(logINFO) << "epsilon = " << *eps;
284
-	Rprintf("epsilon = %g\n", *eps);
301
+	if (*verbosity>=1) Rprintf("HMM: epsilon = %g\n", *eps);
285 302
 	//FILE_LOG(logINFO) << "number of modifications = " << *Nmod;
286
-	Rprintf("number of modifications = %d\n", *Nmod);
303
+	if (*verbosity>=1) Rprintf("HMM: number of modifications = %d\n", *Nmod);
287 304
 
288 305
 	// Flush Rprintf statements to console
289 306
 	R_FlushConsole();
... ...
@@ -303,7 +320,7 @@ void R_multivariate_hmm(int* O, int* T, int* N, int *Nmod, int* comb_states, dou
303 320
 
304 321
 	// Create the HMM
305 322
 	//FILE_LOG(logDEBUG1) << "Creating the multivariate HMM";
306
-	hmm = new ScaleHMM(*T, *N, *Nmod);
323
+	hmm = new ScaleHMM(*T, *N, *Nmod, *verbosity);
307 324
 	// Initialize the transition probabilities and proba
308 325
 	hmm->initialize_transition_probs(initial_A, *use_initial_params);
309 326
 	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
365 382
 	catch (std::exception& e)
366 383
 	{
367 384
 		//FILE_LOG(logERROR) << "Error in Baum-Welch: " << e.what();
368
-		Rprintf("Error in Baum-Welch: %s\n", e.what());
385
+		if (*verbosity>=1) Rprintf("HMM: Error in Baum-Welch: %s\n", e.what());
369 386
 		if (strcmp(e.what(),"nan detected")==0) { *error = 1; }
370 387
 		else { *error = 2; }
371 388
 	}
372 389
 	//FILE_LOG(logDEBUG1) << "Finished with Baum-Welch estimation";
373 390
 	
374
-	// Compute the posteriors and save results directly to the R pointer
391
+	// Get the posteriors and save results directly to the R pointer
375 392
 	if (*keep_posteriors == true)
376 393
 	{
377 394
 		//FILE_LOG(logDEBUG1) << "Recode posteriors into column representation";
378
-		Rprintf("Recoding posteriors ...\n");
395
+		if (*verbosity>=1) Rprintf("HMM: Recoding posteriors ...\n");
379 396
 		R_FlushConsole();
380 397
 		#pragma omp parallel for
381 398
 		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
387 404
 		}
388 405
 	}
389 406
 
407
+	// Get the densities and save results directly to the R pointer
408
+	if (*keep_densities == true)
409
+	{
410
+		//FILE_LOG(logDEBUG1) << "Recode posteriors into column representation";
411
+		if (*verbosity>=1) Rprintf("HMM: Recoding densities ...\n");
412
+		R_FlushConsole();
413
+		#pragma omp parallel for
414
+		for (int iN=0; iN<*N; iN++)
415
+		{
416
+			for (int t=0; t<*T; t++)
417
+			{
418
+				densities[t + iN * (*T)] = hmm->get_density(iN, t);
419
+			}
420
+		}
421
+	}
422
+
390 423
 	// Compute the states from posteriors
391 424
 	//FILE_LOG(logDEBUG1) << "Computing states from posteriors";
392 425
 // 	if (*fdr == -1)
... ...
@@ -401,7 +434,7 @@ void R_multivariate_hmm(int* O, int* T, int* N, int *Nmod, int* comb_states, dou
401 434
 			}
402 435
 			ind_max = std::distance(posterior_per_t.begin(), std::max_element(posterior_per_t.begin(), posterior_per_t.end()));
403 436
 			states[t] = comb_states[ind_max];
404
-// 			Rprintf("comb_states[t=%d, ind_max=%d] = %d\n", t, ind_max, states[t]);
437
+// 			if (*verbosity>=1) Rprintf("HMM: comb_states[t=%d, ind_max=%d] = %d\n", t, ind_max, states[t]);
405 438
 		}
406 439
 // 	}
407 440
 // 	else
Browse code

corrected state calling in multivariate

chakalakka authored on 09/02/2015 11:38:37
Showing1 changed files
... ...
@@ -391,6 +391,7 @@ void R_multivariate_hmm(int* O, int* T, int* N, int *Nmod, int* comb_states, dou
391 391
 	//FILE_LOG(logDEBUG1) << "Computing states from posteriors";
392 392
 // 	if (*fdr == -1)
393 393
 // 	{
394
+		int ind_max;
394 395
 		std::vector<double> posterior_per_t(*N);
395 396
 		for (int t=0; t<*T; t++)
396 397
 		{
... ...
@@ -398,7 +399,9 @@ void R_multivariate_hmm(int* O, int* T, int* N, int *Nmod, int* comb_states, dou
398 399
 			{
399 400
 				posterior_per_t[iN] = hmm->get_posterior(iN, t);
400 401
 			}
401
-			states[t] = comb_states[(int)*std::max_element(posterior_per_t.begin(), posterior_per_t.end())];
402
+			ind_max = std::distance(posterior_per_t.begin(), std::max_element(posterior_per_t.begin(), posterior_per_t.end()));
403
+			states[t] = comb_states[ind_max];
404
+// 			Rprintf("comb_states[t=%d, ind_max=%d] = %d\n", t, ind_max, states[t]);
402 405
 		}
403 406
 // 	}
404 407
 // 	else
Browse code

removed logging capability to pass R CMD check

chakalakka authored on 28/01/2015 14:42:46
Showing1 changed files
... ...
@@ -15,49 +15,49 @@ void R_univariate_hmm(int* O, int* T, int* N, double* size, double* prob, int* m
15 15
 {
16 16
 
17 17
 	// Define logging level
18
-// 	FILE* pFile = fopen("chromStar.log", "w");
19
-// 	Output2FILE::Stream() = pFile;
20
- 	FILELog::ReportingLevel() = FILELog::FromString("ERROR");
21
-//  	FILELog::ReportingLevel() = FILELog::FromString("DEBUG2");
18
+// 	//FILE* pFile = fopen("chromStar.log", "w");
19
+// 	Output2//FILE::Stream() = pFile;
20
+ 	//FILELog::ReportingLevel() = //FILELog::FromString("ERROR");
21
+//  	//FILELog::ReportingLevel() = //FILELog::FromString("DEBUG2");
22 22
 
23
-	FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;
23
+	//FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;
24 24
 	// Parallelization settings
25 25
 	omp_set_num_threads(*num_threads);
26 26
 
27 27
 	// Print some information
28
-	FILE_LOG(logINFO) << "number of states = " << *N;
28
+	//FILE_LOG(logINFO) << "number of states = " << *N;
29 29
 	Rprintf("number of states = %d\n", *N);
30
-	FILE_LOG(logINFO) << "number of bins = " << *T;
30
+	//FILE_LOG(logINFO) << "number of bins = " << *T;
31 31
 	Rprintf("number of bins = %d\n", *T);
32 32
 	if (*maxiter < 0)
33 33
 	{
34
-		FILE_LOG(logINFO) << "maximum number of iterations = none";
34
+		//FILE_LOG(logINFO) << "maximum number of iterations = none";
35 35
 		Rprintf("maximum number of iterations = none\n");
36 36
 	} else {
37
-		FILE_LOG(logINFO) << "maximum number of iterations = " << *maxiter;
37
+		//FILE_LOG(logINFO) << "maximum number of iterations = " << *maxiter;
38 38
 		Rprintf("maximum number of iterations = %d\n", *maxiter);
39 39
 	}
40 40
 	if (*maxtime < 0)
41 41
 	{
42
-		FILE_LOG(logINFO) << "maximum running time = none";
42
+		//FILE_LOG(logINFO) << "maximum running time = none";
43 43
 		Rprintf("maximum running time = none\n");
44 44
 	} else {
45
-		FILE_LOG(logINFO) << "maximum running time = " << *maxtime << " sec";
45
+		//FILE_LOG(logINFO) << "maximum running time = " << *maxtime << " sec";
46 46
 		Rprintf("maximum running time = %d sec\n", *maxtime);
47 47
 	}
48
-	FILE_LOG(logINFO) << "epsilon = " << *eps;
48
+	//FILE_LOG(logINFO) << "epsilon = " << *eps;
49 49
 	Rprintf("epsilon = %g\n", *eps);
50 50
 
51
-	FILE_LOG(logDEBUG3) << "observation vector";
51
+	//FILE_LOG(logDEBUG3) << "observation vector";
52 52
 	for (int t=0; t<50; t++) {
53
-		FILE_LOG(logDEBUG3) << "O["<<t<<"] = " << O[t];
53
+		//FILE_LOG(logDEBUG3) << "O["<<t<<"] = " << O[t];
54 54
 	}
55 55
 
56 56
 	// Flush Rprintf statements to console
57 57
 	R_FlushConsole();
58 58
 
59 59
 	// Create the HMM
60
-	FILE_LOG(logDEBUG1) << "Creating a univariate HMM";
60
+	//FILE_LOG(logDEBUG1) << "Creating a univariate HMM";
61 61
 	hmm = new ScaleHMM(*T, *N);
62 62
 	hmm->set_cutoff(*read_cutoff);
63 63
 	// Initialize the transition probabilities and proba
... ...
@@ -76,7 +76,7 @@ void R_univariate_hmm(int* O, int* T, int* N, double* size, double* prob, int* m
76 76
 		variance+= pow(O[t] - mean, 2);
77 77
 	}
78 78
 	variance = variance / *T;
79
-	FILE_LOG(logINFO) << "data mean = " << mean << ", data variance = " << variance;		
79
+	//FILE_LOG(logINFO) << "data mean = " << mean << ", data variance = " << variance;		
80 80
 	Rprintf("data mean = %g, data variance = %g\n", mean, variance);		
81 81
 	
82 82
 	// Go through all states of the hmm and assign the density functions
... ...
@@ -84,12 +84,12 @@ void R_univariate_hmm(int* O, int* T, int* N, double* size, double* prob, int* m
84 84
 	for (int i_state=0; i_state<*N; i_state++)
85 85
 	{
86 86
 		if (*use_initial_params) {
87
-			FILE_LOG(logINFO) << "Using given parameters for size and prob";
87
+			//FILE_LOG(logINFO) << "Using given parameters for size and prob";
88 88
 			Rprintf("Using given parameters for size and prob\n");
89 89
 			imean = (1-initial_prob[i_state])*initial_size[i_state] / initial_prob[i_state];
90 90
 			ivariance = imean / initial_prob[i_state];
91
-			FILE_LOG(logDEBUG2) << "imean = " << imean;
92
-			FILE_LOG(logDEBUG2) << "ivariance = " << ivariance;
91
+			//FILE_LOG(logDEBUG2) << "imean = " << imean;
92
+			//FILE_LOG(logDEBUG2) << "ivariance = " << ivariance;
93 93
 		} else {
94 94
 
95 95
 			if (*iniproc == 1)
... ...
@@ -97,13 +97,13 @@ void R_univariate_hmm(int* O, int* T, int* N, double* size, double* prob, int* m
97 97
 				// Simple initialization, seems to give the fastest convergence
98 98
 				if (i_state == 1)
99 99
 				{
100
-					FILE_LOG(logDEBUG) << "Initializing size and prob for state 1";
100
+					//FILE_LOG(logDEBUG) << "Initializing size and prob for state 1";
101 101
 					imean = mean;
102 102
 					ivariance = variance;
103 103
 				}
104 104
 				else if (i_state == 2)
105 105
 				{
106
-					FILE_LOG(logDEBUG) << "Initializing size and prob for state 2";
106
+					//FILE_LOG(logDEBUG) << "Initializing size and prob for state 2";
107 107
 					imean = mean+1;
108 108
 					ivariance = variance;
109 109
 				} 
... ...
@@ -116,7 +116,7 @@ void R_univariate_hmm(int* O, int* T, int* N, double* size, double* prob, int* m
116 116
 			else if (*iniproc == 2)
117 117
 			{
118 118
 				// Disturb mean and variance for use as randomized initial parameters
119
-				FILE_LOG(logINFO) << "Using random initialization for size and prob";
119
+				//FILE_LOG(logINFO) << "Using random initialization for size and prob";
120 120
 				Rprintf("Using random initialization for size and prob\n");
121 121
 				srand (clock());
122 122
 				int rand1, rand2;
... ...
@@ -124,25 +124,25 @@ void R_univariate_hmm(int* O, int* T, int* N, double* size, double* prob, int* m
124 124
 				rand2 = rand();
125 125
 				imean = (double)rand1/(double)RAND_MAX * 10*mean;
126 126
 				ivariance = imean + (double)rand2/(double)RAND_MAX * 20*imean; // variance has to be greater than mean, otherwise r will be negative
127
-				FILE_LOG(logDEBUG2) << "RAND_MAX = " << RAND_MAX;
128
-				FILE_LOG(logDEBUG2) << "rand1 = " << rand1;
129
-				FILE_LOG(logDEBUG2) << "rand2 = " << rand2;
130
-				FILE_LOG(logDEBUG2) << "imean = " << imean;
131
-				FILE_LOG(logDEBUG2) << "ivariance = " << ivariance;
127
+				//FILE_LOG(logDEBUG2) << "RAND_MAX = " << RAND_MAX;
128
+				//FILE_LOG(logDEBUG2) << "rand1 = " << rand1;
129
+				//FILE_LOG(logDEBUG2) << "rand2 = " << rand2;
130
+				//FILE_LOG(logDEBUG2) << "imean = " << imean;
131
+				//FILE_LOG(logDEBUG2) << "ivariance = " << ivariance;
132 132
 			}
133 133
 			else if (*iniproc == 3)
134 134
 			{
135 135
 				// Empirical initialization
136 136
 				if (i_state == 1)
137 137
 				{
138
-					FILE_LOG(logINFO) << "Initializing r and p empirically for state 1";
138
+					//FILE_LOG(logINFO) << "Initializing r and p empirically for state 1";
139 139
 					Rprintf("Initializing r and p empirically for state 1\n");
140 140
 					imean = mean/2;
141 141
 					ivariance = imean*2;
142 142
 				}
143 143
 				else if (i_state == 2)
144 144
 				{
145
-					FILE_LOG(logINFO) << "Initializing r and p empirically for state 2";
145
+					//FILE_LOG(logINFO) << "Initializing r and p empirically for state 2";
146 146
 					Rprintf("Initializing r and p empirically for state 2\n");
147 147
 					imean = mean*2;
148 148
 					ivariance = imean*2;
... ...
@@ -157,19 +157,19 @@ void R_univariate_hmm(int* O, int* T, int* N, double* size, double* prob, int* m
157 157
 
158 158
 		if (i_state >= 1)
159 159
 		{
160
-			FILE_LOG(logDEBUG1) << "Using negative binomial for state " << i_state;
160
+			//FILE_LOG(logDEBUG1) << "Using negative binomial for state " << i_state;
161 161
 			NegativeBinomial *d = new NegativeBinomial(O, *T, initial_size[i_state], initial_prob[i_state]); // delete is done inside ~ScaleHMM()
162 162
 			hmm->densityFunctions.push_back(d);
163 163
 		}
164 164
 		else if (i_state == 0)
165 165
 		{
166
-			FILE_LOG(logDEBUG1) << "Using only zeros for state " << i_state;
166
+			//FILE_LOG(logDEBUG1) << "Using only zeros for state " << i_state;
167 167
 			ZeroInflation *d = new ZeroInflation(O, *T); // delete is done inside ~ScaleHMM()
168 168
 			hmm->densityFunctions.push_back(d);
169 169
 		}
170 170
 		else
171 171
 		{
172
-			FILE_LOG(logWARNING) << "Density not specified, using default negative binomial for state " << i_state;
172
+			//FILE_LOG(logWARNING) << "Density not specified, using default negative binomial for state " << i_state;
173 173
 			NegativeBinomial *d = new NegativeBinomial(O, *T, initial_size[i_state], initial_prob[i_state]);
174 174
 			hmm->densityFunctions.push_back(d);
175 175
 		}
... ...
@@ -179,22 +179,22 @@ void R_univariate_hmm(int* O, int* T, int* N, double* size, double* prob, int* m
179 179
 	R_FlushConsole();
180 180
 
181 181
 	// Do the Baum-Welch to estimate the parameters
182
-	FILE_LOG(logDEBUG1) << "Starting Baum-Welch estimation";
182
+	//FILE_LOG(logDEBUG1) << "Starting Baum-Welch estimation";
183 183
 	try
184 184
 	{
185 185
 		hmm->baumWelch(maxiter, maxtime, eps);
186 186
 	}
187 187
 	catch (std::exception& e)
188 188
 	{
189
-		FILE_LOG(logERROR) << "Error in Baum-Welch: " << e.what();
189
+		//FILE_LOG(logERROR) << "Error in Baum-Welch: " << e.what();
190 190
 		Rprintf("Error in Baum-Welch: %s\n", e.what());
191 191
 		if (strcmp(e.what(),"nan detected")==0) { *error = 1; }
192 192
 		else { *error = 2; }
193 193
 	}
194 194
 		
195
-	FILE_LOG(logDEBUG1) << "Finished with Baum-Welch estimation";
195
+	//FILE_LOG(logDEBUG1) << "Finished with Baum-Welch estimation";
196 196
 	// Compute the posteriors and save results directly to the R pointer
197
-	FILE_LOG(logDEBUG1) << "Recode posteriors into column representation";
197
+	//FILE_LOG(logDEBUG1) << "Recode posteriors into column representation";
198 198
 	Rprintf("Recoding posteriors ...\n");
199 199
 	R_FlushConsole();
200 200
 	#pragma omp parallel for
... ...
@@ -206,7 +206,7 @@ void R_univariate_hmm(int* O, int* T, int* N, double* size, double* prob, int* m
206 206
 		}
207 207
 	}
208 208
 
209
-	FILE_LOG(logDEBUG1) << "Return parameters";
209
+	//FILE_LOG(logDEBUG1) << "Return parameters";
210 210
 	// also return the estimated transition matrix and the initial probs
211 211
 	for (int i=0; i<*N; i++)
212 212
 	{
... ...
@@ -236,7 +236,7 @@ void R_univariate_hmm(int* O, int* T, int* N, double* size, double* prob, int* m
236 236
 	*loglik = hmm->get_logP();
237 237
 	hmm->calc_weights(weights);
238 238
 	
239
-	FILE_LOG(logDEBUG1) << "Deleting the hmm";
239
+	//FILE_LOG(logDEBUG1) << "Deleting the hmm";
240 240
 	delete hmm;
241 241
 	hmm = NULL; // assign NULL to defuse the additional delete in on.exit() call
242 242
 }
... ...
@@ -250,39 +250,39 @@ void R_multivariate_hmm(int* O, int* T, int* N, int *Nmod, int* comb_states, dou
250 250
 {
251 251
 
252 252
 	// Define logging level {"ERROR", "WARNING", "INFO", "ITERATION", "DEBUG", "DEBUG1", "DEBUG2", "DEBUG3", "DEBUG4"}
253
-// 	FILE* pFile = fopen("chromStar.log", "w");
254
-// 	Output2FILE::Stream() = pFile;
255
- 	FILELog::ReportingLevel() = FILELog::FromString("ERROR");
256
-//  	FILELog::ReportingLevel() = FILELog::FromString("DEBUG3");
253
+// 	//FILE* pFile = fopen("chromStar.log", "w");
254
+// 	Output2//FILE::Stream() = pFile;
255
+ 	//FILELog::ReportingLevel() = //FILELog::FromString("ERROR");
256
+//  	//FILELog::ReportingLevel() = //FILELog::FromString("DEBUG3");
257 257
 
258
-	FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;
258
+	//FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;
259 259
 	// Parallelization settings
260 260
 	omp_set_num_threads(*num_threads);
261 261
 
262 262
 	// Print some information
263
-	FILE_LOG(logINFO) << "number of states = " << *N;
263
+	//FILE_LOG(logINFO) << "number of states = " << *N;
264 264
 	Rprintf("number of states = %d\n", *N);
265
-	FILE_LOG(logINFO) << "number of bins = " << *T;
265
+	//FILE_LOG(logINFO) << "number of bins = " << *T;
266 266
 	Rprintf("number of bins = %d\n", *T);
267 267
 	if (*maxiter < 0)
268 268
 	{
269
-		FILE_LOG(logINFO) << "maximum number of iterations = none";
269
+		//FILE_LOG(logINFO) << "maximum number of iterations = none";
270 270
 		Rprintf("maximum number of iterations = none\n");
271 271
 	} else {
272
-		FILE_LOG(logINFO) << "maximum number of iterations = " << *maxiter;
272
+		//FILE_LOG(logINFO) << "maximum number of iterations = " << *maxiter;
273 273
 		Rprintf("maximum number of iterations = %d\n", *maxiter);
274 274
 	}
275 275
 	if (*maxtime < 0)
276 276
 	{
277
-		FILE_LOG(logINFO) << "maximum running time = none";
277
+		//FILE_LOG(logINFO) << "maximum running time = none";
278 278
 		Rprintf("maximum running time = none\n");
279 279
 	} else {
280
-		FILE_LOG(logINFO) << "maximum running time = " << *maxtime << " sec";
280
+		//FILE_LOG(logINFO) << "maximum running time = " << *maxtime << " sec";
281 281
 		Rprintf("maximum running time = %d sec\n", *maxtime);
282 282
 	}
283
-	FILE_LOG(logINFO) << "epsilon = " << *eps;
283
+	//FILE_LOG(logINFO) << "epsilon = " << *eps;
284 284
 	Rprintf("epsilon = %g\n", *eps);
285
-	FILE_LOG(logINFO) << "number of modifications = " << *Nmod;
285
+	//FILE_LOG(logINFO) << "number of modifications = " << *Nmod;
286 286
 	Rprintf("number of modifications = %d\n", *Nmod);
287 287
 
288 288
 	// Flush Rprintf statements to console
... ...
@@ -299,10 +299,10 @@ void R_multivariate_hmm(int* O, int* T, int* N, int *Nmod, int* comb_states, dou
299 299
 		}
300 300
 	}
301 301
 // 	dtime = clock() - clocktime;
302
-// 	FILE_LOG(logDEBUG1) << "recoding observation vector to matrix representation: " << dtime << " clicks";
302
+// 	//FILE_LOG(logDEBUG1) << "recoding observation vector to matrix representation: " << dtime << " clicks";
303 303
 
304 304
 	// Create the HMM
305
-	FILE_LOG(logDEBUG1) << "Creating the multivariate HMM";
305
+	//FILE_LOG(logDEBUG1) << "Creating the multivariate HMM";
306 306
 	hmm = new ScaleHMM(*T, *N, *Nmod);
307 307
 	// Initialize the transition probabilities and proba
308 308
 	hmm->initialize_transition_probs(initial_A, *use_initial_params);
... ...
@@ -311,15 +311,15 @@ void R_multivariate_hmm(int* O, int* T, int* N, int *Nmod, int* comb_states, dou
311 311
 	// Print logproba and A
312 312
 // 	for (int iN=0; iN<*N; iN++)
313 313
 // 	{
314
-// 		FILE_LOG(logDEBUG) << "proba["<<iN<<"] = " <<exp(hmm->logproba[iN]);
314
+// 		//FILE_LOG(logDEBUG) << "proba["<<iN<<"] = " <<exp(hmm->logproba[iN]);
315 315
 // 		for (int jN=0; jN<*N; jN++)
316 316
 // 		{
317
-// 			FILE_LOG(logDEBUG) << "A["<<iN<<"]["<<jN<<"] = " << hmm->A[iN][jN];
317
+// 			//FILE_LOG(logDEBUG) << "A["<<iN<<"]["<<jN<<"] = " << hmm->A[iN][jN];
318 318
 // 		}
319 319
 // 	}
320 320
 
321 321
 	// 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)
322
-	FILE_LOG(logDEBUG1) << "Preparing the binary_states vector";
322
+	//FILE_LOG(logDEBUG1) << "Preparing the binary_states vector";
323 323
 	bool **binary_states = CallocBoolMatrix(*N, *Nmod);
324 324
 	for(int iN=0; iN < *N; iN++) //for each comb state considered
325 325
 	{
... ...
@@ -332,7 +332,7 @@ void R_multivariate_hmm(int* O, int* T, int* N, int *Nmod, int* comb_states, dou
332 332
 	}
333 333
 
334 334
 	/* initialize the distributions */
335
-	FILE_LOG(logDEBUG1) << "Initializing the distributions";
335
+	//FILE_LOG(logDEBUG1) << "Initializing the distributions";
336 336
 	for (int iN=0; iN<*N; iN++) //for each combinatorial state
337 337
 	{
338 338
 		std::vector <Density*> tempMarginals;            
... ...
@@ -350,31 +350,31 @@ void R_multivariate_hmm(int* O, int* T, int* N, int *Nmod, int* comb_states, dou
350 350
 			tempMarginals.push_back(d);
351 351
 		}
352 352
 		//MVCopulaApproximation *tempMVdens = new MVCopulaApproximation(O, tempMarginals, &(cor_matrix_inv[iN*Nmod*Nmod]), det[iN]);
353
-		FILE_LOG(logDEBUG1) << "Calling MVCopulaApproximation for state " << iN;
353
+		//FILE_LOG(logDEBUG1) << "Calling MVCopulaApproximation for state " << iN;
354 354
 		MVCopulaApproximation *tempMVdens = new MVCopulaApproximation(multiO, *T, tempMarginals, &(cor_matrix_inv[iN**Nmod**Nmod]), det[iN]); // delete is done inside ~ScaleHMM()
355 355
 		hmm->densityFunctions.push_back(tempMVdens);
356 356
 	}
357 357
 	FreeBoolMatrix(binary_states, *N);
358 358
 	
359 359
 	// Estimate the parameters
360
-	FILE_LOG(logDEBUG1) << "Starting Baum-Welch estimation";
360
+	//FILE_LOG(logDEBUG1) << "Starting Baum-Welch estimation";
361 361
 	try
362 362
 	{
363 363
 		hmm->baumWelch(maxiter, maxtime, eps);
364 364
 	}
365 365
 	catch (std::exception& e)
366 366
 	{
367
-		FILE_LOG(logERROR) << "Error in Baum-Welch: " << e.what();
367
+		//FILE_LOG(logERROR) << "Error in Baum-Welch: " << e.what();
368 368
 		Rprintf("Error in Baum-Welch: %s\n", e.what());
369 369
 		if (strcmp(e.what(),"nan detected")==0) { *error = 1; }
370 370
 		else { *error = 2; }
371 371
 	}
372
-	FILE_LOG(logDEBUG1) << "Finished with Baum-Welch estimation";
372
+	//FILE_LOG(logDEBUG1) << "Finished with Baum-Welch estimation";
373 373
 	
374 374
 	// Compute the posteriors and save results directly to the R pointer
375 375
 	if (*keep_posteriors == true)
376 376
 	{
377
-		FILE_LOG(logDEBUG1) << "Recode posteriors into column representation";
377
+		//FILE_LOG(logDEBUG1) << "Recode posteriors into column representation";
378 378
 		Rprintf("Recoding posteriors ...\n");
379 379
 		R_FlushConsole();
380 380
 		#pragma omp parallel for
... ...
@@ -388,7 +388,7 @@ void R_multivariate_hmm(int* O, int* T, int* N, int *Nmod, int* comb_states, dou
388 388
 	}
389 389
 
390 390
 	// Compute the states from posteriors
391
-	FILE_LOG(logDEBUG1) << "Computing states from posteriors";
391
+	//FILE_LOG(logDEBUG1) << "Computing states from posteriors";
392 392
 // 	if (*fdr == -1)
393 393
 // 	{
394 394
 		std::vector<double> posterior_per_t(*N);
... ...
@@ -416,7 +416,7 @@ void R_multivariate_hmm(int* O, int* T, int* N, int *Nmod, int* comb_states, dou
416 416
 // 		}
417 417
 // 	}
418 418
 	
419
-	FILE_LOG(logDEBUG1) << "Return parameters";
419
+	//FILE_LOG(logDEBUG1) << "Return parameters";
420 420
 	// also return the estimated transition matrix and the initial probs
421 421
 	for (int i=0; i<*N; i++)
422 422
 	{
... ...
@@ -428,7 +428,7 @@ void R_multivariate_hmm(int* O, int* T, int* N, int *Nmod, int* comb_states, dou
428 428
 	}
429 429
 	*loglik = hmm->get_logP();
430 430
 
431
-	FILE_LOG(logDEBUG1) << "Deleting the hmm";
431
+	//FILE_LOG(logDEBUG1) << "Deleting the hmm";
432 432
 	delete hmm;
433 433
 	hmm = NULL; // assign NULL to defuse the additional delete in on.exit() call
434 434
 // 	FreeIntMatrix(multiO, *Nmod); // free on.exit() in R code
... ...
@@ -442,7 +442,7 @@ void R_multivariate_hmm(int* O, int* T, int* N, int *Nmod, int* comb_states, dou
442 442
 extern "C" {
443 443
 void R_univariate_cleanup()
444 444
 {
445
-// 	FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__; // This message will be shown if interrupt happens before start of C-code
445
+// 	//FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__; // This message will be shown if interrupt happens before start of C-code
446 446
 	delete hmm;
447 447
 }
448 448
 }
Browse code

cleaned compiler warnings

chakalakka authored on 26/01/2015 18:39:14
Showing1 changed files
... ...
@@ -2,6 +2,7 @@
2 2
 #include "scalehmm.h"
3 3
 #include <vector> // storing density functions in multivariate
4 4
 #include <omp.h> // parallelization options
5
+#include <string> // strcmp
5 6
 
6 7
 static ScaleHMM* hmm; // declare as static outside the function because we only need one and this enables memory-cleanup on R_CheckUserInterrupt()
7 8
 static int** multiO;
... ...
@@ -79,7 +80,7 @@ void R_univariate_hmm(int* O, int* T, int* N, double* size, double* prob, int* m
79 80
 	Rprintf("data mean = %g, data variance = %g\n", mean, variance);		
80 81
 	
81 82
 	// Go through all states of the hmm and assign the density functions
82
-	double imean, ivariance;
83
+	double imean=0, ivariance=0;
83 84
 	for (int i_state=0; i_state<*N; i_state++)
84 85
 	{
85 86
 		if (*use_initial_params) {
... ...
@@ -187,7 +188,7 @@ void R_univariate_hmm(int* O, int* T, int* N, double* size, double* prob, int* m
187 188
 	{
188 189
 		FILE_LOG(logERROR) << "Error in Baum-Welch: " << e.what();
189 190
 		Rprintf("Error in Baum-Welch: %s\n", e.what());
190
-		if (e.what()=="nan detected") { *error = 1; }
191
+		if (strcmp(e.what(),"nan detected")==0) { *error = 1; }
191 192
 		else { *error = 2; }
192 193
 	}
193 194
 		
... ...
@@ -365,7 +366,7 @@ void R_multivariate_hmm(int* O, int* T, int* N, int *Nmod, int* comb_states, dou
365 366
 	{
366 367
 		FILE_LOG(logERROR) << "Error in Baum-Welch: " << e.what();
367 368
 		Rprintf("Error in Baum-Welch: %s\n", e.what());
368
-		if (e.what()=="nan detected") { *error = 1; }
369
+		if (strcmp(e.what(),"nan detected")==0) { *error = 1; }
369 370
 		else { *error = 2; }
370 371
 	}
371 372
 	FILE_LOG(logDEBUG1) << "Finished with Baum-Welch estimation";
... ...
@@ -390,14 +391,14 @@ void R_multivariate_hmm(int* O, int* T, int* N, int *Nmod, int* comb_states, dou
390 391
 	FILE_LOG(logDEBUG1) << "Computing states from posteriors";
391 392
 // 	if (*fdr == -1)
392 393
 // 	{
393
-		double posterior_per_t [*N];
394
+		std::vector<double> posterior_per_t(*N);
394 395
 		for (int t=0; t<*T; t++)
395 396
 		{
396 397
 			for (int iN=0; iN<*N; iN++)
397 398
 			{
398 399
 				posterior_per_t[iN] = hmm->get_posterior(iN, t);
399 400
 			}
400
-			states[t] = comb_states[argMax(posterior_per_t, *N)];
401
+			states[t] = comb_states[(int)*std::max_element(posterior_per_t.begin(), posterior_per_t.end())];
401 402
 		}
402 403
 // 	}
403 404
 // 	else
Browse code

fixed FDR calculation

chakalakka authored on 26/01/2015 16:42:02
Showing1 changed files
... ...
@@ -388,15 +388,32 @@ void R_multivariate_hmm(int* O, int* T, int* N, int *Nmod, int* comb_states, dou
388 388
 
389 389
 	// Compute the states from posteriors
390 390
 	FILE_LOG(logDEBUG1) << "Computing states from posteriors";
391
-	double posterior_per_t [*N];
392
-	for (int t=0; t<*T; t++)
393
-	{
394
-		for (int iN=0; iN<*N; iN++)
391
+// 	if (*fdr == -1)
392
+// 	{
393
+		double posterior_per_t [*N];
394
+		for (int t=0; t<*T; t++)
395 395
 		{
396
-			posterior_per_t[iN] = hmm->get_posterior(iN, t);
396
+			for (int iN=0; iN<*N; iN++)
397
+			{
398
+				posterior_per_t[iN] = hmm->get_posterior(iN, t);
399
+			}
400
+			states[t] = comb_states[argMax(posterior_per_t, *N)];
397 401
 		}
398
-		states[t] = comb_states[argMax(posterior_per_t, *N)];
399
-	}
402
+// 	}
403
+// 	else
404
+// 	{
405
+// 		double** transformed_posteriors = CallocDoubleMatrix(*T, *Nmod);
406
+// 		for (int t=0; t<*T; t++)
407
+// 		{
408
+// 			for (int iN=0; iN<*N; iN++)
409
+// 			{
410
+// 				for (int iNmod=0; iNmod<*Nmod; iNmod++)
411
+// 				{
412
+// 					transformed_posteriors[t][iNmod] += (double)binary_states[iN][iNmod] * hmm->get_posterior(iN, t);
413
+// 				}
414
+// 			}
415
+// 		}
416
+// 	}
400 417
 	
401 418
 	FILE_LOG(logDEBUG1) << "Return parameters";
402 419
 	// also return the estimated transition matrix and the initial probs
Browse code

posteriors can be kept, FDR can be set

chakalakka authored on 16/01/2015 19:38:49
Showing1 changed files
... ...
@@ -245,7 +245,7 @@ void R_univariate_hmm(int* O, int* T, int* N, double* size, double* prob, int* m
245 245
 // This function takes parameters from R, creates a multivariate HMM object, creates the distributions, runs the Baum-Welch and returns the result to R.
246 246
 // =====================================================================================================================================================
247 247
 extern "C" {
248
-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, int* states, double* A, double* proba, double* loglik, double* initial_A, double* initial_proba, bool* use_initial_params, int* num_threads, int* error)
248
+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)
249 249
 {
250 250
 
251 251
 	// Define logging level {"ERROR", "WARNING", "INFO", "ITERATION", "DEBUG", "DEBUG1", "DEBUG2", "DEBUG3", "DEBUG4"}
... ...
@@ -370,15 +370,21 @@ void R_multivariate_hmm(int* O, int* T, int* N, int *Nmod, int* comb_states, dou
370 370
 	}
371 371
 	FILE_LOG(logDEBUG1) << "Finished with Baum-Welch estimation";
372 372
 	
373
-// 	// Compute the posteriors and save results directly to the R pointer
374
-// 	FILE_LOG(logDEBUG1) << "Recode posteriors into column representation";
375
-// 	for (int iN=0; iN<*N; iN++)
376
-// 	{
377
-// 		for (int t=0; t<*T; t++)
378
-// 		{
379
-// 			posteriors[t + iN * (*T)] = hmm->get_posterior(iN, t);
380
-// 		}
381
-// 	}
373
+	// Compute the posteriors and save results directly to the R pointer
374
+	if (*keep_posteriors == true)
375
+	{
376
+		FILE_LOG(logDEBUG1) << "Recode posteriors into column representation";
377
+		Rprintf("Recoding posteriors ...\n");
378
+		R_FlushConsole();
379
+		#pragma omp parallel for
380
+		for (int iN=0; iN<*N; iN++)
381
+		{
382
+			for (int t=0; t<*T; t++)
383
+			{
384
+				posteriors[t + iN * (*T)] = hmm->get_posterior(iN, t);
385
+			}
386
+		}
387
+	}
382 388
 
383 389
 	// Compute the states from posteriors
384 390
 	FILE_LOG(logDEBUG1) << "Computing states from posteriors";
... ...
@@ -418,7 +424,7 @@ void R_multivariate_hmm(int* O, int* T, int* N, int *Nmod, int* comb_states, dou
418 424
 extern "C" {
419 425
 void R_univariate_cleanup()
420 426
 {
421
-	FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;
427
+// 	FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__; // This message will be shown if interrupt happens before start of C-code
422 428
 	delete hmm;
423 429
 }
424 430
 }
... ...
@@ -426,7 +432,6 @@ void R_univariate_cleanup()
426 432
 extern "C" {
427 433
 void R_multivariate_cleanup(int* Nmod)
428 434
 {
429
-	FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;
430 435