@@ -269,6 +269,20 @@ Chromstar <- function(inputfolder, experiment.table, outputfolder, configfile=NU

     } else {

         combined.model <- loadHmmsFromFiles(savename, check.class=class.combined.multivariate.hmm)[[1]]

     }

+    ## Plot correlations

+    ptm <- startTimedMessage("Plotting read count correlation ...")

+    char.per.cm <- 10

+    legend.cm <- 3

+    savename <- file.path(plotpath, 'read-count-correlation.pdf')

+    ggplt <- plotMultivariateCorrelation(combined.model, cluster=FALSE)

+    width <- length(combined.model$info$ID) + max(sapply(combined.model$info$ID, nchar)) / char.per.cm + legend.cm

+    height <- length(combined.model$info$ID) + max(sapply(combined.model$info$ID, nchar)) / char.per.cm

+    ggsave(savename, plot=ggplt, width=width, height=height, limitsize=FALSE, units='cm')

+    savename <- file.path(plotpath, 'read-count-correlation-clustered.pdf')

+    ggplt <- plotMultivariateCorrelation(combined.model, cluster=TRUE)

+    width <- length(combined.model$info$ID) + max(sapply(combined.model$info$ID, nchar)) / char.per.cm + legend.cm

+    height <- length(combined.model$info$ID) + max(sapply(combined.model$info$ID, nchar)) / char.per.cm

+    ggsave(savename, plot=ggplt, width=width, height=height, limitsize=FALSE, units='cm')

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

     ## Export browser files ##

