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

     ### Define cleanup behaviour ###

     nummod <- dim(bins$counts)[2]

+    offsets <- dimnames(bins$counts)[[3]]

     on.exit(.C("C_multivariate_cleanup", as.integer(nummod)))

     # Prepare input for C function

@@ -229,9 +230,23 @@ runMultivariate <- function(bins, info, comb.states, use.states, distributions,

         startProbs.initial <- rep(1/length(comb.states), length(comb.states))

     }

-    offsets <- dimnames(bins$counts)[[3]]

-    states.list <- list()

-    aposteriors <- array(NA, dim = c(length(bins), length(comb.states), length(offsets)), dimnames = list(bin=NULL, state=comb.states, offset=offsets))

+    ### Arrays for finding maximum posterior for each bin between offsets

+    ## Make bins with offset

+    ptm <- startTimedMessage("Making bins with offsets ...")

+    if (length(offsets) > 1) {

+        stepbins <- suppressMessages( fixedWidthBins(chrom.lengths = seqlengths(bins), binsizes = as.numeric(offsets[2]))[[1]] )

+    } else {

+        stepbins <- bins

+        mcols(stepbins) <- NULL

+    }

+    sbins <- bins

+    mcols(sbins) <- NULL

+    aposteriors.step <- array(0, dim = c(length(stepbins), length(comb.states), 2), dimnames = list(bin=NULL, state=comb.states, offset=c('previousOffsets', 'currentOffset')))

+    acounts.step <- array(0, dim = c(length(stepbins), nummod, 2), dimnames = list(bin=NULL, track=info$ID, offset=c('previousOffsets', 'currentOffset')))

+    amaxPosterior.step <- array(0, dim = c(length(stepbins), 2), dimnames = list(bin=NULL, offset=c('previousOffsets', 'currentOffset')))

+    astates.step <- array(0, dim = c(length(stepbins), 2), dimnames = list(bin=NULL, offset=c('previousOffsets', 'currentOffset')))

+    

+    ### Loop over offsets ###

     for (ioffset in 1:length(offsets)) {

         offset <- offsets[ioffset]

@@ -239,8 +254,8 @@ runMultivariate <- function(bins, info, comb.states, use.states, distributions,

             messageU("Obtaining states for step size = ", offset, overline = '-', underline = NULL)

             ## Run only one iteration (no updating) if we are already over ioffset=1

             max.iter <- 1

-            transitionProbs.initial <- hmm$A

-            startProbs.initial <- hmm$proba

+            transitionProbs.initial <- hmm.A

+            startProbs.initial <- hmm.proba

             verbosity <- 0

         }

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

             densities = double(length=lenDensities), # double* densities

             keep.densities = as.logical(keep.densities), # bool* keep_densities

             states = integer(length=length(bins)), # int* states

+            maxPosterior = double(length=length(bins)), # double* maxPosterior

             A = double(length=length(comb.states)*length(comb.states)), # double* A

             proba = double(length=length(comb.states)), # double* proba

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

@@ -328,134 +344,121 @@ runMultivariate <- function(bins, info, comb.states, use.states, distributions,

                 war <- warning("HMM did not converge!\n")

             }

         }

-        ## Store states and posteriors in list

-        states.list[[offset]] <- hmm$states

-        aposteriors[,, offset] <- hmm$posteriors

+        

+        ## Store counts and posteriors in list

+        dim(hmm$posteriors) <- c(length(bins), length(comb.states))

+        dimnames(hmm$posteriors) <- list(bin=NULL, state=comb.states)

+        dim(hmm$counts) <- c(length(bins), nummod)

+        dimnames(hmm$counts) <- list(bin=NULL, track=info$ID)

+        if (keep.densities) {

+            densities <- hmm$densities

+        }

+        hmm.A <- hmm$A

+        hmm.proba <- hmm$proba

+        

+        ## Inflate posteriors to new offset

+        bins.shift <- suppressWarnings( shift(sbins, shift = as.numeric(offset)) )

+        ind <- findOverlaps(stepbins, bins.shift)

+        aposteriors.step[ind@from, , 'currentOffset'] <- hmm$posteriors[ind@to, , drop=FALSE]

+        acounts.step[ind@from, , 'currentOffset'] <- hmm$counts[ind@to, , drop=FALSE]

+        astates.step[ind@from, 'currentOffset'] <- hmm$states[ind@to]

+        amaxPosterior.step[ind@from, 'currentOffset'] <- hmm$maxPosterior[ind@to]

+        

+        ## Sum counts

