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

     } # loop over offsets

     ### Find maximum posterior for each bin between offsets

-    ptm <- startTimedMessage("Finding maximum posterior 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 {

@@ -360,13 +360,37 @@ runMultivariate <- function(bins, info, comb.states, use.states, distributions,

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

     }

     rm(aposteriors)

+    stopTimedMessage(ptm)

+    

     # Average and normalize counts to RPKM

-    counts.step <- apply(X = acounts.step, MARGIN = c(1,2), FUN = mean, na.rm=TRUE)

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

-    rm(acounts.step)

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

+    # Start stuff to call C code

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

+        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

-    ind <- apply(aposteriors.step, 1, which.max)

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

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

+    ind <- z$ind_max

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

     numstates <- length(comb.states)

     ind <- ceiling(ind / numstates)

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

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

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

             if (verbosity>=1) message("------------------------- Rerunning try ",indexmax," with eps = ",eps," -------------------------")

             hmm <- .C("C_univariate_hmm",

-                counts = as.integer(counts), # double* O

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

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

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

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

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

     } # loop over offsets

     ### Find maximum posterior for each bin between offsets

-    ptm <- startTimedMessage("Finding maximum posterior 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 {

@@ -454,19 +454,45 @@ callPeaksUnivariateAllChr <- function(binned.data, input.data=NULL, eps=0.01, in

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

     }

     rm(aposteriors)

+    stopTimedMessage(ptm)

+    

     # Average and normalize counts to RPKM

-    counts.step <- apply(X = acounts.step, MARGIN = 1, FUN = mean, na.rm=TRUE, drop=FALSE)

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

-    rm(acounts.step)

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

+    # Start stuff to call C code

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

     ## Find offset that maximizes the posteriors for each bin

-    ind <- apply(aposteriors.step, 1, which.max)

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

+    # Start stuff to call C code

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

     }

+    rm(aposteriors.step)

     stepbins$posteriors <- posteriors.step

     stopTimedMessage(ptm)

@@ -36,11 +36,13 @@ dnbinom.variance <- function(size, prob) {

 rpkm.vector <- function(counts, binsize) {

     rpkm <- counts / sum(as.numeric(counts)) * 1e6 * 1000 / binsize

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

     return(rpkm)

 }

 rpkm.matrix <- function(counts, binsize) {

     rpkm <- sweep(counts, MARGIN = 2, STATS = colSums(counts), FUN = '/')

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

     rpkm <- rpkm * 1e6 * 1000 / binsize

     return(rpkm)

 }

\ No newline at end of file

@@ -483,3 +483,69 @@ void multivariate_cleanup(int* Nmod)

 	FreeIntMatrix(multiO, *Nmod);

 }

+

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

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

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

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

+{

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

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

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

+  {

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

+    {

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

+      {

+  			value_per_i0[i1*dim[2]+i2] = array3D[(i1*dim[2]+i2) * dim[0] + i0];

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

+      }

+    }

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

+  }

+	

+}

+

+

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

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

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

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

+{

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

+  double sum=0;

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

+    }

+  }

+	

+}

\ No newline at end of file

@@ -28,3 +28,11 @@ void univariate_cleanup();

 extern "C"

 void multivariate_cleanup(int* Nmod);

+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

@@ -6,12 +6,18 @@

 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 arg4[] = {INTSXP};

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

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

+R_NativePrimitiveArgType arg7[] = {REALSXP, 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_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},

     {NULL, NULL, 0, NULL}

 };