@@ -6,8 +6,8 @@ Date: 2014-05-01

 Author: Aaron Taudt, David Porubsky

 Maintainer: Aaron Taudt <a.s.taudt@umcg.nl>, David Porubsky <d.porubsky@umcg.nl>

 Description: not yet

-Depends: R (>= 3.1.0), GenomicRanges, IRanges, ggplot2, reshape2

-Imports: Rsamtools, Biostrings, GenomicAlignments

+Depends: R (>= 3.1.0), GenomicRanges, IRanges, ggplot2, reshape2, biovizBase

+Imports: Rsamtools, Biostrings, GenomicAlignments, doParallel, foreach, grid

 Suggests:

 License: GPL

 LazyLoad: yes

@@ -36,5 +36,5 @@ mixedorder

 # Import all packages listed as Imports or Depends

 import(

-	GenomicRanges, IRanges, ggplot2, reshape2, Rsamtools, Biostrings, GenomicAlignments

+	GenomicRanges, IRanges, ggplot2, reshape2, Rsamtools, Biostrings, GenomicAlignments, doParallel, foreach, grid, biovizBase

 )

@@ -38,8 +38,6 @@ get.state.table <- function(hmm.list, numCPU=1) {

 	rownames(majorstate.f) <- rownames(majorstate)

 	## Plot the percentages of aneuploidies/states per chromosome

-	library(ggplot2)

-	library(reshape2)

 	df <- melt(majorstate.f, varnames=c("model","chromosome"), value.name="state")

 	df$state <- factor(df$state, levels=state.labels)

 	ggplt <- ggplot(df) + geom_bar(aes(x=chromosome, fill=state, y=..count..)) + theme_bw() + ylab('number of samples') + scale_fill_manual(values=state.colors)

@@ -109,7 +107,6 @@ get.reproducibility.one2many <- function(query, num.query.samples=1, subjects, n

 	## Overlap each subject's states with query

 	cat('overlapping each subject\'s states with query\n')

-	library(doParallel)

 	cl <- makeCluster(numCPU)

 	registerDoParallel(cl)

 	constates <- foreach (subject.gr = subjects.grl, .packages='GenomicRanges', .combine='cbind') %dopar% {

@@ -142,8 +139,6 @@ get.reproducibility.one2many <- function(query, num.query.samples=1, subjects, n

 		means[itest,] <- meanlevel

 	}

 	## Plot the results

-	library(ggplot2)

-	library(reshape2)

 	df <- melt(data.frame(state=mcols(query.gr)$state, meanstate=meanconstates))

 	ggplt <- ggplot(df) + geom_boxplot(aes(x=state, y=value, fill=variable)) + theme_bw() + coord_cartesian(ylim=c(1,length(levels(mcols(query.gr)$state)))) + scale_y_continuous(labels=levels(mcols(query.gr)$state)) + xlab(paste0(num.query.samples,' cell sample')) + ylab(paste0('single cell samples'))

@@ -214,8 +209,6 @@ get.reproducibility.many2many <- function(samples, num.tests=10, size.test=10, n

 	means <- cbind(as.data.frame(means), width=width(consensus))

 	## Plot the results

-	library(ggplot2)

-	library(reshape2)

 	states <- levels(mcols(samples.grl[[1]])$state)

 	ggplt1 <- ggplot(as.data.frame(means)) + geom_point(aes(x=V1, y=V2, alpha=width)) + theme_bw() + coord_cartesian(ylim=c(1,length(states))) + scale_x_continuous(breaks=1:length(states), labels=states) + scale_y_continuous(breaks=1:length(states), labels=states) + xlab(paste0(size.test,' single cell samples')) + ylab(paste0(size.test,' single cell samples'))

@@ -268,7 +268,6 @@ align2binned <- function(file, format, index=file, pairedEndReads=FALSE, chrom.l

 			### GC correction ###

 			if (GC.correction) {

-				library(Biostrings)

 				# Correct seqnames 1->chr1 if necessary

 				if (!grepl('chr',chromosome)) {

 					chrom <- paste0('chr',chromosome)

@@ -276,8 +275,8 @@ align2binned <- function(file, format, index=file, pairedEndReads=FALSE, chrom.l

 					chrom <- chromosome

 				}

 				## Calculating GC for whole bins

-				view.chr <- Views(GC.correction.bsgenome[[chrom]], ranges(i.binned.data))

-				freq <- alphabetFrequency(view.chr, as.prob = T, baseOnly=T)

+				view.chr <- Biostrings::Views(GC.correction.bsgenome[[chrom]], ranges(i.binned.data))

+				freq <- Biostrings::alphabetFrequency(view.chr, as.prob = T, baseOnly=T)

 				if (nrow(freq) > 1) {

 					GC.bin <- rowSums(freq[, c("G","C")])

 				} else {

@@ -14,7 +14,6 @@ ptm <- proc.time()

 	## Cluster the samples

 	cat("clustering the samples ...")

 	# Find states along the consensus template

-	library(foreach)

 	constates <- foreach (gr = grlred, .packages='GenomicRanges', .combine='cbind') %do% {

 		splt <- split(gr, mcols(gr)$state)

 		mind <- as.matrix(findOverlaps(consensus, splt))

@@ -6,7 +6,6 @@ hmmList2GRangesList <- function(hmm.list, reduce=TRUE, numCPU=1, consensus=FALSE

 	## Transform to GRanges

 	cat('transforming to GRanges ...')

 	if (numCPU > 1) {

-		suppressMessages( library(doParallel) )

 		cl <- makeCluster(numCPU)

 		registerDoParallel(cl)

 		cfun <- function(...) { GRangesList(...) }

@@ -53,7 +52,6 @@ hmmList2GRangesList <- function(hmm.list, reduce=TRUE, numCPU=1, consensus=FALSE

 binned2GRanges <- function(binned.data, chrom.length.file=NULL, offset=0) {

-	library(GenomicRanges)

 	gr <- GenomicRanges::GRanges(

 			seqnames = Rle(binned.data$chrom),

 			ranges = IRanges(start=binned.data$start+offset, end=binned.data$end+offset),

@@ -73,7 +71,6 @@ binned2GRanges <- function(binned.data, chrom.length.file=NULL, offset=0) {

 hmm2GRanges <- function(hmm, reduce=TRUE) {

-# 	library(GenomicRanges)

 	### Check user input ###

 	if (check.univariate.model(hmm)!=0 & check.multivariate.model(hmm)!=0) stop("argument 'hmm' expects a univariate or multivariate hmm object (type ?hmm for help)")

 	if (check.logical(reduce)!=0) stop("argument 'reduce' expects TRUE or FALSE")

@@ -87,7 +84,7 @@ hmm2GRanges <- function(hmm, reduce=TRUE) {

 			)

 	seqlengths(gr) <- hmm$seqlengths[names(seqlengths(gr))]

 	# Reorder seqlevels

-	gr <- GenomicRanges::keepSeqlevels(gr, names(hmm$seqlengths))

+	gr <- keepSeqlevels(gr, names(hmm$seqlengths))

 	if (reduce) {

 		# Reduce state by state

@@ -23,6 +23,9 @@ findCNVs <- function(binned.data, ID, method='univariate', eps=0.001, init="stan

 univariate.findCNVs <- function(binned.data, ID, eps=0.001, init="standard", max.time=-1, max.iter=-1, num.trials=1, eps.try=NULL, num.threads=1, read.cutoff.quantile=0.999, use.gc.corrected.reads=TRUE, strand='*') {

+	### Define cleanup behaviour ###

+	on.exit(.C("R_univariate_cleanup"))

+

 	## Intercept user input

 	IDcheck <- ID  #trigger error if not defined

 	if (class(binned.data) != 'GRanges') {

@@ -151,7 +154,6 @@ univariate.findCNVs <- function(binned.data, ID, eps=0.001, init="standard", max

 			loglik = double(length=1), # double* loglik

 			weights = double(length=numstates), # double* weights

 			distr.type = as.integer(state.distributions), # int* distr_type

-			ini.proc = as.integer(iniproc), # int* iniproc

 			size.initial = as.vector(size.initial), # double* initial_size

 			prob.initial = as.vector(prob.initial), # double* initial_prob

 			lambda.initial = as.vector(lambda.initial), # double* initial_lambda

@@ -211,7 +213,6 @@ univariate.findCNVs <- function(binned.data, ID, eps=0.001, init="standard", max

 			loglik = double(length=1), # double* loglik

 			weights = double(length=numstates), # double* weights

 			distr.type = as.integer(state.distributions), # int* distr_type

-			ini.proc = as.integer(iniproc), # int* iniproc

 			size.initial = as.vector(hmm$size), # double* initial_size

 			prob.initial = as.vector(hmm$prob), # double* initial_prob

 			lambda.initial = as.vector(hmm$lambda), # double* initial_lambda

@@ -497,6 +498,8 @@ bivariate.findCNVs <- function(binned.data, ID, eps=0.001, init="standard", max.

 		}

 		message(" done\n", appendLF=F)

+	### Define cleanup behaviour ###

+	on.exit(.C("R_multivariate_cleanup", as.integer(num.comb.states)))

 	### Run the multivariate HMM

 	# Call the C function

@@ -4,10 +4,6 @@

 # ============================================================

 plot.distribution <- function(model, state=NULL, chrom=NULL, start=NULL, end=NULL) {

-	## Load libraries

-	library(ggplot2)

-	library(reshape2)

-

 	# -----------------------------------------

 	# Get right x limit

 	get_rightxlim <- function(histdata, reads) {

@@ -158,8 +154,6 @@ plot.chromosomes.univariate <- function(model, file=NULL) {

 	}

 	## Setup page

-	library(grid)

-	library(ggplot2)

 	nrows <- 2	# rows for plotting chromosomes

 	ncols <- ceiling(num.chroms/nrows)

 	if (!is.null(file)) {

@@ -250,8 +244,6 @@ plot.chromosomes.bivariate <- function(model, file=NULL) {

 	custom.xlim <- model$distributions[[1]]['monosomy','mu'] * 10

 	## Setup page

-	library(grid)

-	library(ggplot2)

 	nrows <- 2	# rows for plotting chromosomes

 	ncols <- ceiling(num.chroms/nrows)

 	if (!is.null(file)) {

@@ -349,8 +341,6 @@ plot.genome.overview <- function(hmm.list, file, bp.per.cm=5e7, chromosome=NULL)

 	uni.hmm.grl <- lapply(hmm.list, '[[', 'bins')

 	## Setup page

-	library(grid)

-	library(ggplot2)

 	nrows <- length(uni.hmm.grl)	# rows for plotting genomes

 	ncols <- 1

 	total.length.bp <- sum(as.numeric(seqlengths(uni.hmm.grl[[1]])))

@@ -442,8 +432,6 @@ plot.genome.summary <- function(hmm.list, file='aneufinder_genome_overview') {

 	flattened_gr_srt <- sort(flattened_gr)

 	## Setup page

-	library(grid)

-	library(ggplot2)

 	nrows <- length(levels(uni.hmm.grl[[1]]$state))	# rows for plotting genomes

 	ncols <- 1

 	pdf(file=paste0(file, '.pdf'), width=ncols*24, height=nrows*2)

@@ -497,7 +485,7 @@ plot.genome.summary <- function(hmm.list, file='aneufinder_genome_overview') {

 			strand = Rle(strand("*")), cov=cov

 		)

-		trans_gr <- transformToGenome(gr, space.skip = 0)

+		trans_gr <- biovizBase::transformToGenome(gr, space.skip = 0)

 		dfplot <- as.data.frame(trans_gr)

deleted file mode 100644

@@ -1,138 +0,0 @@

-simulate.univariate = function(coordinates, transition, emission, initial=1) {

-

-	# Calculate some variables

-	numstates = ncol(transition)

-	if (numstates!=3) {

-		stop("The transition matrix is expected to have 3 columns and 3 rows")

-	}

-	numbins = nrow(coordinates)

-

-	## Make state vector from transition matrix

-	# Generate sample of random numbers

-	rsample = runif(numbins,0,1)

-	# Integrate transition matrix by row and add -1 in front

-	cumtransition = cbind(rep(-1,numstates), t(apply(transition, 1, cumsum)))

-	# Generate the state vector by going through each state

-	cat("Generating states from transition matrix...")

-	states = matrix(rep(NA,numstates*numbins), ncol=numstates)

-	for (irow in 1:numstates) {

-		states[,irow] = findInterval(rsample, cumtransition[irow,], rightmost.closed=TRUE)

-	}

-	statevec = rep(NA,numbins)

-	statevec[1] = initial

-	for (ibin in 2:numbins) {

-		statevec[ibin] = states[ibin,statevec[ibin-1]]

-	}

-	cat(" done\n")

-

-	## Make reads from state vector and distributions

-	# Generate the read counts by drawing from distribution

-	cat("Generating reads from emission parameters and states...")

-	reads = rep(NA, numbins)

-	numbins.in.state = aggregate(rep(1,length(statevec)), list(state=statevec), sum)

-	reads[statevec==1] = 0

-	if (!is.na(numbins.in.state[2,'x'])) {

-		reads[statevec==2] = rnbinom(numbins.in.state[2,'x'], size=emission[2,'r'], prob=emission[2,'p'])

-	}

-	if (!is.na(numbins.in.state[3,'x'])) {

-		reads[statevec==3] = rnbinom(numbins.in.state[3,'x'], size=emission[3,'r'], prob=emission[3,'p'])

-	}

-	cat(" done\n")

-

-	# Return the output

-	out = list(coordinates = coordinates,

-				states = statevec,

-				reads = reads,

-				transition = transition,

-				emission = emission

-				)

-	return(out)

-

-}

-

-

-

-simulate.multivariate = function(coordinates, transition, emissions, weights, sigma, use.states, initial=1) {

-

-	lib = require(mvtnorm)

-	if (lib == FALSE) {

-		install.packages("mvtnorm")

-		library(mvtnorm)

-	}

-

-	# Calculate some variables

-	numstates = ncol(transition)

-	numbins = nrow(coordinates)

-	nummod = length(emissions)

-

-	## Make state vector from transition matrix

-	# Generate sample of random numbers

-	rsample = runif(numbins,0,1)

-	# Integrate transition matrix by row and add -1 in front

-	cumtransition = cbind(rep(-1,numstates), t(apply(transition, 1, cumsum)))

-	# Generate the state vector by going through each state

-	cat("Generating states from transition matrix...")

-	states = matrix(rep(NA,numstates*numbins), ncol=numstates)

-	for (irow in 1:numstates) {

-		states[,irow] = findInterval(rsample, cumtransition[irow,], rightmost.closed=TRUE)

-	}

-	statevec = rep(NA,numbins)

-	statevec[1] = initial

-	for (ibin in 2:numbins) {

-		statevec[ibin] = states[ibin,statevec[ibin-1]]

-	}

-	# Replace the states by combinatorial states

-	statevec = use.states[statevec]

-	cat(" done\n")

-

-	## Make reads from state vector and emission distributions

-	reads = matrix(rep(NA, nummod*numbins), ncol=nummod)

-

-	rs = unlist(lapply(emissions, "[", 2:3, 'r'))

-	ps = unlist(lapply(emissions, "[", 2:3, 'p'))

-	ws = unlist(lapply(weights, "[", 1))

-	for (istate in use.states) {

-		cat("Generating reads for state ",istate,"\n")

-		i = which(use.states==istate)

-		

-		# Convert istate to binary representation

-		binary_state = rev(as.integer(intToBits(istate))[1:nummod])

-		binary_stateTF = unlist(list(c(TRUE,FALSE), c(FALSE,TRUE))[binary_state+1])

-

-		# Draw from the multivariate normal

-		n = length(which(statevec==istate))

-		if (n == 0) next

-		cat("drawing from multivariate normal...             \r")

-		z = matrix(rmvnorm(n, mean=rep(0,nummod), sigma=sigma[,,i]), ncol=nummod)

-		# Transform to uniform space

-		cat("transforming to uniform space...                \r")

-		u = matrix(apply(z, 2, pnorm), ncol=nummod)

-		# Transform to count space using marginals

-		cat("transform to count space...                     \r")

-		irs = rs[binary_stateTF]

-		ips = ps[binary_stateTF]

-		for (imod in 1:nummod) {

-			mask = statevec==istate

-			if (binary_state[imod] == 0) {

-				reads[mask,imod] = qzinbinom(u[,imod], w=ws[imod], size=irs[imod], prob=ips[imod])

-			} else {

-				reads[mask,imod] = qnbinom(u[,imod], size=irs[imod], prob=ips[imod])

-			}

-		}

-		cat("                                                \r")

-			

-	}

-

-	# Return the output

-	out = list(coordinates = coordinates,

-				states = statevec,

-				reads = reads,

-				emissions = emissions,

-				transition = transition,

-				sigma = sigma,

-				use.states = use.states,

-				weights = weights

-				)

-	return(out)

-

-}

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

 #include "utility.h"

 #include "scalehmm.h"

 #include "loghmm.h"

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

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

+#include <string> // strcmp

+

+static ScaleHMM* hmm; // declare as static outside the function because we only need one and this enables memory-cleanup on R_CheckUserInterrupt()

+static double** multiD;

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

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

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

 extern "C" {

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

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

 {

 	// Define logging level

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

 // 	Output2FILE::Stream() = pFile;

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

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

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

 // 	// Parallelization settings

-// 	omp_set_num_threads(*num_threads);

+	omp_set_num_threads(*num_threads);

 	// Print some information

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

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

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

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

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

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

 	if (*maxiter < 0)

 	{

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

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

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

 	} else {

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

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

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

 	}

 	if (*maxtime < 0)

 	{

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

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

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

 	} else {

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

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

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

 	}

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

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

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

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

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

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

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

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

 	}

 	// Flush Rprintf statements to console

 	R_FlushConsole();

 	// Create the HMM

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

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

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

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

 // 	LogHMM* hmm = new LogHMM(*T, *N);

 	hmm->set_cutoff(*read_cutoff);

 	// Initialize the transition probabilities and proba

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

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

 	}

 	variance = variance / *T;

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

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

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

 	// Create the emission densities and initialize

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

 	{

 		if (distr_type[i_state] == 1)

 		{

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

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

 			ZeroInflation *d = new ZeroInflation(O, *T); // delete is done inside ~ScaleHMM()

 			hmm->densityFunctions.push_back(d);

 		}

 		else if (distr_type[i_state] == 2)

 		{

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

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

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

 			hmm->densityFunctions.push_back(d);

 		}

 		else if (distr_type[i_state] == 3)

 		{

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

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

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

 			hmm->densityFunctions.push_back(d);

 		}

 		else if (distr_type[i_state] == 4)

 		{

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

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

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

 			hmm->densityFunctions.push_back(d);

 		}

 		else if (distr_type[i_state] == 5)

 		{

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

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

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

 			hmm->densityFunctions.push_back(d);

 		}

 		else

 		{

-			FILE_LOG(logWARNING) << "Density not specified, using default negative binomial for state " << i_state;

+			//FILE_LOG(logWARNING) << "Density not specified, using default negative binomial for state " << i_state;

 			NegativeBinomial *d = new NegativeBinomial(O, *T, initial_size[i_state], initial_prob[i_state]);

 			hmm->densityFunctions.push_back(d);

 		}

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

 	R_FlushConsole();

 	// Do the Baum-Welch to estimate the parameters

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

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

 	try

 	{

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

 	}

 	catch (std::exception& e)

 	{

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

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

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

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

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

 		else { *error = 2; }

 	}

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

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

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

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

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

 // 	#pragma omp parallel for

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

 // 	{

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

 // 	}

 	// Compute the states from posteriors

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

-	double posterior_per_t [*N];

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

+	int ind_max;

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

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

 	{

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

 		{

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

 		}

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

+		ind_max = std::distance(posterior_per_t.begin(), std::max_element(posterior_per_t.begin(), posterior_per_t.end()));

+		states[t] = state_labels[ind_max];

 	}

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

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

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

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

 	{

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

 	*loglik = hmm->get_logP();

 	hmm->calc_weights(weights);

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

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

 	delete hmm;

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

 }

 } // extern C

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

 // 	Output2FILE::Stream() = pFile;

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

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

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

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

 	// Parallelization settings

-// 	omp_set_num_threads(*num_threads);

+	omp_set_num_threads(*num_threads);

 	// Print some information

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

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

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

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

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

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

 	if (*maxiter < 0)

 	{

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

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

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

 	} else {

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

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

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

 	}

 	if (*maxtime < 0)

 	{

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

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

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

 	} else {

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

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

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

 	}

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

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

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

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

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

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

 	// Flush Rprintf statements to console

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

 	// Recode the densities vector to matrix representation

 // 	clock_t clocktime = clock(), dtime;

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

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

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

 	{

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

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

 		}

 	}

 // 	dtime = clock() - clocktime;

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

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

 	// Create the HMM

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

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

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

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

 	// Initialize the transition probabilities and proba

 	hmm->initialize_transition_probs(initial_A, *use_initial_params);

 	hmm->initialize_proba(initial_proba, *use_initial_params);

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

 	// Print logproba and A

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

 // 	{

-// 		FILE_LOG(logDEBUG) << "proba["<<iN<<"] = " <<exp(hmm->logproba[iN]);

+// 		//FILE_LOG(logDEBUG) << "proba["<<iN<<"] = " <<exp(hmm->logproba[iN]);

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

 // 		{

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

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

 // 		}

 // 	}

 	// Estimate the parameters

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

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

 	try

 	{

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

 	}

 	catch (std::exception& e)

 	{

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

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

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

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

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

 		else { *error = 2; }

 	}

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

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

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

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

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

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

 // 	{

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

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

 // 	}

 	// Compute the states from posteriors

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

-	double posterior_per_t [*N];

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

+	int ind_max;

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

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

 	{

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

 		{

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

 		}

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

+		ind_max = std::distance(posterior_per_t.begin(), std::max_element(posterior_per_t.begin(), posterior_per_t.end()));

+		states[t] = comb_states[ind_max];

 	}

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

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

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

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

 	{

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

 	}

 	*loglik = hmm->get_logP();

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

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

 	delete hmm;

-	freeDoubleMatrix(multiD, *N);

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

+// 	FreeDoubleMatrix(multiD, *N);

 }

 } // extern C

+

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

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

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

+extern "C" {

+void R_univariate_cleanup()

+{

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

+	delete hmm;

+}

+}

+

+extern "C" {

+void R_multivariate_cleanup(int* N)

+{

+	delete hmm;

+	FreeDoubleMatrix(multiD, *N);

+}

+}

+

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

 // Constructor and Destructor ---------------------------------

 Normal::Normal(int* observations, int T, double mean, double variance)

 {

-	FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;

+	//FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;

 	this->obs = observations;

 	this->T = T;

 	this->mean = mean;

@@ -17,13 +17,13 @@ Normal::Normal(int* observations, int T, double mean, double variance)

 Normal::~Normal()

 {

-	FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;

+	//FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;

 }

 // Methods ----------------------------------------------------

 void Normal::calc_densities(double* density)

 {

-	FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;

+	//FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;

 	for (int t=0; t<this->T; t++)

 	{

 		density[t] = dnorm(this->obs[t], this->mean, this->sd, 0);

@@ -32,7 +32,7 @@ void Normal::calc_densities(double* density)

 void Normal::calc_logdensities(double* logdensity)

 {

-	FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;

+	//FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;

 	for (int t=0; t<this->T; t++)

 	{

 		logdensity[t] = dnorm(this->obs[t], this->mean, this->sd, 1);

@@ -41,51 +41,52 @@ void Normal::calc_logdensities(double* logdensity)

 void Normal::update(double* weights)

 {

-	FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;

+	//FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;

+	(void)weights;

 // TODO

 }

 // Getter and Setter ------------------------------------------

 DensityName Normal::get_name()

 {

-	FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;

+	//FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;

 	return(NORMAL);

 }

 double Normal::get_mean()

 {

-	FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;

+	//FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;

 	return(this->mean);

 }

 void Normal::set_mean(double mean)

 {

-	FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;

+	//FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;

 	this->mean = mean;

 }

 double Normal::get_variance()

 {

-	FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;

+	//FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;

 	return(this->variance);

 }

 void Normal::set_variance(double variance)

 {

-	FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;

+	//FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;

 	this->variance = variance;

 	this->sd = sqrt(variance);

 }

 double Normal::get_stdev()

 {

-	FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;

+	//FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;

 	return(this->sd);

 }

 void Normal::set_stdev(double stdev)

 {

-	FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;

+	//FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;

 	this->sd = stdev;

 	this->variance = stdev*stdev;

 }

@@ -98,7 +99,7 @@ void Normal::set_stdev(double stdev)

 // Constructor and Destructor ---------------------------------

 Poisson::Poisson(int* observations, int T, double lambda)

 {

-	FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;

+	//FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;

 	this->obs = observations;

 	this->T = T;

 	this->lambda = lambda;

@@ -107,8 +108,8 @@ Poisson::Poisson(int* observations, int T, double lambda)

 	if (this->obs != NULL)

 	{

 		this->max_obs = intMax(observations, T);

-		this->lxfactorials = (double*) calloc(max_obs+1, sizeof(double));

-		this->lxfactorials[0] = 0.0;	// Not necessary, already 0 because of calloc

+		this->lxfactorials = (double*) Calloc(max_obs+1, double);

+		this->lxfactorials[0] = 0.0;	// Not necessary, already 0 because of Calloc

 		this->lxfactorials[1] = 0.0;

 		for (int j=2; j<=max_obs; j++)

 		{

@@ -119,25 +120,25 @@ Poisson::Poisson(int* observations, int T, double lambda)

 Poisson::~Poisson()

 {

-	FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;

+	//FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;

 	if (this->lxfactorials != NULL)

 	{

-		free(this->lxfactorials);

+		Free(this->lxfactorials);

 	}

 }

 // Methods ----------------------------------------------------

 void Poisson::calc_logdensities(double* logdens)

 {

-	FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;

+	//FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;

 	double logl = log(this->lambda);

 	double l = this->lambda;

 	double lxfactorial;

 	// Select strategy for computing densities

 	if (this->max_obs <= this->T)

 	{

-		FILE_LOG(logDEBUG2) << "Precomputing densities in " << __func__ << " for every obs[t], because max(O)<=T";

-		double logdens_per_read [this->max_obs+1];

+		//FILE_LOG(logDEBUG2) << "Precomputing densities in " << __func__ << " for every obs[t], because max(O)<=T";

+		std::vector<double> logdens_per_read(this->max_obs+1);

 		for (int j=0; j<=this->max_obs; j++)

 		{

 			logdens_per_read[j] = j*logl - l - this->lxfactorials[j];

@@ -145,27 +146,27 @@ void Poisson::calc_logdensities(double* logdens)

 		for (int t=0; t<this->T; t++)

 		{

 			logdens[t] = logdens_per_read[(int) this->obs[t]];

-			FILE_LOG(logDEBUG4) << "logdens["<<t<<"] = " << logdens[t];

+			//FILE_LOG(logDEBUG4) << "logdens["<<t<<"] = " << logdens[t];

 			if (isnan(logdens[t]))

 			{

-				FILE_LOG(logERROR) << __PRETTY_FUNCTION__;

-				FILE_LOG(logERROR) << "logdens["<<t<<"] = "<< logdens[t];

+				//FILE_LOG(logERROR) << __PRETTY_FUNCTION__;

+				//FILE_LOG(logERROR) << "logdens["<<t<<"] = "<< logdens[t];

 				throw nan_detected;

 			}

 		}

 	}

 	else

 	{

-		FILE_LOG(logDEBUG2) << "Computing densities in " << __func__ << " for every t, because max(O)>T";

+		//FILE_LOG(logDEBUG2) << "Computing densities in " << __func__ << " for every t, because max(O)>T";

 		for (int t=0; t<this->T; t++)

 		{

 			lxfactorial = this->lxfactorials[(int) this->obs[t]];

 			logdens[t] = this->obs[t]*logl - l - lxfactorial;

-			FILE_LOG(logDEBUG4) << "logdens["<<t<<"] = " << logdens[t];

+			//FILE_LOG(logDEBUG4) << "logdens["<<t<<"] = " << logdens[t];

 			if (isnan(logdens[t]))

 			{