@@ -14,7 +14,6 @@ plot.distribution,

 univariate2bed,

 binned2wiggle,

 multivariate2bed,

-binned2GRanges,

 hmm2GRanges,

 mixedsort,

@@ -22,6 +21,6 @@ mixedorder

 )

 # Import all packages listed as Imports or Depends

-# import(

-# 	GenomicRanges, IRanges

-# )

+import(

+	GenomicRanges, IRanges

+)

@@ -20,23 +20,25 @@ align2binned <- function(file, format, index=file, chrom.length.file, outputfold

 		dir.create(outputfolder)

 	}

-	## Read in the data

+	### Read in the data

+	## BED (0-based)

 	if (format == "bed") {

 		cat("Reading file",basename(file),"...")

-		# File with chromosome lengths

+		# File with chromosome lengths (1-based)

 		chrom.lengths.df <- read.table(chrom.length.file)

 		chrom.lengths <- chrom.lengths.df[,2]

 		names(chrom.lengths) <- chrom.lengths.df[,1]

-		# File with reads, determine classes first for faster import

+		# File with reads, determine classes first for faster import (0-based)

 		tab5rows <- read.table(file, nrows=5)

 		classes.in.bed <- sapply(tab5rows, class)

 		classes <- rep("NULL",length(classes.in.bed))

 		classes[1:3] <- classes.in.bed[1:3]

 		data <- read.table(file, colClasses=classes)

 		# Convert to GRanges object

-		data <- GenomicRanges::GRanges(seqnames=Rle(data[,1]), ranges=IRanges(start=data[,2], end=data[,3]), strand=Rle(strand("*"), nrow(data)))

+		data <- GenomicRanges::GRanges(seqnames=Rle(data[,1]), ranges=IRanges(start=data[,2]+1, end=data[,3]+1), strand=Rle(strand("*"), nrow(data)))	# +1 to match coordinate systems

 		seqlengths(data) <- as.integer(chrom.lengths[names(seqlengths(data))])

 		chroms.in.data <- seqlevels(data)

+	## BAM (1-based)

 	} else if (format == "bam") {

 		library(Rsamtools) # TODO: put this in Imports in finished package

 		library(GenomicAlignments) # TODO: put this in Imports in finished package

@@ -44,9 +46,10 @@ align2binned <- function(file, format, index=file, chrom.length.file, outputfold

 		file.header <- Rsamtools::scanBamHeader(file)[[1]]

 		chrom.lengths <- file.header$targets

 		chroms.in.data <- names(chrom.lengths)

+	## BEDGraph (0-based)

 	} else if (format == "bedGraph") {

 		cat("Reading file",basename(file),"...")

-		# File with chromosome lengths

+		# File with chromosome lengths (1-based)

 		chrom.lengths.df <- read.table(chrom.length.file)

 		chrom.lengths <- chrom.lengths.df[,2]

 		names(chrom.lengths) <- chrom.lengths.df[,1]

@@ -57,7 +60,7 @@ align2binned <- function(file, format, index=file, chrom.length.file, outputfold

 		classes[1:4] <- classes.in.bed[1:4]

 		data <- read.table(file, colClasses=classes)

 		# Convert to GRanges object

-		data <- GenomicRanges::GRanges(seqnames=Rle(data[,1]), ranges=IRanges(start=data[,2], end=data[,3]), strand=Rle(strand("*"), nrow(data)), signal=data[,4])

+		data <- GenomicRanges::GRanges(seqnames=Rle(data[,1]), ranges=IRanges(start=data[,2]+1, end=data[,3]+1), strand=Rle(strand("*"), nrow(data)), signal=data[,4])	# +1 to match coordinate systems

 		seqlengths(data) <- as.integer(chrom.lengths[names(seqlengths(data))])

 		chroms.in.data <- seqlevels(data)

 	}

@@ -69,7 +72,7 @@ align2binned <- function(file, format, index=file, chrom.length.file, outputfold

 		cat("Binning into binsize",binsize,"\n")

 		### Iterate over all chromosomes

-		binned.data.allchroms <- NULL

+		binned.data <- GenomicRanges::GRanges()

 		if (is.null(chromosomes)) {

 			chromosomes <- chroms.in.data

 		}

@@ -86,38 +89,43 @@ align2binned <- function(file, format, index=file, chrom.length.file, outputfold

 			## Check last incomplete bin

 			incomplete.bin <- chrom.lengths[chromosome] %% binsize > 0

 			if (incomplete.bin) {

-				numbins <- ceiling(chrom.lengths[chromosome]/binsize)

+				numbins <- floor(chrom.lengths[chromosome]/binsize)	# floor: we don't want incomplete bins, ceiling: we want incomplete bins at the end

 			} else {

 				numbins <- chrom.lengths[chromosome]/binsize

 			}

+			if (numbins == 0) {

+				warning("chromosome ",chromosome," is smaller than binsize. Skipped.")

+				next

+			}

 			## Initialize vectors

 			cat("initialize vectors...\r")

 			chroms <- rep(chromosome,numbins)

 			reads <- rep(0,numbins)

-			start <- seq(from=0, by=binsize, length.out=numbins)

-			end <- seq(from=binsize-1, by=binsize, length.out=numbins)

-			end[length(end)] <- chrom.lengths[chromosome]

+			start <- seq(from=1, by=binsize, length.out=numbins)

+			end <- seq(from=binsize+1, by=binsize, length.out=numbins)

+# 			end[length(end)] <- chrom.lengths[chromosome] # last ending coordinate is size of chromosome, only if incomplete bins are desired

 			## Create binned chromosome as GRanges object

 			cat("creating GRanges container...            \r")

-			ichrom <- GenomicRanges::GRanges(seqnames = Rle(chromosome, numbins),

+			i.binned.data <- GenomicRanges::GRanges(seqnames = Rle(chromosome, numbins),

 							ranges = IRanges(start=start, end=end),

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

 							)

+			seqlengths(i.binned.data) <- chrom.lengths[chromosome]

 			if (format=="bam") {

 				cat("reading reads from file...               \r")

-				data <- GenomicAlignments::readGAlignmentsFromBam(file, index=index, param=ScanBamParam(what=c("pos"),which=range(ichrom)))

+				data <- GenomicAlignments::readGAlignmentsFromBam(file, index=index, param=ScanBamParam(what=c("pos"),which=range(i.binned.data)))

 			}

 			## Count overlaps

 			cat("counting overlaps...                     \r")

 			if (format=="bam" | format=="bed") {

-				reads <- GenomicRanges::countOverlaps(ichrom, data[seqnames(data)==chromosome])

+				reads <- GenomicRanges::countOverlaps(i.binned.data, data[seqnames(data)==chromosome])

 			} else if (format=="bedGraph") {

 				# Take the max value from all regions that fall into / overlap a given bin as read count

-				midx <- as.matrix(findOverlaps(ichrom, data[seqnames(data)==chromosome]))

-				reads <- rep(0,length(ichrom))

+				midx <- as.matrix(findOverlaps(i.binned.data, data[seqnames(data)==chromosome]))

+				reads <- rep(0,length(i.binned.data))

 				signal <- mcols(data)$signal

 				rle <- rle(midx[,1])

 				read.idx <- rle$values

@@ -132,10 +140,10 @@ align2binned <- function(file, format, index=file, chrom.length.file, outputfold

 			## Concatenate

 			cat("concatenate...                           \r")

-			binned.data <- data.frame(chroms,start,end,reads)

-			names(binned.data) <- binned.data.names

+			mcols(i.binned.data)$reads <- reads

 			if (separate.chroms==TRUE) {

+				binned.data <- i.binned.data

 				if (save.as.RData==TRUE) {

 					## Print to file

 					filename <- paste(basename(file),"_binsize_",binsize,"_",chromosome,".RData", sep="")

@@ -146,23 +154,20 @@ align2binned <- function(file, format, index=file, chrom.length.file, outputfold

 					return(binned.data)

 				}

 			} else {

-				binned.data.allchroms[[length(binned.data.allchroms)+1]] <- binned.data

+				binned.data <- suppressWarnings(BiocGenerics::append(binned.data, i.binned.data))

 			}

 			cat("                                         \r")

 		}

-		if (separate.chroms!=TRUE) {

-			cat("Concatenating chromosomes ...")

-			binned.data.allchroms <- do.call("rbind",binned.data.allchroms)

-			cat(" done\n")

+		if (separate.chroms==FALSE) {

 			if (save.as.RData==TRUE) {

 				# Print to file

 				filename <- paste(basename(file),"_binsize_",binsize,".RData", sep="")

 				cat("Saving to file ...")

-				save(binned.data.allchroms, file=file.path(outputfolder,filename) )

+				save(binned.data, file=file.path(outputfolder,filename) )

 				cat(" done\n")

 			} else {

-				return(binned.data.allchroms)

+				return(binned.data)

 			}

 		}

@@ -1,19 +1,6 @@

-binned2GRanges <- function(binned.data) {

-

-	library(GenomicRanges)

-	gr <- GenomicRanges::GRanges(

-			seqnames = Rle(binned.data$chrom),

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

-			strand = Rle(strand("*"), nrow(binned.data)),

-			reads = binned.data$reads

-			)

-	return(gr)

-

-}

-

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

-	library(GenomicRanges)

+# 	library(GenomicRanges)

 	### Check user input ###

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

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

@@ -25,10 +12,9 @@ hmm2GRanges <- function(hmm, reduce=TRUE) {

 			ranges = IRanges(start=hmm$coordinates$start, end=hmm$coordinates$end),

 			strand = Rle(strand("*"), nrow(hmm$coordinates))

 			)

-	# Seqlengths gives warnings because our coordinates are 0-based

-# 	for (chrom in unique(hmm$coordinates$chrom)) {

-# 		seqlengths(gr)[names(seqlengths(gr))==chrom] <- max(hmm$coordinates$end[hmm$coordinates$chrom==chrom])

-# 	}

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

+	# Reorder seqlevels

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

 	if (reduce) {

 		# Reduce state by state

@@ -16,25 +16,24 @@ find.aneuploidies <- function(binned.data, use.states=0:3, eps=0.001, init="stan

 	war <- NULL

 	if (is.null(eps.try)) eps.try <- eps

-	names(binned.data) <- binned.data.names # defined globally outside this function

-

 	## Assign variables

 # 	state.labels # assigned globally outside this function

 	use.state.labels <- state.labels[use.states+1]

 	numstates <- length(use.states)

-	numbins <- length(binned.data$reads)

+	numbins <- length(binned.data)

+	reads <- mcols(binned.data)$reads

 	iniproc <- which(init==c("standard","random","empiric")) # transform to int

 	# Check if there are reads in the data, otherwise HMM will blow up

-	if (!any(binned.data$reads!=0)) {

+	if (!any(reads!=0)) {

 		stop("All reads in data are zero. No univariate HMM done.")

 	}

 	# Filter high reads out, makes HMM faster

-	read.cutoff <- as.integer(quantile(binned.data$reads, 0.9999))

+	read.cutoff <- as.integer(quantile(reads, 0.9999))

 	if (filter.reads) {

-		mask <- binned.data$reads > read.cutoff

-		binned.data$reads[mask] <- read.cutoff

+		mask <- reads > read.cutoff

+		reads[mask] <- read.cutoff

 		numfiltered <- length(which(mask))

 		if (numfiltered > 0) {

 			warning(paste("There are very high read counts in your data (probably artificial). Replaced read counts > ",read.cutoff," (99.99% quantile) by ",read.cutoff," in ",numfiltered," bins. Set option 'filter.reads=FALSE' to disable this filtering.", sep=""))

@@ -47,16 +46,16 @@ find.aneuploidies <- function(binned.data, use.states=0:3, eps=0.001, init="stan

 	for (i_try in 1:num.trials) {

 		cat("\n\nTry ",i_try," of ",num.trials," ------------------------------\n")

 		hmm <- .C("R_univariate_hmm",

-			reads = as.integer(binned.data$reads), # int* O

+			reads = as.integer(reads), # int* O

 			num.bins = as.integer(numbins), # int* T

 			num.states = as.integer(numstates), # int* N

-			use.states = as.integer(use.states), # int* states

+			use.states = as.integer(use.states), # int* statelabels

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

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

 			num.iterations = as.integer(max.iter), #  int* maxiter

 			time.sec = as.integer(max.time), # double* maxtime

 			loglik.delta = as.double(eps.try), # double* eps

-			posteriors = double(length=numbins * numstates), # double* posteriors

+			states = integer(length=numbins), # int* states

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

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

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

@@ -89,7 +88,6 @@ find.aneuploidies <- function(binned.data, use.states=0:3, eps=0.001, init="stan

 				warning("HMM did not converge in trial run ",i_try,"!\n")

 			}

 			# Store model in list

-			hmm$posteriors <- NULL

 			hmm$reads <- NULL

 			modellist[[i_try]] <- hmm

 		}

@@ -104,16 +102,16 @@ find.aneuploidies <- function(binned.data, use.states=0:3, eps=0.001, init="stan

 		# Rerun the HMM with different epsilon and initial parameters from trial run

 		cat("\n\nRerunning try ",indexmax," with eps =",eps,"--------------------\n")

 		hmm <- .C("R_univariate_hmm",

-			reads = as.integer(binned.data$reads), # int* O

+			reads = as.integer(reads), # int* O

 			num.bins = as.integer(numbins), # int* T

 			num.states = as.integer(numstates), # int* N

-			use.states = as.integer(use.states), # int* states

+			use.states = as.integer(use.states), # int* statelabels

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

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

 			num.iterations = as.integer(max.iter), #  int* maxiter

 			time.sec = as.integer(max.time), # double* maxtime

 			loglik.delta = as.double(eps), # double* eps

-			posteriors = double(length=numbins * numstates), # double* posteriors

+			states = integer(length=numbins), # int* states

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

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

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

@@ -132,11 +130,11 @@ find.aneuploidies <- function(binned.data, use.states=0:3, eps=0.001, init="stan

 	# Add useful entries

 	names(hmm$weights) <- use.state.labels

-	hmm$coordinates <- binned.data[,coordinate.names]

-	hmm$posteriors <- matrix(hmm$posteriors, ncol=hmm$num.states)

-	colnames(hmm$posteriors) <- paste("P(",use.state.labels,")", sep="")

+	hmm$coordinates <- data.frame(as.character(seqnames(binned.data)), start(ranges(binned.data)), end(ranges(binned.data)))

+	names(hmm$coordinates) <- coordinate.names

+	hmm$seqlengths <- seqlengths(binned.data)

 	class(hmm) <- class.aneufinder.univariate

-	hmm$states <- use.state.labels[apply(hmm$posteriors, 1, which.max)]

+	hmm$states <- use.state.labels[hmm$states+1]

 	hmm$eps <- eps

 	hmm$A <- matrix(hmm$A, ncol=hmm$num.states, byrow=TRUE)

 	rownames(hmm$A) <- use.state.labels

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

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

 extern "C" {

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

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

 {

 	// Define logging level

@@ -89,8 +89,8 @@ void R_univariate_hmm(int* O, int* T, int* N, int* states, double* size, double*

 			if (*iniproc == 1)

 			{

 				// Simple initialization based on data mean, assumed to be the disomic mean

-				imean = mean/2 * states[i_state];

-				ivariance = variance/2 * states[i_state];

+				imean = mean/2 * statelabels[i_state];

+				ivariance = variance/2 * statelabels[i_state];

 			}

 			// Calculate r and p from mean and variance

@@ -137,15 +137,28 @@ void R_univariate_hmm(int* O, int* T, int* N, int* states, double* size, double*

 	}

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

-	#pragma omp parallel for

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

+

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

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

+// 	#pragma omp parallel for

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

+// 	{

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

+// 		{

+// 			posteriors[t + iN * (*T)] = hmm->get_posterior(iN, t);

+// 		}

+// 	}

+

+	// Compute the states from posteriors

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

+	double posterior_per_t [*N];

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

 	{

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

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

 		{

-			posteriors[t + iN * (*T)] = hmm->get_posterior(iN, t);

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

 		}

+		states[t] = statelabels[argMax(posterior_per_t, *N)];

 	}

 	FILE_LOG(logDEBUG1) << "Return parameters";

@@ -17,7 +17,7 @@ ScaleHMM::ScaleHMM(int T, int N)

 	this->scalealpha = allocDoubleMatrix(T, N);

 	this->scalebeta = allocDoubleMatrix(T, N);

 	this->densities = allocDoubleMatrix(N, T);

-	this->tdensities = allocDoubleMatrix(T, N);

+// 	this->tdensities = allocDoubleMatrix(T, N);

 	this->proba = (double*) calloc(N, sizeof(double));

 	this->gamma = allocDoubleMatrix(N, T);

 	this->sumgamma = (double*) calloc(N, sizeof(double));

@@ -26,7 +26,7 @@ ScaleHMM::ScaleHMM(int T, int N)

 	this->dlogP = INFINITY;

 	this->sumdiff_state_last = 0;

 	this->sumdiff_posterior = 0.0;

-	this->use_tdens = false;

+// 	this->use_tdens = false;

 }

@@ -38,7 +38,7 @@ ScaleHMM::~ScaleHMM()

 	freeDoubleMatrix(this->scalealpha, this->T);

 	freeDoubleMatrix(this->scalebeta, this->T);

 	freeDoubleMatrix(this->densities, this->N);

-	freeDoubleMatrix(this->tdensities, this->T);

+// 	freeDoubleMatrix(this->tdensities, this->T);

 	freeDoubleMatrix(this->gamma, this->N);

 	freeDoubleMatrix(this->sumxi, this->N);

 	free(this->proba);

@@ -213,7 +213,7 @@ void ScaleHMM::baumWelch(int* maxiter, int* maxtime, double* eps)

 		this->print_uni_iteration(iteration);

 		// Check convergence

-		if(fabs(this->dlogP) < *eps) //it has converged

+		if(this->dlogP < *eps) //it has converged

 		{

 			FILE_LOG(logINFO) << "\nConvergence reached!\n";

 			Rprintf("\nConvergence reached!\n\n");

@@ -371,8 +371,8 @@ void ScaleHMM::forward()

 	FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;

 	clock_t time = clock(), dtime;

-	if (not this->use_tdens)

-	{

+// 	if (not this->use_tdens)

+// 	{

 		double alpha [this->N];

 		// Initialization

@@ -424,61 +424,61 @@ void ScaleHMM::forward()

 			}

 		}

-	}

-	else if (this->use_tdens)

-	{

-

-		double alpha [this->N];

-		// Initialization

-		this->scalefactoralpha[0] = 0.0;

-		for (int iN=0; iN<this->N; iN++)

-		{

-			alpha[iN] = this->proba[iN] * this->tdensities[0][iN];

-			FILE_LOG(logDEBUG4) << "alpha["<<iN<<"] = " << alpha[iN];

-			this->scalefactoralpha[0] += alpha[iN];

-		}

-		FILE_LOG(logDEBUG4) << "scalefactoralpha["<<0<<"] = " << scalefactoralpha[0];

-		for (int iN=0; iN<this->N; iN++)

-		{

-			this->scalealpha[0][iN] = alpha[iN] / this->scalefactoralpha[0];

-			FILE_LOG(logDEBUG4) << "scalealpha["<<0<<"]["<<iN<<"] = " << scalealpha[0][iN];

-		}

-		// Induction

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

-		{

-			this->scalefactoralpha[t] = 0.0;

-			for (int iN=0; iN<this->N; iN++)

-			{

-				double helpsum = 0.0;

-				for (int jN=0; jN<this->N; jN++)

-				{

-					helpsum += this->scalealpha[t-1][jN] * this->A[jN][iN];

-				}

-				alpha[iN] = helpsum * this->tdensities[t][iN];

-				FILE_LOG(logDEBUG4) << "alpha["<<iN<<"] = " << alpha[iN];

-				this->scalefactoralpha[t] += alpha[iN];

-			}

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

-			for (int iN=0; iN<this->N; iN++)

-			{

-				this->scalealpha[t][iN] = alpha[iN] / this->scalefactoralpha[t];

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

-				if(isnan(this->scalealpha[t][iN]))

-				{

-					FILE_LOG(logWARNING) << __PRETTY_FUNCTION__;

-					for (int jN=0; jN<this->N; jN++)

-					{

-						FILE_LOG(logWARNING) << "scalealpha["<<t-1<<"]["<<jN<<"] = " << scalealpha[t-1][jN];

-						FILE_LOG(logWARNING) << "A["<<iN<<"]["<<jN<<"] = " << A[iN][jN];

-					}

-					FILE_LOG(logWARNING) << "scalefactoralpha["<<t<<"] = "<<scalefactoralpha[t] << ", tdensities = "<<tdensities[t][iN];

-					FILE_LOG(logWARNING) << "scalealpha["<<t<<"]["<<iN<<"] = " << scalealpha[t][iN];

-					exit(1);

-				}

-			}

-		}

-

-	}

+// 	}

+// 	else if (this->use_tdens)

+// 	{

+// 

+// 		double alpha [this->N];

+// 		// Initialization

+// 		this->scalefactoralpha[0] = 0.0;

+// 		for (int iN=0; iN<this->N; iN++)

+// 		{

+// 			alpha[iN] = this->proba[iN] * this->tdensities[0][iN];

+// 			FILE_LOG(logDEBUG4) << "alpha["<<iN<<"] = " << alpha[iN];

+// 			this->scalefactoralpha[0] += alpha[iN];

+// 		}

+// 		FILE_LOG(logDEBUG4) << "scalefactoralpha["<<0<<"] = " << scalefactoralpha[0];

+// 		for (int iN=0; iN<this->N; iN++)

+// 		{

+// 			this->scalealpha[0][iN] = alpha[iN] / this->scalefactoralpha[0];

+// 			FILE_LOG(logDEBUG4) << "scalealpha["<<0<<"]["<<iN<<"] = " << scalealpha[0][iN];

+// 		}

+// 		// Induction

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

+// 		{

+// 			this->scalefactoralpha[t] = 0.0;

+// 			for (int iN=0; iN<this->N; iN++)

+// 			{

+// 				double helpsum = 0.0;

+// 				for (int jN=0; jN<this->N; jN++)

+// 				{

+// 					helpsum += this->scalealpha[t-1][jN] * this->A[jN][iN];

+// 				}

+// 				alpha[iN] = helpsum * this->tdensities[t][iN];

+// 				FILE_LOG(logDEBUG4) << "alpha["<<iN<<"] = " << alpha[iN];

+// 				this->scalefactoralpha[t] += alpha[iN];

+// 			}

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

+// 			for (int iN=0; iN<this->N; iN++)

+// 			{

+// 				this->scalealpha[t][iN] = alpha[iN] / this->scalefactoralpha[t];

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

+// 				if(isnan(this->scalealpha[t][iN]))

+// 				{

+// 					FILE_LOG(logWARNING) << __PRETTY_FUNCTION__;

+// 					for (int jN=0; jN<this->N; jN++)

+// 					{

+// 						FILE_LOG(logWARNING) << "scalealpha["<<t-1<<"]["<<jN<<"] = " << scalealpha[t-1][jN];

+// 						FILE_LOG(logWARNING) << "A["<<iN<<"]["<<jN<<"] = " << A[iN][jN];

+// 					}

+// 					FILE_LOG(logWARNING) << "scalefactoralpha["<<t<<"] = "<<scalefactoralpha[t] << ", tdensities = "<<tdensities[t][iN];

+// 					FILE_LOG(logWARNING) << "scalealpha["<<t<<"]["<<iN<<"] = " << scalealpha[t][iN];

+// 					exit(1);

+// 				}

+// 			}

+// 		}

+// 

+// 	}

 	dtime = clock() - time;

 	FILE_LOG(logDEBUG) << "forward(): " << dtime << " clicks";

@@ -489,8 +489,8 @@ void ScaleHMM::backward()

 	FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;

 	clock_t time = clock(), dtime;

-	if (not this->use_tdens)

-	{

+// 	if (not this->use_tdens)

+// 	{

 		double beta [this->N];

 		// Initialization

@@ -538,57 +538,57 @@ void ScaleHMM::backward()

 			}

 		}

-	}

-	else if (this->use_tdens)

-	{

-

-		double beta [this->N];

-		// Initialization

-		for (int iN=0; iN<this->N; iN++)

-		{

-			beta[iN] = 1.0;

-			FILE_LOG(logDEBUG4) << "beta["<<iN<<"] = " << beta[iN];

-		}

-		FILE_LOG(logDEBUG4) << "scalefactoralpha["<<T-1<<"] = " << scalefactoralpha[T-1];

-		for (int iN=0; iN<this->N; iN++)

-		{

-			this->scalebeta[T-1][iN] = beta[iN] / this->scalefactoralpha[T-1];

-			FILE_LOG(logDEBUG4) << "scalebeta["<<T-1<<"]["<<iN<<"] = " << scalebeta[T-1][iN];

-		}

-		// Induction

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

-		{

-			for (int iN=0; iN<this->N; iN++)

-			{

-				FILE_LOG(logDEBUG4) << "Calculating backward variable for state " << iN;

-				beta[iN] = 0.0;

-				for(int jN=0; jN<this->N; jN++)

-				{

-					beta[iN] += this->A[iN][jN] * this->tdensities[t+1][jN] * this->scalebeta[t+1][jN];

-				}

-			}

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

-			for (int iN=0; iN<this->N; iN++)

-			{

-				this->scalebeta[t][iN] = beta[iN] / this->scalefactoralpha[t];

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

-				if (isnan(this->scalebeta[t][iN]))

-				{

-					FILE_LOG(logWARNING) << __PRETTY_FUNCTION__;

-					for (int jN=0; jN<this->N; jN++)

-					{

-						FILE_LOG(logWARNING) << "scalebeta["<<jN<<"]["<<t+1<<"] = " << scalebeta[t+1][jN];

-						FILE_LOG(logWARNING) << "A["<<iN<<"]["<<jN<<"] = " << A[iN][jN];

-						FILE_LOG(logWARNING) << "tdensities["<<t+1<<"]["<<jN<<"] = " << tdensities[t+1][jN];

-					}

-					FILE_LOG(logWARNING) << "this->scalefactoralpha[t]["<<t<<"] = "<<this->scalefactoralpha[t] << ", tdensities = "<<densities[t][iN];

-					FILE_LOG(logWARNING) << "scalebeta["<<iN<<"]["<<t<<"] = " << scalebeta[t][iN];

-					exit(1);

-				}

-			}

-		}

-

-	}

+// 	}

+// 	else if (this->use_tdens)

+// 	{

+// 

+// 		double beta [this->N];

+// 		// Initialization

+// 		for (int iN=0; iN<this->N; iN++)

+// 		{

+// 			beta[iN] = 1.0;

+// 			FILE_LOG(logDEBUG4) << "beta["<<iN<<"] = " << beta[iN];

+// 		}

+// 		FILE_LOG(logDEBUG4) << "scalefactoralpha["<<T-1<<"] = " << scalefactoralpha[T-1];

+// 		for (int iN=0; iN<this->N; iN++)

+// 		{

+// 			this->scalebeta[T-1][iN] = beta[iN] / this->scalefactoralpha[T-1];

+// 			FILE_LOG(logDEBUG4) << "scalebeta["<<T-1<<"]["<<iN<<"] = " << scalebeta[T-1][iN];

+// 		}

+// 		// Induction

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

+// 		{

+// 			for (int iN=0; iN<this->N; iN++)

+// 			{

+// 				FILE_LOG(logDEBUG4) << "Calculating backward variable for state " << iN;

+// 				beta[iN] = 0.0;

+// 				for(int jN=0; jN<this->N; jN++)

+// 				{

+// 					beta[iN] += this->A[iN][jN] * this->tdensities[t+1][jN] * this->scalebeta[t+1][jN];

+// 				}

+// 			}

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

+// 			for (int iN=0; iN<this->N; iN++)

+// 			{

+// 				this->scalebeta[t][iN] = beta[iN] / this->scalefactoralpha[t];

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

+// 				if (isnan(this->scalebeta[t][iN]))

+// 				{

+// 					FILE_LOG(logWARNING) << __PRETTY_FUNCTION__;

+// 					for (int jN=0; jN<this->N; jN++)

+// 					{

+// 						FILE_LOG(logWARNING) << "scalebeta["<<jN<<"]["<<t+1<<"] = " << scalebeta[t+1][jN];

+// 						FILE_LOG(logWARNING) << "A["<<iN<<"]["<<jN<<"] = " << A[iN][jN];

+// 						FILE_LOG(logWARNING) << "tdensities["<<t+1<<"]["<<jN<<"] = " << tdensities[t+1][jN];

+// 					}

+// 					FILE_LOG(logWARNING) << "this->scalefactoralpha[t]["<<t<<"] = "<<this->scalefactoralpha[t] << ", tdensities = "<<densities[t][iN];

+// 					FILE_LOG(logWARNING) << "scalebeta["<<iN<<"]["<<t<<"] = " << scalebeta[t][iN];

+// 					exit(1);

+// 				}

+// 			}

+// 		}

+// 

+// 	}

 	dtime = clock() - time;

 	FILE_LOG(logDEBUG) << "backward(): " << dtime << " clicks";

@@ -641,8 +641,8 @@ void ScaleHMM::calc_sumxi()

 	}	

 	// Compute the sumxi

-	if (not this->use_tdens)

-	{

+// 	if (not this->use_tdens)

+// 	{

 		#pragma omp parallel for

 		for (int iN=0; iN<this->N; iN++)

@@ -658,25 +658,25 @@ void ScaleHMM::calc_sumxi()

 			}

 		}

-	}

-	else if (this->use_tdens)

-	{

-

-		#pragma omp parallel for

-		for (int iN=0; iN<this->N; iN++)

-		{

-			FILE_LOG(logDEBUG3) << "Calculating sumxi["<<iN<<"][jN]";

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

-			{

-				for (int jN=0; jN<this->N; jN++)

-				{

-					xi = this->scalealpha[t][iN] * this->A[iN][jN] * this->tdensities[t+1][jN] * this->scalebeta[t+1][jN];

-					this->sumxi[iN][jN] += xi;

-				}

-			}

-		}

-

-	}

+// 	}

+// 	else if (this->use_tdens)

+// 	{

+// 

+// 		#pragma omp parallel for

+// 		for (int iN=0; iN<this->N; iN++)

+// 		{

+// 			FILE_LOG(logDEBUG3) << "Calculating sumxi["<<iN<<"][jN]";

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

+// 			{

+// 				for (int jN=0; jN<this->N; jN++)

+// 				{

+// 					xi = this->scalealpha[t][iN] * this->A[iN][jN] * this->tdensities[t+1][jN] * this->scalebeta[t+1][jN];

+// 					this->sumxi[iN][jN] += xi;

+// 				}

+// 			}

+// 		}

+// 

+// 	}

 	dtime = clock() - time;

 	FILE_LOG(logDEBUG) << "calc_sumxi(): " << dtime << " clicks";

@@ -708,18 +708,18 @@ void ScaleHMM::calc_densities()

 		this->densityFunctions[iN]->calc_densities(this->densities[iN]);

 	}

-	if (this->use_tdens)

-	{

-

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

-		{

-			for (int iN=0; iN<this->N; iN++)

-			{

-				this->tdensities[t][iN] = this->densities[iN][t];

-			}

-		}

-

-	}

+// 	if (this->use_tdens)

+// 	{

+// 

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

+// 		{

+// 			for (int iN=0; iN<this->N; iN++)

+// 			{

+// 				this->tdensities[t][iN] = this->densities[iN][t];

+// 			}

+// 		}

+// 

+// 	}

 	dtime = clock() - time;

 	FILE_LOG(logDEBUG) << "calc_densities(): " << dtime << " clicks";

@@ -46,7 +46,7 @@ class ScaleHMM  {

 		double** scalealpha; ///< matrix [T x N] of forward probabilities

 		double** scalebeta; ///<  matrix [T x N] of backward probabilities

 		double** densities; ///< matrix [N x T] of density values

-		double** tdensities; ///< matrix [T x N] of density values, for use in multivariate

+// 		double** tdensities; ///< matrix [T x N] of density values, for use in multivariate !increases speed, but on cost of RAM usage and that seems to be limiting

 		double* sumgamma; ///< vector[N] of sum of posteriors (gamma values)

 		double** sumxi; ///< matrix[N x N] of xi values

 		double** gamma; ///< matrix[N x T] of posteriors

@@ -55,7 +55,7 @@ class ScaleHMM  {

 		int baumWelchTime_real; ///< elapsed time from start of the 0th iteration

 		int sumdiff_state_last; ///< sum of the difference in the state 1 assignments from one iteration to the next

 		double sumdiff_posterior; ///< sum of the difference in posterior (gamma) values from one iteration to the next

-		bool use_tdens; ///< switch for using the tdensities in the calculations

+// 		bool use_tdens; ///< switch for using the tdensities in the calculations

 		// Methods

 		void forward(); ///< calculate forward variables (alpha)