+        for (i3 in 1:dim(acounts.step)[2]) {

+            acounts.step[, i3, 'previousOffsets'] <- acounts.step[, i3, 'previousOffsets'] + acounts.step[, i3, 'currentOffset']

+        }

+        

+        ## Find offset that maximizes the posteriors for each bin

+        # Start stuff to call C code

+            # Work with changing dimensions to avoid copies being made

+            dim_amaxPosterior.step <- dim(amaxPosterior.step)

+            dimnames_amaxPosterior.step <- dimnames(amaxPosterior.step)

+            dim(amaxPosterior.step) <- NULL

+            z <- .C("C_array2D_which_max",

+                    array2D = amaxPosterior.step,

+                    dim = as.integer(dim_amaxPosterior.step),

+                    ind_max = integer(dim_amaxPosterior.step[1]),

+                    value_max = double(dim_amaxPosterior.step[1]))

+            dim(amaxPosterior.step) <- dim_amaxPosterior.step

+            dimnames(amaxPosterior.step) <- dimnames_amaxPosterior.step

+            ind <- z$ind_max

+        # End stuff to call C code

+        for (i1 in 1:2) {

+            mask <- ind == i1

+            aposteriors.step[mask, , 'previousOffsets'] <- aposteriors.step[mask,,i1, drop=FALSE]

+            astates.step[mask, 'previousOffsets'] <- astates.step[mask,i1, drop=FALSE]

+            amaxPosterior.step[mask, 'previousOffsets'] <- amaxPosterior.step[mask,i1, drop=FALSE]

+        }

         message("Time spent in multivariate HMM: ", appendLF=FALSE)

         stopTimedMessage(ptm)

         ptm <- proc.time()

+        rm(hmm)

     } # loop over offsets

+    states.step <- astates.step[, 'previousOffsets']

+    rm(amaxPosterior.step, astates.step)

-    ### Find maximum posterior for each bin between offsets

-    ## Make bins with offset

-    ptm <- startTimedMessage("Making bins with offsets ...")