@@ -297,17 +311,6 @@ Chromstar <- function(inputfolder, experiment.table, outputfolder, configfile=NU

         ptm <- startTimedMessage("Making plots ...")

         char.per.cm <- 10

         legend.cm <- 3

-        savename1 <- paste0(savename, '_correlation.pdf')

-        ggplt <- suppressMessages( graphics::plot(multimodel, type='correlation', cluster=FALSE) )

-        width <- length(multimodel$info$ID) + max(sapply(multimodel$info$ID, nchar)) / char.per.cm + legend.cm

-        height <- length(multimodel$info$ID) + max(sapply(multimodel$info$ID, nchar)) / char.per.cm

-        ggsave(savename1, plot=ggplt, width=width, height=height, limitsize=FALSE, units='cm')

-

-        savename1.1 <- paste0(savename, '_correlation-clustered.pdf')

-        ggplt <- suppressMessages( graphics::plot(multimodel, type='correlation', cluster=TRUE) )

-        width <- length(multimodel$info$ID) + max(sapply(multimodel$info$ID, nchar)) / char.per.cm + legend.cm

-        height <- length(multimodel$info$ID) + max(sapply(multimodel$info$ID, nchar)) / char.per.cm

-        ggsave(savename1.1, plot=ggplt, width=width, height=height, limitsize=FALSE, units='cm')

         savename2 <- paste0(savename, '_transitionMatrix.pdf')

         ggplt <- suppressMessages( graphics::plot(multimodel, type='transitionMatrix') )

@@ -232,7 +232,7 @@ runMultivariate <- function(bins, info, comb.states, use.states, distributions,

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

         max.states = as.integer(length(comb.states)), # int* N

         num.modifications = as.integer(ncol(bins$counts)), # int* Nmod

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

+        comb.states = as.numeric(comb.states), # double* comb_states

         size = as.double(rs), # double* size

         prob = as.double(ps), # double* prob

         w = as.double(ws), # double* w

@@ -359,6 +359,9 @@ runMultivariate <- function(bins, info, comb.states, use.states, distributions,

 prepareMultivariate = function(hmms, use.states=NULL, max.states=NULL, chromosomes=NULL) {

     nummod <- length(hmms)

+    if (nummod > 53) {

+        stop("We can't handle more than 53 'hmms' Please decrease the number of input 'hmms'.")

+    }

     ## Load first HMM for coordinates

     ptm <- startTimedMessage("Getting coordinates ...")

@@ -410,7 +413,7 @@ prepareMultivariate = function(hmms, use.states=NULL, max.states=NULL, chromosom

     remove(hmm)

     ## Transform binary to decimal

-    ptm <- startTimedMessage("Getting combinatorial states")

+    ptm <- startTimedMessage("Getting combinatorial states ...")

     decimal_states <- rep(0,length(bins))

     for (imod in 1:nummod) {

         decimal_states <- decimal_states + 2^(nummod-imod) * bins$binary_statesmatrix[,imod]

@@ -431,7 +434,7 @@ prepareMultivariate = function(hmms, use.states=NULL, max.states=NULL, chromosom

     stopTimedMessage(ptm)

     ## We pre-compute the z-values for each number of counts

-    ptm <- startTimedMessage("Computing pre z-matrix...")

+    ptm <- startTimedMessage("Computing pre z-matrix ...")

     z.per.read <- array(NA, dim=c(maxcounts+1, nummod, 2))

     xcounts = 0:maxcounts

     for (imod in 1:nummod) {

@@ -463,7 +466,7 @@ prepareMultivariate = function(hmms, use.states=NULL, max.states=NULL, chromosom

     stopTimedMessage(ptm)

     ## Compute the z matrix

-    ptm <- startTimedMessage("Transfering values into z-matrix...")

+    ptm <- startTimedMessage("Transfering values into z-matrix ...")

     z.per.bin = array(NA, dim=c(length(bins), nummod, 2), dimnames=list(bin=1:length(bins), track=info$ID, state.labels[2:3]))

     for (imod in 1:nummod) {

         for (i1 in 1:2) {

@@ -476,18 +479,30 @@ prepareMultivariate = function(hmms, use.states=NULL, max.states=NULL, chromosom

     stopTimedMessage(ptm)

     ### Calculate correlation matrix

-    ptm <- startTimedMessage("Computing inverse of correlation matrix...")

+    ptm <- startTimedMessage("Computing inverse of correlation matrix ...")

     correlationMatrix = array(NA, dim=c(nummod,nummod,length(comb.states)), dimnames=list(track=info$ID, track=info$ID, comb.state=comb.states))

     correlationMatrixInverse = array(NA, dim=c(nummod,nummod,length(comb.states)), dimnames=list(track=info$ID, track=info$ID, comb.state=comb.states))

     determinant = rep(NA, length(comb.states))

     names(determinant) <- comb.states

+    numericToBin <- function(x, num.digits) {

+        bin <- logical(length=num.digits)

+        xi <- x

+        i1 <- num.digits

+        while (xi > 0) {

+            bin[i1] <- xi %% 2

+            xi <- xi %/% 2

+            i1 <- i1 - 1

+        }

+        return(bin)

+    }

+    

     ## Calculate correlation matrix serial

     for (state in comb.states) {

         istate = which(comb.states==state)

-        mask = which(bins$state==as.integer(state))

+        mask = which(bins$state==as.numeric(state))

         # Convert state to binary representation

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

+        binary_state <- numericToBin(as.numeric(state), nummod)

         # Subselect z

         z.temp <- matrix(NA, ncol=nummod, nrow=length(mask))

         for (i1 in 1:length(binary_state)) {

@@ -79,11 +79,13 @@ combineMultivariates <- function(hmms, mode) {

         hmm <- suppressMessages( loadHmmsFromFiles(hmms[[1]], check.class=class.multivariate.hmm)[[1]] )

         bins <- hmm$bins

         mcols(bins) <- NULL

-        ## Add combinatorial states and posteriors

+        ## Add combinatorial states, counts and posteriors

         infos <- list()

         infos[[1]] <- hmm$info

         combs <- list()

         combs[[1]] <- hmm$bins$combination

+        counts <- list()

+        counts[[1]] <- hmm$bins$counts

         posteriors <- list()

         posteriors[[1]] <- hmm$bins$posteriors

         binstates <- list()

@@ -96,11 +98,13 @@ combineMultivariates <- function(hmms, mode) {

                 hmm <- suppressMessages( loadHmmsFromFiles(hmms[[i1]], check.class=class.multivariate.hmm)[[1]] )

                 infos[[i1]] <- hmm$info

                 combs[[i1]] <- hmm$bins$combination

+                counts[[i1]] <- hmm$bins$counts

                 posteriors[[i1]] <- hmm$bins$posteriors

                 binstates[[i1]] <- dec2bin(hmm$bins$state, colnames=hmm$info$ID)

                 stopTimedMessage(ptm)

             }

         }

+        counts <- do.call(cbind, counts)

         posteriors <- do.call(cbind, posteriors)

         infos <- do.call(rbind, infos)

         conditions <- unique(infos$condition)

@@ -112,16 +116,19 @@ combineMultivariates <- function(hmms, mode) {

     } else if (mode == 'condition') {

         ### Get posteriors and binary states

         infos <- list()

+        counts <- list()

         posteriors <- list()

         binstates <- list()

         for (i1 in 1:length(hmms)) {

             ptm <- startTimedMessage("Processing HMM ",i1," ...")

             hmm <- suppressMessages( loadHmmsFromFiles(hmms[[i1]], check.class=class.multivariate.hmm)[[1]] )

             infos[[i1]] <- hmm$info

+            counts[[i1]] <- hmm$bins$counts

             posteriors[[i1]] <- hmm$bins$posteriors

             binstates[[i1]] <- dec2bin(hmm$bins$state, colnames=hmm$info$ID)

             stopTimedMessage(ptm)

         }

+        counts <- do.call(cbind, counts)

         posteriors <- do.call(cbind, posteriors)

         infos <- do.call(rbind, infos)

         conditions <- unique(infos$condition)

@@ -157,6 +164,7 @@ combineMultivariates <- function(hmms, mode) {

         mcols(bins) <- NULL

         infos <- hmm$info

         conditions <- unique(infos$condition)

+        counts <- hmm$bins$counts

         posteriors <- hmm$bins$posteriors

         states <- hmm$bins$state

         combs <- list()

@@ -177,9 +185,11 @@ combineMultivariates <- function(hmms, mode) {

         hmm <- suppressMessages( loadHmmsFromFiles(hmms[[1]], check.class=class.multivariate.hmm)[[1]] )

         bins <- hmm$bins

         mcols(bins) <- NULL

-        ## Add combinatorial states and posteriors

+        ## Add combinatorial states, counts and posteriors

         infos <- list()

         infos[[1]] <- hmm$info

+        counts <- list()

+        counts[[1]] <- hmm$bins$counts

         posteriors <- list()

         posteriors[[1]] <- hmm$bins$posteriors

         binstates <- list()

@@ -191,11 +201,13 @@ combineMultivariates <- function(hmms, mode) {

                 ptm <- startTimedMessage("Processing mark-condition ",i1," ...")

                 hmm <- suppressMessages( loadHmmsFromFiles(hmms[[i1]], check.class=class.multivariate.hmm)[[1]] )

                 infos[[i1]] <- hmm$info

+                counts[[i1]] <- hmm$bins$counts

                 posteriors[[i1]] <- hmm$bins$posteriors

                 binstates[[i1]] <- dec2bin(hmm$bins$state, colnames=hmm$info$ID)

                 stopTimedMessage(ptm)

             }

         }

+        counts <- do.call(cbind, counts)

         posteriors <- do.call(cbind, posteriors)

         binstates <- do.call(cbind, binstates)

         infos <- do.call(rbind, infos)

@@ -230,7 +242,8 @@ combineMultivariates <- function(hmms, mode) {

     bins$state <- states

     mcols(bins)[names(combs.df)] <- combs.df

-    ## Transferring posteriors

+    ## Transferring counts and posteriors

+    bins$counts <- counts

     bins$posteriors <- posteriors

     ## Add differential score ##

@@ -239,7 +252,7 @@ combineMultivariates <- function(hmms, mode) {

     ### Redo the segmentation for all conditions combined

     ptm <- startTimedMessage("Redoing segmentation for all conditions combined ...")

     segments <- suppressMessages( multivariateSegmentation(bins, column2collapseBy='state') )

-    names(mcols(segments)) <- setdiff(names(mcols(bins)), 'posteriors')

+    names(mcols(segments)) <- setdiff(names(mcols(bins)), c('posteriors','counts'))

     stopTimedMessage(ptm)

     ### Redo the segmentation for each condition separately

@@ -51,7 +51,7 @@ unis2pseudomulti <- function(uni.hmm.list) {

     weights = lapply(uni.hmm.list,"[[","weights")

     # Extract the reads

-    ptm <- startTimedMessage("Extracting read counts from uni.hmm.list...")

+    ptm <- startTimedMessage("Extracting read counts from uni.hmm.list ...")

     reads = matrix(NA, ncol=nummod, nrow=numbins)

     colnames(reads) <- info$ID

     for (imod in 1:nummod) {

@@ -62,7 +62,7 @@ unis2pseudomulti <- function(uni.hmm.list) {

     stopTimedMessage(ptm)

     ## Get combinatorial states

-    ptm <- startTimedMessage("Getting combinatorial states")

+    ptm <- startTimedMessage("Getting combinatorial states ...")

     combstates.per.bin = combinatorialStates(uni.hmm.list)

     comb.states.table = table(combstates.per.bin)

     comb.states = as.numeric(names(sort(comb.states.table, decreasing=TRUE)))

@@ -84,7 +84,7 @@ unis2pseudomulti <- function(uni.hmm.list) {

     stopTimedMessage(ptm)

     ## Calculate transition matrix

-    ptm <- startTimedMessage("Estimating transition matrix...")

+    ptm <- startTimedMessage("Estimating transition matrix ...")

     A.estimated = matrix(0, ncol=2^nummod, nrow=2^nummod)

     colnames(A.estimated) = 1:2^nummod-1

     rownames(A.estimated) = 1:2^nummod-1

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

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

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

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

-void multivariate_hmm(int* O, int* T, int* N, int *Nmod, int* comb_states, double* size, double* prob, double* w, double* cor_matrix_inv, double* det, int* maxiter, int* maxtime, double* eps, double* posteriors, bool* keep_posteriors, double* densities, bool* keep_densities, int* states, double* A, double* proba, double* loglik, double* initial_A, double* initial_proba, bool* use_initial_params, int* num_threads, int* error, int* verbosity)

+void multivariate_hmm(int* O, int* T, int* N, int *Nmod, double* comb_states, double* size, double* prob, double* w, double* cor_matrix_inv, double* det, int* maxiter, int* maxtime, double* eps, double* posteriors, bool* keep_posteriors, double* densities, bool* keep_densities, int* states, double* A, double* proba, double* loglik, double* initial_A, double* initial_proba, bool* use_initial_params, int* num_threads, int* error, int* verbosity)

 {

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

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

 	// Prepare the binary_states (univariate) vector: binary_states[N][Nmod], e.g., binary_states[iN][imod] tells me at state comb_states[iN], modification imod is non-enriched (0) or enriched (1)

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

+	double res;

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

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

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

 	{

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

+		res = comb_states[iN];

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

 		{

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

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

-				binary_states[iN][imod] = 1;

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

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

 		}

 	}

@@ -20,7 +20,7 @@ extern "C"

 void univariate_hmm(int* O, int* T, int* N, double* size, double* prob, int* maxiter, int* maxtime, double* eps, double* posteriors, double* densities, bool* keep_densities, double* A, double* proba, double* loglik, double* weights, int* iniproc, double* initial_size, double* initial_prob, double* initial_A, double* initial_proba, bool* use_initial_params, int* num_threads, int* error, int* read_cutoff, int* verbosity);

 extern "C"

-void multivariate_hmm(int* O, int* T, int* N, int *Nmod, int* comb_states, double* size, double* prob, double* w, double* cor_matrix_inv, double* det, int* maxiter, int* maxtime, double* eps, double* posteriors, bool* keep_posteriors, double* densities, bool* keep_densities, int* states, double* A, double* proba, double* loglik, double* initial_A, double* initial_proba, bool* use_initial_params, int* num_threads, int* error, int* verbosity);

+void multivariate_hmm(int* O, int* T, int* N, int *Nmod, double* comb_states, double* size, double* prob, double* w, double* cor_matrix_inv, double* det, int* maxiter, int* maxtime, double* eps, double* posteriors, bool* keep_posteriors, double* densities, bool* keep_densities, int* states, double* A, double* proba, double* loglik, double* initial_A, double* initial_proba, bool* use_initial_params, int* num_threads, int* error, int* verbosity);

 extern "C"

 void univariate_cleanup();

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

 R_NativePrimitiveArgType arg1[] = {INTSXP, INTSXP, INTSXP, REALSXP, REALSXP, INTSXP, INTSXP, REALSXP, REALSXP, REALSXP, LGLSXP, REALSXP, REALSXP, REALSXP, REALSXP, INTSXP, REALSXP, REALSXP, REALSXP, REALSXP, LGLSXP, INTSXP, INTSXP, INTSXP, INTSXP};

-R_NativePrimitiveArgType arg2[] = {INTSXP, INTSXP, INTSXP, INTSXP, INTSXP, REALSXP, REALSXP, REALSXP, REALSXP, REALSXP, INTSXP, INTSXP, REALSXP, REALSXP, LGLSXP, REALSXP, LGLSXP, INTSXP, REALSXP, REALSXP, REALSXP, REALSXP, REALSXP, LGLSXP, INTSXP, INTSXP, INTSXP};

+R_NativePrimitiveArgType arg2[] = {INTSXP, INTSXP, INTSXP, INTSXP, REALSXP, REALSXP, REALSXP, REALSXP, REALSXP, REALSXP, INTSXP, INTSXP, REALSXP, REALSXP, LGLSXP, REALSXP, LGLSXP, INTSXP, REALSXP, REALSXP, REALSXP, REALSXP, REALSXP, LGLSXP, INTSXP, INTSXP, INTSXP};

 R_NativePrimitiveArgType arg4[] = {INTSXP};

 static const R_CMethodDef CEntries[]  = {