Browse code

Aneufinder with stepsize

chakalakka authored on 17/05/2017 12:38:47
Showing 8 changed files

... ...
@@ -41,7 +41,7 @@
41 41
 #'## The following call produces plots and genome browser files for all BAM files in "my-data-folder"
42 42
 #'Aneufinder(inputfolder="my-data-folder", outputfolder="my-output-folder")}
43 43
 #'
44
-Aneufinder <- function(inputfolder, outputfolder, configfile=NULL, numCPU=1, reuse.existing.files=TRUE, binsizes=1e6, variable.width.reference=NULL, reads.per.bin=NULL, pairedEndReads=FALSE, assembly=NULL, chromosomes=NULL, remove.duplicate.reads=TRUE, min.mapq=10, blacklist=NULL, use.bamsignals=FALSE, reads.store=FALSE, correction.method=NULL, GC.BSgenome=NULL, method=c('dnacopy','HMM'), strandseq=FALSE, eps=0.1, max.time=60, max.iter=5000, num.trials=15, states=c('zero-inflation',paste0(0:10,'-somy')), most.frequent.state='2-somy', most.frequent.state.strandseq='1-somy', resolution=c(3,6), min.segwidth=2, bw=4*binsizes[1], pval=1e-8, cluster.plots=TRUE) {
44
+Aneufinder <- function(inputfolder, outputfolder, configfile=NULL, numCPU=1, reuse.existing.files=TRUE, binsizes=1e6, stepsizes=NULL, variable.width.reference=NULL, reads.per.bin=NULL, pairedEndReads=FALSE, assembly=NULL, chromosomes=NULL, remove.duplicate.reads=TRUE, min.mapq=10, blacklist=NULL, use.bamsignals=FALSE, reads.store=FALSE, correction.method=NULL, GC.BSgenome=NULL, method=c('dnacopy','HMM'), strandseq=FALSE, eps=0.1, max.time=60, max.iter=5000, num.trials=15, states=c('zero-inflation',paste0(0:10,'-somy')), most.frequent.state='2-somy', most.frequent.state.strandseq='1-somy', resolution=c(3,6), min.segwidth=2, bw=4*binsizes[1], pval=1e-8, cluster.plots=TRUE) {
45 45
 
46 46
 #=======================
47 47
 ### Helper functions ###
... ...
@@ -87,7 +87,7 @@ if (class(GC.BSgenome)=='BSgenome') {
87 87
 numCPU <- as.numeric(numCPU)
88 88
 
89 89
 ## Put options into list and merge with conf
90
-params <- list(numCPU=numCPU, reuse.existing.files=reuse.existing.files, binsizes=binsizes, variable.width.reference=variable.width.reference, reads.per.bin=reads.per.bin, pairedEndReads=pairedEndReads, assembly=assembly, chromosomes=chromosomes, remove.duplicate.reads=remove.duplicate.reads, min.mapq=min.mapq, blacklist=blacklist, reads.store=reads.store, use.bamsignals=use.bamsignals, correction.method=correction.method, GC.BSgenome=GC.BSgenome, method=method, strandseq=strandseq, eps=eps, max.time=max.time, max.iter=max.iter, num.trials=num.trials, states=states, most.frequent.state=most.frequent.state, most.frequent.state.strandseq=most.frequent.state.strandseq, min.segwidth=min.segwidth, resolution=resolution, bw=bw, pval=pval, cluster.plots=cluster.plots)
90
+params <- list(numCPU=numCPU, reuse.existing.files=reuse.existing.files, binsizes=binsizes, stepsizes=stepsizes, variable.width.reference=variable.width.reference, reads.per.bin=reads.per.bin, pairedEndReads=pairedEndReads, assembly=assembly, chromosomes=chromosomes, remove.duplicate.reads=remove.duplicate.reads, min.mapq=min.mapq, blacklist=blacklist, reads.store=reads.store, use.bamsignals=use.bamsignals, correction.method=correction.method, GC.BSgenome=GC.BSgenome, method=method, strandseq=strandseq, eps=eps, max.time=max.time, max.iter=max.iter, num.trials=num.trials, states=states, most.frequent.state=most.frequent.state, most.frequent.state.strandseq=most.frequent.state.strandseq, min.segwidth=min.segwidth, resolution=resolution, bw=bw, pval=pval, cluster.plots=cluster.plots)
91 91
 conf <- c(conf, params[setdiff(names(params),names(conf))])
92 92
 
93 93
 ## Check user input
... ...
@@ -108,6 +108,7 @@ if (any(formats == 'bed') & is.null(conf[['assembly']])) {
108 108
 
109 109
 ## Helpers
110 110
 binsizes <- conf[['binsizes']]
111
+stepsizes <- conf[['stepsizes']]
111 112
 reads.per.bins <- conf[['reads.per.bin']]
112 113
 patterns <- c(paste0('reads.per.bin_',reads.per.bins,'_'), paste0('binsize_',format(binsizes, scientific=TRUE, trim=TRUE),'_'))
113 114
 patterns <- setdiff(patterns, c('reads.per.bin__','binsize__'))
... ...
@@ -206,9 +207,9 @@ if (!is.null(conf[['variable.width.reference']])) {
206 207
 	} else if (format == 'bed') {
207 208
 		reads <- bed2GRanges(conf[['variable.width.reference']], assembly=chrom.lengths.df, chromosomes=conf[['chromosomes']], remove.duplicate.reads=conf[['remove.duplicate.reads']], min.mapq=conf[['min.mapq']], blacklist=conf[['blacklist']])
208 209
 	}
209
-	bins <- variableWidthBins(reads, binsizes=conf[['binsizes']], chromosomes=conf[['chromosomes']])
210
+	bins <- variableWidthBins(reads, binsizes=conf[['binsizes']], stepsizes=conf[['stepsizes']], chromosomes=conf[['chromosomes']])
210 211
 } else {
211
-  bins <- fixedWidthBins(chrom.lengths=chrom.lengths, chromosomes=conf[['chromosomes']], binsizes=conf[['binsizes']])
212
+  bins <- fixedWidthBins(chrom.lengths=chrom.lengths, chromosomes=conf[['chromosomes']], binsizes=conf[['binsizes']], stepsizes=conf[['stepsizes']])
212 213
 }
213 214
 message("==| Finished making bins.")
214 215
 
... ...
@@ -315,6 +315,7 @@ binReads <- function(file, assembly, ID=basename(file), bamindex=file, chromosom
315 315
 
316 316
 			### ID ###
317 317
 			attr(bins, 'ID') <- ID
318
+			stopTimedMessage(ptm)
318 319
 
319 320
 			### Save or return the bins ###
320 321
 			if (save.as.RData==TRUE) {
... ...
@@ -326,7 +327,6 @@ binReads <- function(file, assembly, ID=basename(file), bamindex=file, chromosom
326 327
 			} else {
327 328
 				bins.list[[ibinsize]] <- bins
328 329
 			}
329
-			stopTimedMessage(ptm)
330 330
 
331 331
 	} ### end loop binsizes ###
332 332
 
... ...
@@ -23,7 +23,7 @@
23 23
 #'## Check the fit
24 24
 #'plot(model, type='histogram')
25 25
 #'
26
-findCNVs <- function(binned.data, ID=NULL, eps=0.1, init="standard", max.time=-1, max.iter=1000, num.trials=15, eps.try=10*eps, num.threads=1, count.cutoff.quantile=0.999, strand='*', states=c("zero-inflation",paste0(0:10,"-somy")), most.frequent.state="2-somy", method="HMM", algorithm="EM", initial.params=NULL) {
26
+findCNVs <- function(binned.data, ID=NULL, eps=0.1, init="standard", max.time=-1, max.iter=1000, num.trials=15, eps.try=10*eps, num.threads=1, count.cutoff.quantile=0.999, strand='*', states=c("zero-inflation",paste0(0:10,"-somy")), most.frequent.state="2-somy", method="HMM", algorithm="EM", initial.params=NULL, verbosity=1) {
27 27
 
28 28
 	## Intercept user input
29 29
   binned.data <- loadFromFiles(binned.data, check.class=c('GRanges', 'GRangesList'))[[1]]
... ...
@@ -46,7 +46,7 @@ findCNVs <- function(binned.data, ID=NULL, eps=0.1, init="standard", max.time=-1
46 46
 	message("Method = ", method)
47 47
 
48 48
 	if (method == 'HMM') {
49
-		model <- univariate.findCNVs(binned.data, ID, eps=eps, init=init, max.time=max.time, max.iter=max.iter, num.trials=num.trials, eps.try=eps.try, num.threads=num.threads, count.cutoff.quantile=count.cutoff.quantile, strand=strand, states=states, most.frequent.state=most.frequent.state, algorithm=algorithm, initial.params=initial.params)
49
+		model <- univariate.findCNVs(binned.data, ID, eps=eps, init=init, max.time=max.time, max.iter=max.iter, num.trials=num.trials, eps.try=eps.try, num.threads=num.threads, count.cutoff.quantile=count.cutoff.quantile, strand=strand, states=states, most.frequent.state=most.frequent.state, algorithm=algorithm, initial.params=initial.params, verbosity=verbosity)
50 50
 	} else if (method == 'dnacopy') {
51 51
 	  model <- DNAcopy.findCNVs(binned.data, ID, CNgrid.start=1.5, count.cutoff.quantile=count.cutoff.quantile, strand=strand)
52 52
 	}
... ...
@@ -82,9 +82,10 @@ findCNVs <- function(binned.data, ID=NULL, eps=0.1, init="standard", max.time=-1
82 82
 #' @param most.frequent.state One of the states that were given in \code{states}. The specified state is assumed to be the most frequent one. This can help the fitting procedure to converge into the correct fit.
83 83
 #' @param algorithm One of \code{c('baumWelch','EM')}. The expectation maximization (\code{'EM'}) will find the most likely states and fit the best parameters to the data, the \code{'baumWelch'} will find the most likely states using the initial parameters.
84 84
 #' @param initial.params A \code{\link{aneuHMM}} object or file containing such an object from which initial starting parameters will be extracted.
85
+#' @param verbosity Integer specifying the verbosity of printed messages.
85 86
 #' @return An \code{\link{aneuHMM}} object.
86 87
 #' @importFrom stats runif
87
-univariate.findCNVs <- function(binned.data, ID=NULL, eps=0.1, init="standard", max.time=-1, max.iter=-1, num.trials=1, eps.try=NULL, num.threads=1, count.cutoff.quantile=0.999, strand='*', states=c("zero-inflation",paste0(0:10,"-somy")), most.frequent.state="2-somy", algorithm="EM", initial.params=NULL) {
88
+univariate.findCNVs <- function(binned.data, ID=NULL, eps=0.1, init="standard", max.time=-1, max.iter=-1, num.trials=1, eps.try=NULL, num.threads=1, count.cutoff.quantile=0.999, strand='*', states=c("zero-inflation",paste0(0:10,"-somy")), most.frequent.state="2-somy", algorithm="EM", initial.params=NULL, verbosity=1) {
88 89
 
89 90
 	### Define cleanup behaviour ###
90 91
 	on.exit(.C("C_univariate_cleanup", PACKAGE = 'AneuFinder'))
... ...
@@ -183,6 +184,7 @@ univariate.findCNVs <- function(binned.data, ID=NULL, eps=0.1, init="standard",
183 184
       init <- 'initial.params'
184 185
     	algorithm <- factor('baumWelch', levels=c('baumWelch','viterbi','EM'))
185 186
       num.trials <- 1
187
+      verbosity <- 0
186 188
     }
187 189
 
188 190
   	# Check if there are counts in the data, otherwise HMM will blow up
... ...
@@ -213,7 +215,7 @@ univariate.findCNVs <- function(binned.data, ID=NULL, eps=0.1, init="standard",
213 215
   	## Call univariate in a for loop to enable multiple trials
214 216
   	modellist <- list()
215 217
   	for (i_try in 1:num.trials) {
216
-  		message(paste0("Trial ",i_try," / ",num.trials))
218
+  		if (verbosity >= 1) { message(paste0("Trial ",i_try," / ",num.trials)) }
217 219
   
218 220
   		## Initial parameters
219 221
   		if (init == 'initial.params') {
... ...
@@ -305,6 +307,7 @@ univariate.findCNVs <- function(binned.data, ID=NULL, eps=0.1, init="standard",
305 307
   			error = as.integer(0), # int* error (error handling)
306 308
   			count.cutoff = as.integer(count.cutoff), # int* count.cutoff
307 309
   			algorithm = as.integer(algorithm), # int* algorithm
310
+  			verbosity = as.integer(verbosity), # int* verbosity
308 311
   			PACKAGE = 'AneuFinder'
309 312
   		)
310 313
   
... ...
@@ -412,6 +415,7 @@ univariate.findCNVs <- function(binned.data, ID=NULL, eps=0.1, init="standard",
412 415
   				error = as.integer(0), # int* error (error handling)
413 416
   				count.cutoff = as.integer(count.cutoff), # int* count.cutoff
414 417
   				algorithm = as.integer(algorithm), # int* algorithm
418
+    			verbosity = as.integer(verbosity), # int* verbosity
415 419
   				PACKAGE = 'AneuFinder'
416 420
   			)
417 421
   		}
... ...
@@ -558,7 +562,7 @@ univariate.findCNVs <- function(binned.data, ID=NULL, eps=0.1, init="standard",
558 562
 #' @inheritParams findCNVs
559 563
 #' @return An \code{\link{aneuBiHMM}} object.
560 564
 #' @importFrom stats pgeom pnbinom qnorm
561
-bivariate.findCNVs <- function(binned.data, ID=NULL, eps=0.1, init="standard", max.time=-1, max.iter=-1, num.trials=1, eps.try=NULL, num.threads=1, count.cutoff.quantile=0.999, states=c("zero-inflation",paste0(0:10,"-somy")), most.frequent.state="1-somy", algorithm='EM', initial.params=NULL) {
565
+bivariate.findCNVs <- function(binned.data, ID=NULL, eps=0.1, init="standard", max.time=-1, max.iter=-1, num.trials=1, eps.try=NULL, num.threads=1, count.cutoff.quantile=0.999, states=c("zero-inflation",paste0(0:10,"-somy")), most.frequent.state="1-somy", algorithm='EM', initial.params=NULL, verbosity=1) {
562 566
 
563 567
 	## Intercept user input
564 568
   binned.data <- loadFromFiles(binned.data, check.class='GRanges')[[1]]
... ...
@@ -864,6 +868,7 @@ bivariate.findCNVs <- function(binned.data, ID=NULL, eps=0.1, init="standard", m
864 868
 		num.threads = as.integer(num.threads), # int* num_threads
865 869
 		error = as.integer(0), # error handling
866 870
 		algorithm = as.integer(algorithm), # int* algorithm
871
+		verbosity = as.integer(verbosity), # int* verbosity
867 872
 		PACKAGE = 'AneuFinder'
868 873
 		)
869 874
 			
... ...
@@ -72,7 +72,7 @@ writeConfig <- function(conf, configfile) {
72 72
 		cat(i1," = ",formatstring(conf[[i1]]),"\n", file=f)
73 73
 	}
74 74
 	cat("\n[Binning]\n", file=f)
75
-	for (i1 in c('binsizes', 'variable.width.reference', 'reads.per.bin', 'pairedEndReads', 'assembly', 'chromosomes', 'remove.duplicate.reads', 'min.mapq', 'blacklist', 'reads.store', 'use.bamsignals')) {
75
+	for (i1 in c('binsizes', 'stepsizes', 'variable.width.reference', 'reads.per.bin', 'pairedEndReads', 'assembly', 'chromosomes', 'remove.duplicate.reads', 'min.mapq', 'blacklist', 'reads.store', 'use.bamsignals')) {
76 76
 		cat(i1," = ",formatstring(conf[[i1]]),"\n", file=f)
77 77
 	}
78 78
 	cat("\n[Correction]\n", file=f)
... ...
@@ -5,7 +5,7 @@
5 5
 \title{Wrapper function for the \code{\link{AneuFinder}} package}
6 6
 \usage{
7 7
 Aneufinder(inputfolder, outputfolder, configfile = NULL, numCPU = 1,
8
-  reuse.existing.files = TRUE, binsizes = 1e+06,
8
+  reuse.existing.files = TRUE, binsizes = 1e+06, stepsizes = NULL,
9 9
   variable.width.reference = NULL, reads.per.bin = NULL,
10 10
   pairedEndReads = FALSE, assembly = NULL, chromosomes = NULL,
11 11
   remove.duplicate.reads = TRUE, min.mapq = 10, blacklist = NULL,
... ...
@@ -30,6 +30,8 @@ Aneufinder(inputfolder, outputfolder, configfile = NULL, numCPU = 1,
30 30
 
31 31
 \item{binsizes}{An integer vector with bin sizes. If more than one value is given, output files will be produced for each bin size.}
32 32
 
33
+\item{stepsizes}{A vector of step sizes the same length as \code{binsizes}.}
34
+
33 35
 \item{variable.width.reference}{A BAM file that is used as reference to produce variable width bins. See \code{\link{variableWidthBins}} for details.}
34 36
 
35 37
 \item{reads.per.bin}{Approximate number of desired reads per bin. The bin size will be selected accordingly. Output files are produced for each value.}
... ...
@@ -9,7 +9,7 @@ static double** multiD;
9 9
 // ===================================================================================================================================================
10 10
 // This function takes parameters from R, creates a univariate HMM object, creates the distributions, runs the EM and returns the result to R.
11 11
 // ===================================================================================================================================================
12
-void univariate_hmm(int* O, int* T, int* N, int* state_labels, double* size, double* prob, int* maxiter, int* maxtime, double* eps, double* maxPosterior, int* states, double* A, double* proba, double* loglik, double* weights, int* distr_type, double* initial_size, double* initial_prob, double* initial_A, double* initial_proba, bool* use_initial_params, int* num_threads, int* error, int* read_cutoff, int* algorithm)
12
+void univariate_hmm(int* O, int* T, int* N, int* state_labels, double* size, double* prob, int* maxiter, int* maxtime, double* eps, double* maxPosterior, int* states, double* A, double* proba, double* loglik, double* weights, int* distr_type, double* initial_size, double* initial_prob, double* initial_A, double* initial_proba, bool* use_initial_params, int* num_threads, int* error, int* read_cutoff, int* algorithm, int* verbosity)
13 13
 {
14 14
 
15 15
 	// Define logging level
... ...
@@ -23,34 +23,34 @@ void univariate_hmm(int* O, int* T, int* N, int* state_labels, double* size, dou
23 23
 
24 24
 	// Print some information
25 25
 	//FILE_LOG(logINFO) << "number of states = " << *N;
26
-	Rprintf("number of states = %d\n", *N);
26
+	if (*verbosity>=1) Rprintf("number of states = %d\n", *N);
27 27
 	//FILE_LOG(logINFO) << "number of bins = " << *T;
28
-	Rprintf("number of bins = %d\n", *T);
28
+	if (*verbosity>=1) Rprintf("number of bins = %d\n", *T);
29 29
 	if (*maxiter < 0)
30 30
 	{
31 31
 		//FILE_LOG(logINFO) << "maximum number of iterations = none";
32
-		Rprintf("maximum number of iterations = none\n");
32
+		if (*verbosity>=1) Rprintf("maximum number of iterations = none\n");
33 33
 	} else {
34 34
 		//FILE_LOG(logINFO) << "maximum number of iterations = " << *maxiter;
35
-		Rprintf("maximum number of iterations = %d\n", *maxiter);
35
+		if (*verbosity>=1) Rprintf("maximum number of iterations = %d\n", *maxiter);
36 36
 	}
37 37
 	if (*maxtime < 0)
38 38
 	{
39 39
 		//FILE_LOG(logINFO) << "maximum running time = none";
40
-		Rprintf("maximum running time = none\n");
40
+		if (*verbosity>=1) Rprintf("maximum running time = none\n");
41 41
 	} else {
42 42
 		//FILE_LOG(logINFO) << "maximum running time = " << *maxtime << " sec";
43
-		Rprintf("maximum running time = %d sec\n", *maxtime);
43
+		if (*verbosity>=1) Rprintf("maximum running time = %d sec\n", *maxtime);
44 44
 	}
45 45
 	//FILE_LOG(logINFO) << "epsilon = " << *eps;
46
-	Rprintf("epsilon = %g\n", *eps);
46
+	if (*verbosity>=1) Rprintf("epsilon = %g\n", *eps);
47 47
 
48 48
 	//FILE_LOG(logDEBUG3) << "observation vector";
49 49
 	for (int t=0; t<50; t++) {
50 50
 		//FILE_LOG(logDEBUG3) << "O["<<t<<"] = " << O[t];
51 51
 	}
52 52
 
53
-	// Flush Rprintf statements to console
53
+	// Flush if (*verbosity>=1) Rprintf statements to console
54 54
 	R_FlushConsole();
55 55
 
56 56
 	// Create the HMM
... ...
@@ -75,7 +75,7 @@ void univariate_hmm(int* O, int* T, int* N, int* state_labels, double* size, dou
75 75
 	}
76 76
 	variance = variance / *T;
77 77
 	//FILE_LOG(logINFO) << "data mean = " << mean << ", data variance = " << variance;		
78
-	Rprintf("data mean = %g, data variance = %g\n", mean, variance);		
78
+	if (*verbosity>=1) Rprintf("data mean = %g, data variance = %g\n", mean, variance);		
79 79
 	
80 80
 	// Create the emission densities and initialize
81 81
 	for (int i_state=0; i_state<*N; i_state++)
... ...
@@ -112,7 +112,7 @@ void univariate_hmm(int* O, int* T, int* N, int* state_labels, double* size, dou
112 112
 		}
113 113
 	}
114 114
 
115
-	// Flush Rprintf statements to console
115
+	// Flush if (*verbosity>=1) Rprintf statements to console
116 116
 	R_FlushConsole();
117 117
 
118 118
 	// Do the EM to estimate the parameters
... ...
@@ -132,7 +132,7 @@ void univariate_hmm(int* O, int* T, int* N, int* state_labels, double* size, dou
132 132
 	catch (std::exception& e)
133 133
 	{
134 134
 		//FILE_LOG(logERROR) << "Error in EM/baumWelch: " << e.what();
135
-		Rprintf("Error in EM/baumWelch: %s\n", e.what());
135
+		if (*verbosity>=1) Rprintf("Error in EM/baumWelch: %s\n", e.what());
136 136
 		if (strcmp(e.what(),"nan detected")==0) { *error = 1; }
137 137
 		else { *error = 2; }
138 138
 	}
... ...
@@ -214,7 +214,7 @@ void univariate_hmm(int* O, int* T, int* N, int* state_labels, double* size, dou
214 214
 // =====================================================================================================================================================
215 215
 // This function takes parameters from R, creates a multivariate HMM object, runs the EM and returns the result to R.
216 216
 // =====================================================================================================================================================
217
-void multivariate_hmm(double* D, int* T, int* N, int *Nmod, int* comb_states, int* maxiter, int* maxtime, double* eps, int* states, double* A, double* proba, double* loglik, double* initial_A, double* initial_proba, bool* use_initial_params, int* num_threads, int* error, int* algorithm)
217
+void multivariate_hmm(double* D, int* T, int* N, int *Nmod, int* comb_states, int* maxiter, int* maxtime, double* eps, int* states, double* A, double* proba, double* loglik, double* initial_A, double* initial_proba, bool* use_initial_params, int* num_threads, int* error, int* algorithm, int* verbosity)
218 218
 {
219 219
 
220 220
 	// Define logging level {"ERROR", "WARNING", "INFO", "ITERATION", "DEBUG", "DEBUG1", "DEBUG2", "DEBUG3", "DEBUG4"}
... ...
@@ -229,31 +229,31 @@ void multivariate_hmm(double* D, int* T, int* N, int *Nmod, int* comb_states, in
229 229
 
230 230
 	// Print some information
231 231
 	//FILE_LOG(logINFO) << "number of states = " << *N;
232
-	Rprintf("number of states = %d\n", *N);
232
+	if (*verbosity>=1) Rprintf("number of states = %d\n", *N);
233 233
 	//FILE_LOG(logINFO) << "number of bins = " << *T;
234
-	Rprintf("number of bins = %d\n", *T);
234
+	if (*verbosity>=1) Rprintf("number of bins = %d\n", *T);
235 235
 	if (*maxiter < 0)
236 236
 	{
237 237
 		//FILE_LOG(logINFO) << "maximum number of iterations = none";
238
-		Rprintf("maximum number of iterations = none\n");
238
+		if (*verbosity>=1) Rprintf("maximum number of iterations = none\n");
239 239
 	} else {
240 240
 		//FILE_LOG(logINFO) << "maximum number of iterations = " << *maxiter;
241
-		Rprintf("maximum number of iterations = %d\n", *maxiter);
241
+		if (*verbosity>=1) Rprintf("maximum number of iterations = %d\n", *maxiter);
242 242
 	}
243 243
 	if (*maxtime < 0)
244 244
 	{
245 245
 		//FILE_LOG(logINFO) << "maximum running time = none";
246
-		Rprintf("maximum running time = none\n");
246
+		if (*verbosity>=1) Rprintf("maximum running time = none\n");
247 247
 	} else {
248 248
 		//FILE_LOG(logINFO) << "maximum running time = " << *maxtime << " sec";
249
-		Rprintf("maximum running time = %d sec\n", *maxtime);
249
+		if (*verbosity>=1) Rprintf("maximum running time = %d sec\n", *maxtime);
250 250
 	}
251 251
 	//FILE_LOG(logINFO) << "epsilon = " << *eps;
252
-	Rprintf("epsilon = %g\n", *eps);
252
+	if (*verbosity>=1) Rprintf("epsilon = %g\n", *eps);
253 253
 	//FILE_LOG(logINFO) << "number of modifications = " << *Nmod;
254
-	Rprintf("number of modifications = %d\n", *Nmod);
254
+	if (*verbosity>=1) Rprintf("number of modifications = %d\n", *Nmod);
255 255
 
256
-	// Flush Rprintf statements to console
256
+	// Flush if (*verbosity>=1) Rprintf statements to console
257 257
 	R_FlushConsole();
258 258
 
259 259
 	// Recode the densities vector to matrix representation
... ...
@@ -303,7 +303,7 @@ void multivariate_hmm(double* D, int* T, int* N, int *Nmod, int* comb_states, in
303 303
 	catch (std::exception& e)
304 304
 	{
305 305
 		//FILE_LOG(logERROR) << "Error in EM/baumWelch: " << e.what();
306
-		Rprintf("Error in EM/baumWelch: %s\n", e.what());
306
+		if (*verbosity>=1) Rprintf("Error in EM/baumWelch: %s\n", e.what());
307 307
 		if (strcmp(e.what(),"nan detected")==0) { *error = 1; }
308 308
 		else { *error = 2; }
309 309
 	}
... ...
@@ -379,10 +379,10 @@ void array2D_which_max(double* array2D, int* dim, int* ind_max, double* value_ma
379 379
     for (int i1=0; i1<dim[1]; i1++)
380 380
     {
381 381
 			value_per_i0[i1] = array2D[i1 * dim[0] + i0];
382
-      // Rprintf("i0=%d, i1=%d, value_per_i0[%d] = %g\n", i0, i1, i1, value_per_i0[i1]);
382
+      // if (*verbosity>=1) Rprintf("i0=%d, i1=%d, value_per_i0[%d] = %g\n", i0, i1, i1, value_per_i0[i1]);
383 383
     }
384 384
 		ind_max[i0] = 1 + std::distance(value_per_i0.begin(), std::max_element(value_per_i0.begin(), value_per_i0.end()));
385 385
     value_max[i0] = *std::max_element(value_per_i0.begin(), value_per_i0.end());
386 386
   }
387 387
 	
388
-}
389 388
\ No newline at end of file
389
+}
... ...
@@ -13,10 +13,10 @@
13 13
 // #endif
14 14
 
15 15
 extern "C"
16
-void univariate_hmm(int* O, int* T, int* N, int* state_labels, double* size, double* prob, int* maxiter, int* maxtime, double* eps, double* maxPosterior, int* states, double* A, double* proba, double* loglik, double* weights, int* distr_type, double* initial_size, double* initial_prob, double* initial_A, double* initial_proba, bool* use_initial_params, int* num_threads, int* error, int* read_cutoff, int* algorithm);
16
+void univariate_hmm(int* O, int* T, int* N, int* state_labels, double* size, double* prob, int* maxiter, int* maxtime, double* eps, double* maxPosterior, int* states, double* A, double* proba, double* loglik, double* weights, int* distr_type, double* initial_size, double* initial_prob, double* initial_A, double* initial_proba, bool* use_initial_params, int* num_threads, int* error, int* read_cutoff, int* algorithm, int* verbosity);
17 17
 
18 18
 extern "C"
19
-void multivariate_hmm(double* D, int* T, int* N, int *Nmod, int* comb_states, int* maxiter, int* maxtime, double* eps, int* states, double* A, double* proba, double* loglik, double* initial_A, double* initial_proba, bool* use_initial_params, int* num_threads, int* error, int* algorithm);
19
+void multivariate_hmm(double* D, int* T, int* N, int *Nmod, int* comb_states, int* maxiter, int* maxtime, double* eps, int* states, double* A, double* proba, double* loglik, double* initial_A, double* initial_proba, bool* use_initial_params, int* num_threads, int* error, int* algorithm, int* verbosity);
20 20
 
21 21
 extern "C"
22 22
 void univariate_cleanup();
... ...
@@ -3,14 +3,14 @@
3 3
 #include "R_interface.h"
4 4
 
5 5
 
6
-R_NativePrimitiveArgType arg1[] = {INTSXP, INTSXP, INTSXP, INTSXP, REALSXP, REALSXP, INTSXP, INTSXP, REALSXP, REALSXP, INTSXP, REALSXP, REALSXP, REALSXP, REALSXP, INTSXP, REALSXP, REALSXP, REALSXP, REALSXP, LGLSXP, INTSXP, INTSXP, INTSXP, INTSXP};
7
-R_NativePrimitiveArgType arg2[] = {REALSXP, INTSXP, INTSXP, INTSXP, INTSXP, INTSXP, INTSXP, REALSXP, INTSXP, REALSXP, REALSXP, REALSXP, REALSXP, REALSXP, LGLSXP, INTSXP, INTSXP, INTSXP};
6
+R_NativePrimitiveArgType arg1[] = {INTSXP, INTSXP, INTSXP, INTSXP, REALSXP, REALSXP, INTSXP, INTSXP, REALSXP, REALSXP, INTSXP, REALSXP, REALSXP, REALSXP, REALSXP, INTSXP, REALSXP, REALSXP, REALSXP, REALSXP, LGLSXP, INTSXP, INTSXP, INTSXP, INTSXP, INTSXP};
7
+R_NativePrimitiveArgType arg2[] = {REALSXP, INTSXP, INTSXP, INTSXP, INTSXP, INTSXP, INTSXP, REALSXP, INTSXP, REALSXP, REALSXP, REALSXP, REALSXP, REALSXP, LGLSXP, INTSXP, INTSXP, INTSXP, INTSXP};
8 8
 R_NativePrimitiveArgType arg4[] = {INTSXP};
9 9
 R_NativePrimitiveArgType arg5[] = {REALSXP, INTSXP, INTSXP, REALSXP};
10 10
 
11 11
 static const R_CMethodDef CEntries[]  = {
12
-    {"C_univariate_hmm", (DL_FUNC) &univariate_hmm, 25, arg1},
13
-    {"C_multivariate_hmm", (DL_FUNC) &multivariate_hmm, 18, arg2},
12
+    {"C_univariate_hmm", (DL_FUNC) &univariate_hmm, 26, arg1},
13
+    {"C_multivariate_hmm", (DL_FUNC) &multivariate_hmm, 19, arg2},
14 14
     {"C_univariate_cleanup", (DL_FUNC) &univariate_cleanup, 0, NULL},
15 15
     {"C_multivariate_cleanup", (DL_FUNC) &multivariate_cleanup, 1, arg4},
16 16
     {"C_array2D_which_max", (DL_FUNC) &array2D_which_max, 4, arg5},