-    if (length(offsets) > 1) {

-        stepbins <- suppressMessages( fixedWidthBins(chrom.lengths = seqlengths(bins), binsizes = as.numeric(offsets[2]))[[1]] )

+    # Average and normalize counts to RPKM

+    ptm <- startTimedMessage("Collecting counts and posteriors over offsets ...")

+    counts.step <- acounts.step[,, 'previousOffsets'] / dim(acounts.step)[2]

+    rm(acounts.step)

+    counts.step <- rpkm.matrix(counts.step, binsize = width(bins)[1])

+    stepbins$posteriors <- aposteriors.step[,,'previousOffsets']

+    stopTimedMessage(ptm)

+    

+    ## Counts ##

+    result$bincounts <- bins

+    result$bincounts$state <- NULL

+    

+    ## Bin coordinates, posteriors and states ##

+    result$bins <- GRanges(seqnames=seqnames(stepbins), ranges=ranges(stepbins))

+    result$bins$counts.rpkm <- counts.step

+    if (!is.null(use.states$state)) {

+        state.levels <- levels(use.states$state)

     } else {

-        stepbins <- bins

-        mcols(stepbins) <- NULL

+        state.levels <- comb.states

     }

-    aposteriors.step <- array(0, dim = c(length(stepbins), length(comb.states), length(offsets)), dimnames = list(bin=NULL, state=comb.states, offset=offsets))

-    acounts.step <- array(0, dim = c(length(stepbins), nummod, length(offsets)), dimnames = list(bin=NULL, track=result$info$ID, offset=offsets))

-    astates.step <- array(0, dim = c(length(stepbins), length(offsets)), dimnames = list(bin=NULL, offset=offsets))

-    ## Inflate posteriors to new offset

-    sbins <- bins

-    mcols(sbins) <- NULL

-    for (offset in offsets) {

-        bins.shift <- suppressWarnings( shift(sbins, shift = as.numeric(offset)) )

-        ind <- findOverlaps(stepbins, bins.shift)

-        aposteriors.step[ind@from, , offset] <- aposteriors[ind@to, , offset, drop=FALSE]

-        acounts.step[ind@from, , offset] <- bins$counts[ind@to, , offset, drop=FALSE]

-        astates.step[ind@from, offset] <- states.list[[offset]][ind@to]

+    result$bins$state <- factor(states.step, levels=state.levels)

+    if (get.posteriors) {

+        ptm <- startTimedMessage("Transforming posteriors to `per sample` representation ...")

+        binstates <- dec2bin(comb.states, ndigits=nummod)

+        post.per.track <- stepbins$posteriors %*% binstates

+        colnames(post.per.track) <- result$info$ID

+        result$bins$posteriors <- post.per.track

+        result$bins$peakScores <- getPeakScores(result$bins)

+        result$bins$differential.score <- differentialScoreSum(result$bins$peakScores, result$info)

+        stopTimedMessage(ptm)

+    }

+    if (keep.densities) {

+        result$bincounts$densities <- matrix(densities, ncol=length(comb.states))

+        colnames(result$bincounts$densities) <- comb.states

     }

-    rm(aposteriors)

-    stopTimedMessage(ptm)

-    

-    # Average and normalize counts to RPKM

-    ptm <- startTimedMessage("Averaging counts between offsets ...")

-    # Start stuff to call C code

-        # Work with changing dimensions to avoid copies being made

-        dim_acounts.step <- dim(acounts.step)

-        dimnames_acounts.step <- dimnames(acounts.step)

-        dim(acounts.step) <- NULL

-        z <- .C("C_array3D_mean",

-                array3D = acounts.step,

-                dim = as.integer(dim_acounts.step),

-                mean = double(dim_acounts.step[1]*dim_acounts.step[2]))

-        # dim(acounts.step) <- dim_acounts.step

-        rm(acounts.step)

-        dim(z$mean) <- dim_acounts.step[1:2]

-        dimnames(z$mean) <- dimnames_acounts.step[1:2]

-        counts.step <- z$mean

-        counts.step <- rpkm.matrix(counts.step, binsize = width(bins)[1])

-    # End stuff to call C code

-    stopTimedMessage(ptm)

-    ## Find offset that maximizes the posteriors for each bin

-    ptm <- startTimedMessage("Finding maximum posterior between offsets ...")

-    # Start stuff to call C code

-        # Work with changing dimensions to avoid copies being made

-        dim_aposteriors.step <- dim(aposteriors.step)

-        dimnames_aposteriors.step <- dimnames(aposteriors.step)

-        dim(aposteriors.step) <- NULL

-        z <- .C("C_array3D_which_max",

-                array3D = aposteriors.step,

-                dim = as.integer(dim_aposteriors.step),

-                ind_max = integer(dim_aposteriors.step[1]))

-        dim(aposteriors.step) <- dim_aposteriors.step

-        ind <- z$ind_max

-    # End stuff to call C code

-    numstates <- length(comb.states)

-    ind <- ceiling(ind / numstates)

-    posteriors.step <- array(0, dim = c(length(stepbins), numstates), dimnames = list(bin=NULL, state=comb.states))

-    states.step <- numeric(length=length(stepbins))

-    for (i1 in 1:length(offsets)) {

-        mask <- ind == i1

-        posteriors.step[mask,] <- aposteriors.step[mask,,i1, drop=FALSE]

-        states.step[mask] <- astates.step[mask, i1, drop=FALSE]

+    ## Add combinations ##

+    if (!is.null(use.states)) {

+        result$bins$combination <- factor(mapping[as.character(result$bins$state)], levels=levels(use.states$combination))

+    } else {

+        result$bins$combination <- result$bins$state

     }

-    stepbins$posteriors <- posteriors.step

-    stopTimedMessage(ptm)

-    ## Counts

-        result$bincounts <- bins

-        result$bincounts$state <- NULL

-    ## Bin coordinates, posteriors and states

-        result$bins <- GRanges(seqnames=seqnames(stepbins), ranges=ranges(stepbins))

-        result$bins$counts.rpkm <- counts.step

-        if (!is.null(use.states$state)) {

-            state.levels <- levels(use.states$state)

-        } else {

-            state.levels <- hmm$comb.states

-        }

-        result$bins$state <- factor(states.step, levels=state.levels)

-        if (get.posteriors) {

-            ptm <- startTimedMessage("Transforming posteriors to `per sample` representation ...")

-            binstates <- dec2bin(hmm$comb.states, ndigits=hmm$num.modifications)

-            post.per.track <- stepbins$posteriors %*% binstates

-            colnames(post.per.track) <- result$info$ID

-            result$bins$posteriors <- post.per.track

-            result$bins$peakScores <- getPeakScores(result$bins)

-            result$bins$differential.score <- differentialScoreSum(result$bins$peakScores, result$info)

-            stopTimedMessage(ptm)

-        }

-        if (keep.densities) {

-            result$bins$densities <- matrix(hmm$densities, ncol=hmm$max.states)

-            colnames(result$bins$densities) <- hmm$comb.states

-        }

-    ## Add combinations

-        if (!is.null(use.states)) {

-            result$bins$combination <- factor(mapping[as.character(result$bins$state)], levels=levels(use.states$combination))

-        } else {

-            result$bins$combination <- result$bins$state

-        }

-    ## Segmentation

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

-        if (!keep.posteriors) {

-            result$bins$posteriors <- NULL

-        }

-    ## Peaks

-        result$peaks <- list()

-        for (i1 in 1:ncol(result$segments$peakScores)) {

-            mask <- result$segments$peakScores[,i1] > 0

-            peaks <- result$segments[mask]

-            mcols(peaks) <- NULL

-            peaks$peakScores <- result$segments$peakScores[mask,i1]

-            result$peaks[[i1]] <- peaks

-        }

-        names(result$peaks) <- colnames(result$segments$peakScores)

+    ## Segmentation ##

+    result$segments <- multivariateSegmentation(result$bins, column2collapseBy='state')

+    if (!keep.posteriors) {

+        result$bins$posteriors <- NULL

+    }

+        

+    ## Peaks ##

+    result$peaks <- list()

+    for (i1 in 1:ncol(result$segments$peakScores)) {

+        mask <- result$segments$peakScores[,i1] > 0

+        peaks <- result$segments[mask]

+        mcols(peaks) <- NULL

+        peaks$peakScores <- result$segments$peakScores[mask,i1]

+        result$peaks[[i1]] <- peaks

+    }

+    names(result$peaks) <- colnames(result$segments$peakScores)

     return(result)

@@ -170,6 +170,7 @@ callPeaksUnivariateAllChr <- function(binned.data, input.data=NULL, eps=0.01, in

     }

     numstates <- length(state.labels)

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

+    offsets <- dimnames(binned.data$counts)[[2]]

     ### Assign initial parameters ###

     if (class(binned.data) == class.univariate.hmm) {

@@ -199,9 +200,23 @@ callPeaksUnivariateAllChr <- function(binned.data, input.data=NULL, eps=0.01, in

     binsize <- width(binned.data)[1]

     if (keep.densities) { lenDensities <- numbins * numstates } else { lenDensities <- 1 }

-    offsets <- dimnames(binned.data$counts)[[2]]

-    counts.list <- list()

-    aposteriors <- array(NA, dim = c(numbins, numstates, length(offsets)), dimnames = list(bin=NULL, state=state.labels, offset=offsets))

+    ### Arrays for finding maximum posterior for each bin between offsets

+    ## Make bins with offset

+    ptm <- startTimedMessage("Making bins with offsets ...")

+    if (length(offsets) > 1) {

+        stepbins <- suppressMessages( fixedWidthBins(chrom.lengths = seqlengths(binned.data), binsizes = as.numeric(offsets[2]))[[1]] )

+    } else {

+        stepbins <- binned.data

+        mcols(stepbins) <- NULL

+    }

+    bins <- binned.data

+    mcols(bins) <- NULL

+    aposteriors.step <- array(0, dim = c(length(stepbins), numstates, 2), dimnames = list(bin=NULL, state=state.labels, offset=c('previousOffsets', 'currentOffset')))

+    acounts.step <- array(0, dim = c(length(stepbins), 2), dimnames = list(bin=NULL, offset=c('previousOffsets', 'currentOffset')))

+    amaxPosterior.step <- array(0, dim = c(length(stepbins), 2), dimnames = list(bin=NULL, offset=c('previousOffsets', 'currentOffset')))

+    astates.step <- array(0, dim = c(length(stepbins), 2), dimnames = list(bin=NULL, offset=c('previousOffsets', 'currentOffset')))

+    

+    ### Loop over offsets ###

     for (ioffset in 1:length(offsets)) {

         offset <- offsets[ioffset]

         counts <- binned.data$counts[,offset, drop=FALSE]

@@ -302,6 +317,8 @@ callPeaksUnivariateAllChr <- function(binned.data, input.data=NULL, eps=0.01, in

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

                 densities = double(length=lenDensities), # double* densities

                 keep.densities = as.logical(keep.densities), # bool* keep_densities

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

+                maxPosterior = double(length=numbins), # double* maxPosterior

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

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

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

@@ -357,6 +374,8 @@ callPeaksUnivariateAllChr <- function(binned.data, input.data=NULL, eps=0.01, in

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

                 densities = double(length=lenDensities), # double* densities

                 keep.densities = as.logical(keep.densities), # bool* keep_densities

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

+                maxPosterior = double(length=numbins), # double* maxPosterior

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

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

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

@@ -426,128 +445,114 @@ callPeaksUnivariateAllChr <- function(binned.data, input.data=NULL, eps=0.01, in

             ## Add class

                 class(result) <- class.univariate.hmm

         }

+        

+        ptm <- startTimedMessage("Maximizing over offsets ...")

         ## Store counts and posteriors in list

-        counts.list[[offset]] <- hmm$counts

-        aposteriors[,, offset] <- hmm$posteriors

+        dim(hmm$posteriors) <- c(numbins, numstates)

+        dimnames(hmm$posteriors) <- list(bin=NULL, state=state.labels)

+        binned.data$counts[,offset] <- hmm$counts

+        if (keep.densities) {

+            densities <- hmm$densities

+        }

+        ## Inflate posteriors to new offset

+        bins.shift <- suppressWarnings( shift(bins, shift = as.numeric(offset)) )

+        ind <- findOverlaps(stepbins, bins.shift)

+        aposteriors.step[ind@from, , 'currentOffset'] <- hmm$posteriors[ind@to, , drop=FALSE]

+        acounts.step[ind@from, 'currentOffset'] <- hmm$counts[ind@to]

+        astates.step[ind@from, 'currentOffset'] <- hmm$states[ind@to]

+        amaxPosterior.step[ind@from, 'currentOffset'] <- hmm$maxPosterior[ind@to]

+        

+        ## Sum counts

+        acounts.step[, 'previousOffsets'] <- rowSums(acounts.step)

+        ## Find offset that maximizes the posteriors for each bin

+        # Start stuff to call C code

+            # Work with changing dimensions to avoid copies being made

+            dim_amaxPosterior.step <- dim(amaxPosterior.step)

+            dimnames_amaxPosterior.step <- dimnames(amaxPosterior.step)

+            dim(amaxPosterior.step) <- NULL

+            z <- .C("C_array2D_which_max",

+                    array2D = amaxPosterior.step,

+                    dim = as.integer(dim_amaxPosterior.step),

+                    ind_max = integer(dim_amaxPosterior.step[1]),

+                    value_max = double(dim_amaxPosterior.step[1]))

+            dim(amaxPosterior.step) <- dim_amaxPosterior.step

+            dimnames(amaxPosterior.step) <- dimnames_amaxPosterior.step

+            ind <- z$ind_max

+        # End stuff to call C code

+        for (i1 in 1:2) {

+            mask <- ind == i1

+            aposteriors.step[mask, , 'previousOffsets'] <- aposteriors.step[mask,,i1, drop=FALSE]

+            astates.step[mask, 'previousOffsets'] <- astates.step[mask,i1, drop=FALSE]

+            amaxPosterior.step[mask, 'previousOffsets'] <- amaxPosterior.step[mask,i1, drop=FALSE]

+        }

+        stopTimedMessage(ptm)

+        

+        rm(hmm)

     } # loop over offsets

+    rm(amaxPosterior.step, astates.step)

-    ### Find maximum posterior for each bin between offsets

-    ## Make bins with offset

-    ptm <- startTimedMessage("Making bins with offsets ...")

-    if (length(offsets) > 1) {

-        stepbins <- suppressMessages( fixedWidthBins(chrom.lengths = seqlengths(binned.data), binsizes = as.numeric(offsets[2]))[[1]] )

+    # Average and normalize counts to RPKM

+    ptm <- startTimedMessage("Collecting counts and posteriors over offsets ...")

+    counts.step <- acounts.step[, 'previousOffsets'] / ncol(acounts.step)

+    rm(acounts.step)

+    counts.step <- rpkm.vector(counts.step, binsize = binsize)

+    stepbins$posteriors <- aposteriors.step[,,'previousOffsets']

+    rm(aposteriors.step)

+    stopTimedMessage(ptm)

+    

+    ## Get states ##

+    ptm <- startTimedMessage("Calculating states from posteriors ...")

+    posteriors <- stepbins$posteriors

+    threshold <- 1-post.cutoff

+    if (control) {

+        states <- rep(1, ncol(posteriors))

+        states[ posteriors[,2] >= posteriors[,1] ] <- 2

     } else {

-        stepbins <- binned.data

-        mcols(stepbins) <- NULL

+        states <- rep(NA, ncol(posteriors))

+        states[ posteriors[,3]<threshold & posteriors[,2]<=posteriors[,1] ] <- 1

+        states[ posteriors[,3]<threshold & posteriors[,2]>posteriors[,1] ] <- 2

+        states[ posteriors[,3]>=threshold ] <- 3

     }

-    aposteriors.step <- array(0, dim = c(length(stepbins), numstates, length(offsets)), dimnames = list(bin=NULL, state=state.labels, offset=offsets))

-    acounts.step <- array(0, dim = c(length(stepbins), length(offsets)), dimnames = list(bin=NULL, offset=offsets))

-    ## Inflate posteriors to new offset

-    bins <- binned.data

-    mcols(bins) <- NULL

-    for (offset in offsets) {

-        bins.shift <- suppressWarnings( shift(bins, shift = as.numeric(offset)) )

-        ind <- findOverlaps(stepbins, bins.shift)

-        aposteriors.step[ind@from, , offset] <- aposteriors[ind@to, , offset, drop=FALSE]

-        acounts.step[ind@from, offset] <- counts.list[[offset]][ind@to]

+    states <- state.labels[states]

+    

+    ## Counts ##

+    result$bincounts <- binned.data

+    

+    ## Bin coordinates, posteriors and states ##

+    result$bins <- GRanges(seqnames=seqnames(stepbins), ranges=ranges(stepbins))

+    result$bins$counts.rpkm <- counts.step

+    mcols(result$bins)[, 'state'] <- states

+    result$bins$posteriors <- posteriors

+    if (!control) {

+        result$bins$posterior.modified <- posteriors[,'modified']

     }

-    rm(aposteriors)

+    if (keep.densities) {

+        result$bincounts$densities <- matrix(densities, ncol=numstates)

+        colnames(result$bincounts$densities) <- state.labels

+    }

+    seqlengths(result$bins) <- seqlengths(stepbins)

     stopTimedMessage(ptm)

-    # Average and normalize counts to RPKM

-    ptm <- startTimedMessage("Averaging counts between offsets ...")

-    # Start stuff to call C code

-        # Work with changing dimensions to avoid copies being made

-        dim_acounts.step <- dim(acounts.step)

-        dimnames_acounts.step <- dimnames(acounts.step)

-        dim(acounts.step) <- NULL

-        z <- .C("C_array2D_mean",

-                array2D = acounts.step,

-                dim = as.integer(dim_acounts.step),

-                mean = double(dim_acounts.step[1]))

-        # dim(acounts.step) <- dim_acounts.step

-        rm(acounts.step)

-        counts.step <- z$mean

-        counts.step <- rpkm.vector(counts.step, binsize = binsize)

-    # End stuff to call C code

-    stopTimedMessage(ptm)

+    ## Peak score as maximum posterior in that peak ##

+    result$bins$peakScores <- getPeakScore.univariate(result$bins)

-    ## Find offset that maximizes the posteriors for each bin

-    ptm <- startTimedMessage("Finding maximum posterior between offsets ...")

-    # Start stuff to call C code

-        # Work with changing dimensions to avoid copies being made

-        dim_aposteriors.step <- dim(aposteriors.step)

-        dimnames_aposteriors.step <- dimnames(aposteriors.step)

-        dim(aposteriors.step) <- NULL

-        z <- .C("C_array3D_which_max",

-                array3D = aposteriors.step,

-                dim = as.integer(dim_aposteriors.step),

-                ind_max = integer(dim_aposteriors.step[1]))

-        dim(aposteriors.step) <- dim_aposteriors.step

-        ind <- z$ind_max

-    # End stuff to call C code

-    ind <- ceiling(ind / numstates)

-    posteriors.step <- array(0, dim = c(length(stepbins), numstates), dimnames = list(bin=NULL, state=state.labels))

-    for (i1 in 1:length(offsets)) {

-        mask <- ind == i1

-        posteriors.step[mask,] <- aposteriors.step[mask,,i1, drop=FALSE]

+    ## Segmentation ##

+    ptm <- startTimedMessage("Making segmentation ...")

+    df <- as.data.frame(result$bins)

+    red.df <- suppressMessages(collapseBins(df, column2collapseBy='state', columns2drop=c('width',grep('posterior', names(df), value=TRUE), 'counts.rpkm')))

+    result$segments <- methods::as(red.df, 'GRanges')

+    seqlengths(result$segments) <- seqlengths(binned.data)[seqlevels(result$segments)]

+    if (!keep.posteriors) {

+        result$bins$posteriors <- NULL

     }

-    rm(aposteriors.step)

-    stepbins$posteriors <- posteriors.step

     stopTimedMessage(ptm)

-    ## Get states

-        ptm <- startTimedMessage("Calculating states from posteriors ...")

-        hmm$posteriors <- stepbins$posteriors

-        threshold <- 1-post.cutoff

-        if (control) {

-            states <- rep(1, ncol(hmm$posteriors))

-            states[ hmm$posteriors[,2] >= hmm$posteriors[,1] ] <- 2

-        } else {

-            states <- rep(NA, ncol(hmm$posteriors))

-            states[ hmm$posteriors[,3]<threshold & hmm$posteriors[,2]<=hmm$posteriors[,1] ] <- 1

-            states[ hmm$posteriors[,3]<threshold & hmm$posteriors[,2]>hmm$posteriors[,1] ] <- 2

-            states[ hmm$posteriors[,3]>=threshold ] <- 3

-        }

-        states <- state.labels[states]

-    ## Counts

-        result$bincounts <- binned.data

-        for (offset in offsets) {

-            result$bincounts$counts[, offset] <- counts.list[[offset]]

-        }

-        rm(counts.list)

-        seqlengths(result$bincounts) <- seqlengths(binned.data)

-    ## Bin coordinates, posteriors and states

-        result$bins <- GRanges(seqnames=seqnames(stepbins), ranges=ranges(stepbins))

-        result$bins$counts.rpkm <- counts.step

-        mcols(result$bins)[, 'state'] <- states

-        result$bins$posteriors <- hmm$posteriors

-        if (!control) {

-            result$bins$posterior.modified <- hmm$posteriors[,'modified']

-        }

-        if (keep.densities) {

-            result$bins$densities <- matrix(hmm$densities, ncol=hmm$num.states)

-        }

-        seqlengths(result$bins) <- seqlengths(stepbins)

-        stopTimedMessage(ptm)

-    ## Peak score as maximum posterior in that peak

-        result$bins$peakScores <- getPeakScore.univariate(result$bins)

-    ## Segmentation

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

-        df <- as.data.frame(result$bins)

-        red.df <- suppressMessages(collapseBins(df, column2collapseBy='state', columns2drop=c('width',grep('posterior', names(df), value=TRUE), 'counts.rpkm')))

-        result$segments <- methods::as(red.df, 'GRanges')

-        seqlengths(result$segments) <- seqlengths(binned.data)[seqlevels(result$segments)]

-        if (!keep.posteriors) {

-            result$bins$posteriors <- NULL

-        }

-        stopTimedMessage(ptm)

-    ## Peaks

-        result$peaks <- result$segments[result$segments$state == 'modified']

-        result$peaks$state <- NULL

-        result$segments <- NULL

+    ## Peaks ##

+    result$peaks <- result$segments[result$segments$state == 'modified']

+    result$peaks$state <- NULL

+    result$segments <- NULL

     # Return results

@@ -45,4 +45,4 @@ rpkm.matrix <- function(counts, binsize) {

     rpkm[is.nan(rpkm)] <- 0

     rpkm <- rpkm * 1e6 * 1000 / binsize

     return(rpkm)

-}

\ No newline at end of file

+}

@@ -6,7 +6,7 @@ static int** multiO;

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

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

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

-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)

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

 {

 	// Define logging level

@@ -220,6 +220,21 @@ void univariate_hmm(int* O, int* T, int* N, double* size, double* prob, int* max

 		}

 	}

+	// Compute the states from posteriors

+	//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);

+		}

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

+		states[t] = ind_max + 1;

+		maxPosterior[t] = posterior_per_t[ind_max];

+	}

+		

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

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

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

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

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

 // 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, 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)

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

 {

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

@@ -432,6 +447,7 @@ void multivariate_hmm(int* O, int* T, int* N, int *Nmod, double* comb_states, do

 			}

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

+			maxPosterior[t] = posterior_per_t[ind_max];

 		}

 // 	}

 // 	else

@@ -506,46 +522,22 @@ void array3D_which_max(double* array3D, int* dim, int* ind_max)

 }

-

-// ====================================================================================

-// C version of apply(array2D, MARGIN = 1, FUN = mean)

-// ====================================================================================

-void array2D_mean(double* array2D, int* dim, double* mean)

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

+// C version of apply(array2D, 1, which.max) and apply(array2D, 1, max)

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

+void array2D_which_max(double* array2D, int* dim, int* ind_max, double* value_max)

 {

-  // array2D is actually a vector, but is intended to originate from a matrix in R

-  double sum=0;

+  // array2D is actually a vector, but is intended to originate from a 2D array in R

+	std::vector<double> value_per_i0(dim[1]);

   for (int i0=0; i0<dim[0]; i0++)

   {

-    sum = 0;

     for (int i1=0; i1<dim[1]; i1++)

     {

-      sum += array2D[i1*dim[0] + i0];

-    }

-    mean[i0] = sum / dim[1];

-  }

-	

-}

-

-

-// ====================================================================================

-// C version of apply(array3D, MARGIN = c(1,2), FUN = mean)

-// ====================================================================================

-void array3D_mean(double* array3D, int* dim, double* mean)

-{

-  // array3D is actually a vector, but is intended to originate from a 3D array in R

-  double sum=0;

-  for (int i0=0; i0<dim[0]; i0++)

-  {

-    for (int i1=0; i1<dim[1]; i1++)

-    {

-      sum = 0;

-      for (int i2=0; i2<dim[2]; i2++)

-      {

-        sum += array3D[(i2*dim[1] + i1)*dim[0] + i0];

-        // Rprintf("i0=%d, i1=%d, i2=%d, array3D[(i2*dim[1] + i1)*dim[0] + i0 = %d] = %g\n", i0, i1, i2, (i2*dim[1] + i1)*dim[0] + i0, array3D[(i2*dim[1] + i1)*dim[0] + i0]);

-      }

-      mean[i1*dim[0] + i0] = sum / dim[2];

+			value_per_i0[i1] = array2D[i1 * dim[0] + i0];

+      // Rprintf("i0=%d, i1=%d, value_per_i0[%d] = %g\n", i0, i1, i1, value_per_i0[i1]);

     }

+		ind_max[i0] = 1 + std::distance(value_per_i0.begin(), std::max_element(value_per_i0.begin(), value_per_i0.end()));

+    value_max[i0] = *std::max_element(value_per_i0.begin(), value_per_i0.end());

   }

 }

\ No newline at end of file

@@ -17,10 +17,10 @@

 // #endif

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

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

 extern "C"

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

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

 extern "C"

 void univariate_cleanup();

@@ -32,7 +32,4 @@ extern "C"

 void array3D_which_max(double* array3D, int* dim, int* ind_max);

 extern "C"

-void array2D_mean(double* array2D, int* dim, double* mean);

-

-extern "C"

-void array3D_mean(double* array3D, int* dim, double* mean);

\ No newline at end of file

+void array2D_which_max(double* array3D, int* dim, int* ind_max, double* value_max);

\ No newline at end of file

@@ -3,21 +3,19 @@

 #include "R_interface.h"

-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, 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 arg1[] = {INTSXP, INTSXP, INTSXP, REALSXP, REALSXP, INTSXP, INTSXP, REALSXP, REALSXP, REALSXP, LGLSXP, INTSXP, REALSXP, REALSXP, REALSXP, REALSXP, REALSXP, INTSXP, REALSXP, REALSXP, REALSXP, REALSXP, LGLSXP, INTSXP, 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, REALSXP, LGLSXP, INTSXP, INTSXP, INTSXP};

 R_NativePrimitiveArgType arg4[] = {INTSXP};

 R_NativePrimitiveArgType arg5[] = {REALSXP, INTSXP, INTSXP};

-R_NativePrimitiveArgType arg6[] = {REALSXP, INTSXP, REALSXP};

-R_NativePrimitiveArgType arg7[] = {REALSXP, INTSXP, REALSXP};

+R_NativePrimitiveArgType arg6[] = {REALSXP, INTSXP, INTSXP, REALSXP};

 static const R_CMethodDef CEntries[]  = {

-    {"C_univariate_hmm", (DL_FUNC) &univariate_hmm, 25, arg1},

-    {"C_multivariate_hmm", (DL_FUNC) &multivariate_hmm, 27, arg2},

+    {"C_univariate_hmm", (DL_FUNC) &univariate_hmm, 27, arg1},

+    {"C_multivariate_hmm", (DL_FUNC) &multivariate_hmm, 28, arg2},

     {"C_univariate_cleanup", (DL_FUNC) &univariate_cleanup, 0, NULL},

     {"C_multivariate_cleanup", (DL_FUNC) &multivariate_cleanup, 1, arg4},

     {"C_array3D_which_max", (DL_FUNC) &array3D_which_max, 3, arg5},

-    {"C_array2D_mean", (DL_FUNC) &array2D_mean, 3, arg6},

-    {"C_array3D_mean", (DL_FUNC) &array3D_mean, 3, arg7},

+    {"C_array2D_which_max", (DL_FUNC) &array2D_which_max, 4, arg6},

     {NULL, NULL, 0, NULL}

 };