@@ -7,12 +7,12 @@ export(binaryA)

 export(calcCoGAPSStat)

 export(calcGeneGSStat)

 export(calcZ)

+export(computeGeneGSProb)

 export(createGWCoGAPSSets)

 export(displayBuildReport)

 export(gapsMapRun)

 export(gapsRun)

 export(generateSeeds)

-export(getParam)

 export(patternMarkers)

 export(patternMatch4Parallel)

 export(patternMatcher)

@@ -27,8 +27,6 @@ export(reOrderBySet)

 export(reconstructGene)

 export(reorderByPatternMatch)

 export(residuals)

-export(setParam)

-exportClasses(GapsParams)

 import(doParallel)

 import(foreach)

 import(ggplot2)

@@ -56,12 +54,8 @@ importFrom(graphics,points)

 importFrom(graphics,screen)

 importFrom(graphics,split.screen)

 importFrom(graphics,title)

-importFrom(methods,"slot<-")

-importFrom(methods,callNextMethod)

 importFrom(methods,is)

 importFrom(methods,new)

-importFrom(methods,slot)

-importFrom(methods,validObject)

 importFrom(reshape2,melt)

 importFrom(stats,D)

 importFrom(stats,as.dist)

@@ -8,41 +8,41 @@

 #' @return list with A and P matrix estimates

 #' @importFrom methods new

 #' @export

-CoGAPS <- function(D, S, ...)

+CoGAPS <- function(D, S, nFactor=7, nEquil=1000, nSample=1000, nOutputs=1000,

+nSnapshots=0, alphaA=0.01, alphaP=0.01, maxGibbmassA=100, maxGibbmassP=100,

+seed=-1, messages=TRUE, singleCellRNASeq=FALSE, whichMatrixFixed = 'N',

+fixedPatterns = matrix(0), checkpointInterval=0, ...)

 {

-    # process v2 style arguments for backwards compatibility

-    extraArgs <- list(...)

-    params <- oldParams(new('GapsParams', 7, 1000, 1000), extraArgs)

+    # get v2 arguments

+    oldArgs <- list(...)

+    if (!is.null(oldArgs$nOutR))

+        nOutputs <- oldArgs$nOutR

+    if (!is.null(oldArgs$max_gibbmass_paraA))

+        maxGibbmassA <- oldArgs$max_gibbmass_paraA

+    if (!is.null(oldArgs$max_gibbmass_paraP))

+        maxGibbmassP <- oldArgs$max_gibbmass_paraP

+    if (!is.null(oldArgs$sampleSnapshots) & is.null(oldArgs$numSnapshots))

+        nSnapshots <- 100

+    if (!is.null(oldArgs$sampleSnapshots) & !is.null(oldArgs$numSnapshots))

+        nSnapshots <- oldArgs$numSnapshots

+    if (missing(D) & !is.null(oldArgs$data))

+        D <- oldArgs$data

+    if (missing(S) & !is.null(oldArgs$unc))

+        S <- oldArgs$unc

-    # check data

-    if (missing(D) & length(extraArgs$data))

-        D <- extraArgs$data

-    if (missing(S) & length(extraArgs$unc))

-        S <- extraArgs$unc

+    # check arguments

+    if (class(D) != "matrix" | class(S) != "matrix")

+        stop('D and S must be matrices')

     if (any(D < 0) | any(S < 0))

         stop('D and S matrix must be non-negative')

-    checkParamsWithData(params, nrow(D), ncol(D), nrow(S), ncol(S))

     # run algorithm with call to C++ code

-    result <- cogaps_cpp(as.matrix(D), as.matrix(S),

-        as.matrix(params$fixedPatterns), params)

+    result <- cogaps_cpp(D, S, nFactor, nEquil, nEquil/10, nSample, nOutputs, nSnapshots,

+        alphaA, alphaP, maxGibbmassA, maxGibbmassP, seed, messages,

+        singleCellRNASeq, whichMatrixFixed, fixedPatterns, checkpointInterval)

     # backwards compatible with v2

-    if (!is.null(extraArgs$GStoGenes))

-    {

-        #warning('the GStoGenes argument is deprecated with v3.0')

-        if (is.null(extraArgs$plot) | extraArgs$plot)

-            plotGAPS(result$Amean, result$Pmean)

-        if (is.null(extraArgs$nPerm))

-            extraArgs$nPerm <- 500

-        GSP <- calcCoGAPSStat(result$Amean, result$Asd, extraArgs$GStoGenes,

-            extraArgs$nPerm)

-        result <- list(meanChi2=result$meanChi2, D=D, Sigma=S,

-            Amean=result$Amean, Asd=result$Asd, Pmean=result$Pmean,

-            Psd=result$Psd, GSUpreg=GSP$GSUpreg, GSDownreg=GSP$GSDownreg,

-            GSActEst=GSP$GSActEst)

-    }

-    return(result)

+    return(v2CoGAPS(result, ...))

 }

 #' Restart CoGAPS from Checkpoint File

@@ -59,11 +59,6 @@ CoGapsFromCheckpoint <- function(D, S, path)

     cogapsFromCheckpoint_cpp(D, S, path)

 }

-GWCoGapsFromCheckpoint <- function(fname)

-{

-    #TODO

-}

-

 #' Display Information About Package Compilation

 #' @export

 displayBuildReport <- function()

@@ -75,15 +70,20 @@ displayBuildReport <- function()

 #'

 #' @param D data matrix

 #' @param S uncertainty matrix

-#' @param ... v2 style parameters

 #' @return list with A and P matrix estimates

 #' @importFrom methods new

 #' @export

-gapsRun <- function(D, S, ...)

+gapsRun <- function(D, S, ABins=data.frame(), PBins=data.frame(), nFactor=7,

+simulation_id="simulation", nEquil=1000, nSample=1000, nOutR=1000,

+output_atomic=FALSE, fixedBinProbs=FALSE, fixedDomain="N", sampleSnapshots=TRUE,

+numSnapshots=100, alphaA=0.01, nMaxA=100000, max_gibbmass_paraA=100.0,

+alphaP=0.01, nMaxP=100000, max_gibbmass_paraP=100.0, seed=-1, messages=TRUE)

 {

     #warning('gapsRun is deprecated with v3.0, use CoGAPS')

-    params <- oldParams(new('GapsParams', 7, 1000, 1000), list(...))

-    CoGAPS(D, S, params)

+    CoGAPS(D, S, nFactor=nFactor, nEquil=nEquil, nSample=nSample, nOutputs=nOutR,

+        nSnapshots=ifelse(sampleSnapshots,numSnapshots,0), alphaA=alphaA,

+        alphaP=alphaP, maxGibbmassA=max_gibbmass_paraA, messages=messages,

+        maxGibbmassP=max_gibbmass_paraP, seed=seed)

 }

 #' Backwards Compatibility with v2

@@ -95,34 +95,26 @@ gapsRun <- function(D, S, ...)

 #' @return list with A and P matrix estimates

 #' @importFrom methods new

 #' @export

-gapsMapRun <- function(D, S, FP, ...)

+gapsMapRun <- function(D, S, FP, ABins=data.frame(), PBins=data.frame(),

+nFactor=5, simulation_id="simulation", nEquil=1000, nSample=1000, nOutR=1000,

+output_atomic=FALSE, fixedMatrix="P", fixedBinProbs=FALSE, fixedDomain="N",

+sampleSnapshots=TRUE, numSnapshots=100, alphaA=0.01, nMaxA=100000,

+max_gibbmass_paraA=100.0, alphaP=0.01, nMaxP=100000, max_gibbmass_paraP=100.0,

+seed=-1, messages=TRUE)

 {

     #warning('gapsMapRun is deprecated with v3.0, use CoGaps')

-    params <- oldParams(new('GapsParams', 7, 1000, 1000), list(...))

-    CoGAPS(D, S, params)

+    CoGAPS(D, S, nFactor=nFactor, nEquil=nEquil, nSample=nSample, nOutputs=nOutR,

+        nSnapshots=ifelse(sampleSnapshots,numSnapshots,0), alphaA=alphaA,

+        alphaP=alphaP, maxGibbmassA=max_gibbmass_paraA, messages=messages,

+        maxGibbmassP=max_gibbmass_paraP, seed=seed, whichMatrixFixed='P',

+        fixedPatterns=as.matrix(FP))

 }

-# helper function to process v2 parameters

-oldParams <- function(params, args)

+v2CoGAPS <- function(result, ...)

 {

-    # standard params

-    if (length(args$nFactor))      params@nFactor    <- as.integer(args$nFactor)

-    if (length(args$nEquil))       params@nEquil     <- as.integer(args$nEquil)

-    if (length(args$nSample))      params@nSample    <- as.integer(args$nSample)

-    if (length(args$alphaA))       params@alphaA     <- args$alphaA

-    if (length(args$alphaP))       params@alphaP     <- args$alphaP

-    if (length(args$seed))         params@seed       <- as.integer(args$seed)

-    if (length(args$messages))     params@messages   <- args$messages

-    if (length(args$numSnapshots))

-        params@nSnapshots <- as.integer(args$numSnapshots)

-

-    # gapsMap params

-    if (length(args$FP))

+    if (!is.null(list(...)$GStoGenes))

     {

-        params@fixedPatterns <- as.matrix(args$FP)

-        params@whichMatrixFixed <- 'P'

-    }

-    # return v3 style GapsParams

-    return(params)    

+    }

+    return(result)

 }

\ No newline at end of file

@@ -93,4 +93,9 @@ minNS=NA, ...)

     class(AP.fixed) <- append(class(AP.fixed),"CoGAPS")

     save(AP.fixed, file=paste(fname,".Rda",sep=""))

     message(paste(fname,".Rda",sep=""))

+}

+

+GWCoGapsFromCheckpoint <- function(fname)

+{

+    #TODO

 }

\ No newline at end of file

@@ -1,8 +1,8 @@

 # Generated by using Rcpp::compileAttributes() -> do not edit by hand

 # Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393

-cogaps_cpp <- function(DMatrix, SMatrix, fixedPatterns, params) {

-    .Call('_CoGAPS_cogaps_cpp', PACKAGE = 'CoGAPS', DMatrix, SMatrix, fixedPatterns, params)

+cogaps_cpp <- function(D, S, nFactor, nEquil, nEquilCool, nSample, nOutputs, nSnapshots, alphaA, alphaP, maxGibbmassA, maxGibbmassP, seed, messages, singleCellRNASeq, whichMatrixFixed, FP, checkpointInterval) {

+    .Call('_CoGAPS_cogaps_cpp', PACKAGE = 'CoGAPS', D, S, nFactor, nEquil, nEquilCool, nSample, nOutputs, nSnapshots, alphaA, alphaP, maxGibbmassA, maxGibbmassP, seed, messages, singleCellRNASeq, whichMatrixFixed, FP, checkpointInterval)

 }

 cogapsFromCheckpoint_cpp <- function(D, S, fileName) {

@@ -48,12 +48,20 @@ nullGenes=FALSE)

 #' Compute Gene Probability

 #'

+#' @details Computes the p-value for gene set membership using the CoGAPS-based

+#'  statistics developed in Fertig et al. (2012).  This statistic refines set

+#'  membership for each candidate gene in a set specified in \code{GSGenes} by

+#'  comparing the inferred activity of that gene to the average activity of the

+#'  set.

 #' @param Amean A matrix mean values

 #' @param Asd A matrix standard deviations

 #' @param GSGenes data.frame or list with gene sets

 #' @param Pw weight on genes

 #' @param numPerm number of permutations for null

 #' @param PwNull - logical indicating gene adjustment

+#' @return A vector of length GSGenes containing the p-values of set membership

+#'  for each gene containined in the set specified in GSGenes.

+#' @export

 computeGeneGSProb <- function(Amean, Asd, GSGenes, Pw=rep(1,ncol(Amean)),

 numPerm=500, PwNull=FALSE)

 {

deleted file mode 100644

@@ -1,120 +0,0 @@

-#' @title GapsParams

-#' @description Parameters for running CoGAPS

-#'

-#' @slot nFactor number of patterns (basis vectors, metagenes)

-#' @slot nEquil number of iterations for burn-in

-#' @slot nSample number of iterations for sampling

-#' @slot nOutput how often (number of iterations) to print status in R console

-#' @slot numSnapshots the number of individual samples to capture

-#' @slot alphaA sparsity parameter for A domain

-#' @slot alphaP sparsity parameter for P domain

-#' @slot maxGibbmassA limit truncated normal to max size

-#' @slot maxGibbmassP limit truncated normal to max size

-#' @slot seed positive values are kept, negative values will be overwritten

-#'  by a seed generated from the current time

-#' @slot messages display messages during the run

-#' @slot singleCelLRNASeq indicate if the data is single cell RNA-seq data

-#' @slot fixedPatterns matrix of fixed values in either A or P matrix

-#' @slot whichMatrixFixed characted indicating whether A or P contains

-#'  the fixed patterns

-#' @slot checkpointInterval how often (number of seconds) to create a checkpoint file

-#' @slot nCores max number of cores to run on

-#' @export

-setClass('GapsParams', slots = c(

-    nFactor = 'integer',

-    nEquil = 'integer',

-    nEquilCool = 'integer',

-    nSample = 'integer',

-    nOutput = 'integer',

-    nSnapshots = 'integer',

-    alphaA = 'numeric',

-    alphaP = 'numeric',

-    maxGibbmassA = 'numeric',

-    maxGibbmassP = 'numeric',

-    seed = 'integer',

-    messages  = 'logical',

-    singleCellRNASeq = 'logical',

-    fixedPatterns = 'matrix',

-    whichMatrixFixed = 'character',

-    checkpointInterval = 'integer',

-    nCores = 'integer'

-))

-

-#' @importFrom methods callNextMethod

-setMethod('initialize', 'GapsParams', 

-    function(.Object, nFactor, nEquil, nSample, ...)

-    {

-        .Object@nFactor <- as.integer(nFactor)

-        .Object@nEquil <- as.integer(nEquil)

-        .Object@nEquil <- as.integer(floor(nEquil/10))

-        .Object@nSample <- as.integer(nSample)

-        .Object@nOutput <- as.integer(1000)

-        .Object@nSnapshots <- as.integer(0)

-        .Object@alphaA <- 0.01

-        .Object@alphaP <- 0.01

-        .Object@maxGibbmassA <- 100.0

-        .Object@maxGibbmassP <- 100.0

-        .Object@seed <- as.integer(-1)

-        .Object@messages <- TRUE

-        .Object@singleCellRNASeq <- FALSE

-        .Object@fixedPatterns <- matrix(0)

-        .Object@whichMatrixFixed <- 'N'

-        .Object@checkpointInterval <- as.integer(0)

-        .Object@nCores <- as.integer(1)

-

-        .Object <- callNextMethod(.Object, ...)

-        .Object

-    }

-)

-

-setValidity('GapsParams',

-    function(object)

-    {

-        if (object@whichMatrixFixed == 'P' & nrow(object@fixedPatterns) > object@nFactor)

-            'number of fixed patterns greater than nFactor'

-        if (object@whichMatrixFixed == 'A' & ncol(object@fixedPatterns) > object@nFactor)

-            'number of fixed patterns greater than nFactor'

-    }

-)

-

-##################### Generics ###################

-

-#' @export

-setGeneric('getParam', function(object, name)

-    {standardGeneric('getParam')})

-

-#' @export

-setGeneric('setParam', function(object, name, value)

-    {standardGeneric('setParam')})

-

-##################### Methods ####################

-

-#' @importFrom methods slot

-setMethod("getParam", "GapsParams", function(object, name)

-{

-    slot(object, name)

-})

-

-#' @importFrom methods slot<- validObject

-setMethod("setParam", "GapsParams", function(object, name, value)

-{

-    slot(object, name) <- value

-    validObject(object)

-    return(object)

-})

-

-setMethod("show", "GapsParams", function(object)

-{

-    cat("CoGAPS Parameters Object\n")

-    cat("Number of patterns to estimate: ", object@nFactor)

-})

-

-checkParamsWithData <- function(params, nRD, nCD, nRS, nCS)

-{

-    if (nRD != nRS | nCD != nCS)

-        stop('D and S matrices have mismached dimensions')

-    if (params@whichMatrixFixed == 'A' & nrow(params@fixedPatterns) != nRD)

-        stop('invalid fixed pattern length')

-    if (params@whichMatrixFixed == 'P' & ncol(params@fixedPatterns) != nCD)

-        stop('invalid fixed pattern length')

-}

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

 \alias{CoGAPS}

 \title{CoGAPS Matrix Factorization Algorithm}

 \usage{

-CoGAPS(D, S, ...)

+CoGAPS(D, S, nFactor = 7, nEquil = 1000, nSample = 1000,

+  nOutputs = 1000, nSnapshots = 0, alphaA = 0.01, alphaP = 0.01,

+  maxGibbmassA = 100, maxGibbmassP = 100, seed = -1, messages = TRUE,

+  singleCellRNASeq = FALSE, whichMatrixFixed = "N",

+  fixedPatterns = matrix(0), checkpointInterval = 0, ...)

 }

 \arguments{

 \item{D}{data matrix}

deleted file mode 100644

@@ -1,47 +0,0 @@

-% Generated by roxygen2: do not edit by hand

-% Please edit documentation in R/class-GapsParams.R

-\docType{class}

-\name{GapsParams-class}

-\alias{GapsParams-class}

-\title{GapsParams}

-\description{

-Parameters for running CoGAPS

-}

-\section{Slots}{

-

-\describe{

-\item{\code{nFactor}}{number of patterns (basis vectors, metagenes)}

-

-\item{\code{nEquil}}{number of iterations for burn-in}

-

-\item{\code{nSample}}{number of iterations for sampling}

-

-\item{\code{nOutput}}{how often (number of iterations) to print status in R console}

-

-\item{\code{numSnapshots}}{the number of individual samples to capture}

-

-\item{\code{alphaA}}{sparsity parameter for A domain}

-

-\item{\code{alphaP}}{sparsity parameter for P domain}

-

-\item{\code{maxGibbmassA}}{limit truncated normal to max size}

-

-\item{\code{maxGibbmassP}}{limit truncated normal to max size}

-

-\item{\code{seed}}{positive values are kept, negative values will be overwritten

-by a seed generated from the current time}

-

-\item{\code{messages}}{display messages during the run}

-

-\item{\code{singleCelLRNASeq}}{indicate if the data is single cell RNA-seq data}

-

-\item{\code{fixedPatterns}}{matrix of fixed values in either A or P matrix}

-

-\item{\code{whichMatrixFixed}}{characted indicating whether A or P contains

-the fixed patterns}

-

-\item{\code{checkpointInterval}}{how often (number of seconds) to create a checkpoint file}

-

-\item{\code{nCores}}{max number of cores to run on}

-}}

-

@@ -1,67 +1,37 @@

+% Generated by roxygen2: do not edit by hand

+% Please edit documentation in R/calcGeneGSStat.R

 \name{computeGeneGSProb}

 \alias{computeGeneGSProb}

-\alias{geneGSProb}

-\title{CoGAPS gene membership statistic}

-

-\description{Computes the p-value for gene set membership using the CoGAPS-based statistics developed in Fertig et al. (2012).  This statistic refines set membership for each candidate gene in a set specified in \code{GSGenes} by comparing the inferred activity of that gene to the average activity of the set.  Specifically, we compute the following summary statistic for each gene \eqn{g} that is a candidate member of gene set \eqn{G}:

-\deqn{S_{g,G} = (\sum_{p} -log(Pr_{G,p})Pw[p](A_{gp}/\sigma_{gp})) / \sum_{p}-log(Pr_{G,p})Pw[p],}

-where \eqn{p} indexes each of the patterns, \eqn{Pr_{G,p}} is the probability that gene set \eqn{G} is upregulated computed with \code{\link{calcCoGAPSStat}}, \eqn{A_{gp}} is the mean amplitude matrix from the \code{GAPS} matrix factorization, Pw[p] is a prior weighting for each pattern based upon the context to which that pattern relates, and \eqn{\sigma_{gp}} is the standard deviation of the amplitude matrix. P-values are formulated from a permutation test comparing the value of \eqn{S_{g,G}} for genes in \code{GSGenes} relative to the value of \eqn{S_{g,G}} \code{numPerm} random gene sets with the same number of targets.

+\title{Compute Gene Probability}

+\usage{

+computeGeneGSProb(Amean, Asd, GSGenes, Pw = rep(1, ncol(Amean)),

+  numPerm = 500, PwNull = FALSE)

 }

-

-\usage{computeGeneGSProb(Amean, Asd, GSGenes, Pw=rep(1,ncol(Amean)),numPerm=500,PwNull=F)}

-

 \arguments{

-\item{Amean}{Sampled mean value of the amplitude matrix \eqn{{\bf{A}}}.  \code{row.names(Amean)} must correspond to the gene names contained in GSGenes.}

-\item{Asd}{Sampled standard deviation of the amplitude matrix \eqn{{\bf{A}}}.}

-\item{GSGenes}{Vector containing the prior estimate of members of the gene set of interest.}

-\item{Pw}{Vector containing the weight to assign each pattern in the gene statistic assumed to be computed from the association of the pattern with samples in a given context (optional: default=1 giving all patterns equal weight).}

-\item{numPerm}{Number of permuations used for the null distribution in the gene set statistic. (optional; default=500)}

-\item{PwNull}{Logical value. If TRUE, use pattern weighting in Pw when computing the null distribution for the statistic.  If FALSE, do not use the pattern weighting so that the null is context independent. (optional; default=F)}

-}

-

-\value{

-  A vector of length GSGenes containing the p-values of set membership for each gene containined in the set specified in GSGenes.

-}

-

-\examples{

-\dontrun{

+\item{Amean}{A matrix mean values}

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

-# Results for GIST data in Fertig et al. (2012) #

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

+\item{Asd}{A matrix standard deviations}

-# load the data

-data('GIST_TS_20084')

-data('TFGSList')

+\item{GSGenes}{data.frame or list with gene sets}

-# define transcription factors of interest based on Ochs et al. (2009)

-TFs <- c("c.Jun", 'NF.kappaB', 'Smad4', "STAT3", "Elk.1", "c.Myc", "E2F.1", 

-         "AP.1", "CREB", "FOXO", "p53", "Sp1")

+\item{Pw}{weight on genes}

-# run the GAPS matrix factorization

-nIter <- 10000

-results <- CoGAPS(GIST.D, GIST.S, tf2ugFC,

-                  nFactor=5,

-                  nEquil=nIter, nSample=nIter,

-                  plot=FALSE)

+\item{numPerm}{number of permutations for null}

-# set membership statistics

-permTFStats <- list()

-for (tf in TFs) {

-     genes <- levels(tf2ugFC[,tf])

-     genes <- genes[2:length(genes)]

-     permTFStats[[tf]] <- computeGeneTFProb(Amean = GISTResults$Amean, 

-                                            Asd = GistResults$Asd, genes) 

+\item{PwNull}{- logical indicating gene adjustment}

 }

-

+\value{

+A vector of length GSGenes containing the p-values of set membership

+ for each gene containined in the set specified in GSGenes.

 }

+\description{

+Compute Gene Probability

 }

-

-\author{Elana J. Fertig \email{ejfertig@jhmi.edu}}

-

-\references{

-E.J. Fertig, A.V. Favorov, and M.F. Ochs (2013) Identifying context-specific transcription factor targets from prior knowledge and gene expression data. 2012 IEEE Nanobiosciences. 

+\details{

+Computes the p-value for gene set membership using the CoGAPS-based

+ statistics developed in Fertig et al. (2012).  This statistic refines set

+ membership for each candidate gene in a set specified in \code{GSGenes} by

+ comparing the inferred activity of that gene to the average activity of the

+ set.

 }

-\seealso{\code{\link{calcCoGAPSStat}}}

-\keyword{misc}

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

 \alias{gapsMapRun}

 \title{Backwards Compatibility with v2}

 \usage{

-gapsMapRun(D, S, FP, ...)

+gapsMapRun(D, S, FP, ABins = data.frame(), PBins = data.frame(),

+  nFactor = 5, simulation_id = "simulation", nEquil = 1000,

+  nSample = 1000, nOutR = 1000, output_atomic = FALSE,

+  fixedMatrix = "P", fixedBinProbs = FALSE, fixedDomain = "N",

+  sampleSnapshots = TRUE, numSnapshots = 100, alphaA = 0.01,

+  nMaxA = 1e+05, max_gibbmass_paraA = 100, alphaP = 0.01, nMaxP = 1e+05,

+  max_gibbmass_paraP = 100, seed = -1, messages = TRUE)

 }

 \arguments{

 \item{D}{data matrix}

@@ -4,14 +4,17 @@

 \alias{gapsRun}

 \title{Backwards Compatibility with v2}

 \usage{

-gapsRun(D, S, ...)

+gapsRun(D, S, ABins = data.frame(), PBins = data.frame(), nFactor = 7,

+  simulation_id = "simulation", nEquil = 1000, nSample = 1000,

+  nOutR = 1000, output_atomic = FALSE, fixedBinProbs = FALSE,

+  fixedDomain = "N", sampleSnapshots = TRUE, numSnapshots = 100,

+  alphaA = 0.01, nMaxA = 1e+05, max_gibbmass_paraA = 100, alphaP = 0.01,

+  nMaxP = 1e+05, max_gibbmass_paraP = 100, seed = -1, messages = TRUE)

 }

 \arguments{

 \item{D}{data matrix}

 \item{S}{uncertainty matrix}

-

-\item{...}{v2 style parameters}

 }

 \value{

 list with A and P matrix estimates

@@ -102,7 +102,7 @@ const float *AP, const float *other, float delta)

     const gaps::simd::packedFloat pDelta = delta, two = 2.f;

     gaps::simd::packedFloat d, pOther, pD, pAP, pS, partialSum = 0.f;

     gaps::simd::Index i = 0;

-    /*for (; i <= size - i.increment(); ++i)

+    for (; i <= size - i.increment(); ++i)

     {   

         pOther.load(other + i);

         pD.load(D + i);

@@ -110,7 +110,7 @@ const float *AP, const float *other, float delta)

         pS.load(S + i);

         d = pDelta * pOther;

         partialSum += (d * (two * (pD - pAP) - d)) / (two * pS * pS);

-    }*/

+    }

     float fd, delLL = partialSum.scalar();

     for (unsigned j = i.value(); j < size; ++j)

     {

@@ -127,7 +127,7 @@ float delta2)

     const gaps::simd::packedFloat pDelta1 = delta1, pDelta2 = delta2, two = 2.f;

     gaps::simd::packedFloat d, pOther1, pOther2, pD, pAP, pS, partialSum = 0.f;

     gaps::simd::Index i = 0;

-    /*for (; i <= size - i.increment(); ++i)

+    for (; i <= size - i.increment(); ++i)

     {   

         pOther1.load(other1 + i);

         pOther2.load(other2 + i);

@@ -136,7 +136,7 @@ float delta2)

         pS.load(S + i);

         d = pDelta1 * pOther1 + pDelta2 * pOther2;

         partialSum += (d * (two * (pD - pAP) - d)) / (two * pS * pS);

-    }*/

+    }

     float fd, delLL = partialSum.scalar();

     for (unsigned j = i.value(); j < size; ++j)

     {

@@ -1,11 +1,16 @@

 #ifndef __COGAPS_ARCHIVE_H__

 #define __COGAPS_ARCHIVE_H__

+#include <Rcpp.h>

 #include <fstream>

+// flags for opening an archive

 #define ARCHIVE_READ  std::ios::in

 #define ARCHIVE_WRITE std::ios::out | std::ios::trunc

+// magic number written to beginning of archive files

+#define ARCHIVE_MAGIC_NUM 0xCE45D32A

+

 class Archive

 {

 private:

@@ -16,7 +21,21 @@ public:

     Archive(const std::string &path, std::ios_base::openmode flags)

         : mStream(path.c_str(), std::ios::binary | flags)

-    {}

+    {

+        /*if (flags == ARCHIVE_WRITE)

+        {

+            *this << ARCHIVE_MAGIC_NUM;

+        }

+        else if (flags == ARCHIVE_READ)

+        {

+            uint32_t magicNum = 0;

+            *this >> magicNum;

+            if (magicNum != ARCHIVE_MAGIC_NUM)

+            {

+                Rcpp::Rcout << "warning: invalid checkpoint file" << std::endl;

+            }

+        }*/

+    }

     void close() {mStream.close();}

@@ -16,13 +16,16 @@

 #define SSTR(x) static_cast<std::ostringstream&>( \

     (std::ostringstream() << std::dec << x)).str()

+// used to convert defined macro values into strings

 #define STR_HELPER(x) #x

 #define STR(x) STR_HELPER(x)

-#define ARCHIVE_MAGIC_NUM 0xCE45D32A

-

+// boost time helpers

 namespace bpt = boost::posix_time;

-static bpt::ptime lastCheckpoint; // keep track of when checkpoints are made

+#define bpt_now() bpt::microsec_clock::local_time()

+

+// keeps track of when checkpoints are made

+static bpt::ptime lastCheckpoint; 

 // save the current internal state to a file

 static void createCheckpoint(GapsInternalState &state)

@@ -31,82 +34,112 @@ static void createCheckpoint(GapsInternalState &state)

     state.numCheckpoints++;

     // record starting time

-    bpt::ptime start = bpt::microsec_clock::local_time();

+    bpt::ptime start = bpt_now();

     // save state to file, write magic number at beginning

-    Archive ar("gaps_checkpoint_" + SSTR(state.numCheckpoints) + ".out",

-        ARCHIVE_WRITE);

-    ar << ARCHIVE_MAGIC_NUM;

+    std::string fname("gaps_checkpoint_" + SSTR(state.numCheckpoints) + ".out");

+    Archive ar(fname, ARCHIVE_WRITE);

     gaps::random::save(ar);

     ar << state.sampler.nFactor() << state.nEquil << state.nSample << state;

     ar.close();

     // display time it took to create checkpoint

-    bpt::time_duration diff = bpt::microsec_clock::local_time() - start;

+    bpt::time_duration diff = bpt_now() - start;

     double elapsed = diff.total_milliseconds() / 1000.;

     Rcpp::Rcout << " finished in " << elapsed << " seconds\n";

 }

-static void runGibbsSampler(GapsInternalState &state, unsigned nIterTotal,

-Vector &chi2Vec, Vector &aAtomVec, Vector &pAtomVec)

+static void updateSampler(GapsInternalState &state)

+{

+    state.nUpdatesA += state.nIterA;

+    for (unsigned j = 0; j < state.nIterA; ++j)

+    {

+        state.sampler.update('A');

+    }

+

+    state.nUpdatesP += state.nIterP;

+    for (unsigned j = 0; j < state.nIterP; ++j)

+    {

+        state.sampler.update('P');

+    }

+}

+

+static void makeCheckpointIfNeeded(GapsInternalState &state)

+{

+    bpt::time_duration diff = bpt_now() - lastCheckpoint;

+    int diff_sec = diff.total_milliseconds() / 1000;

+    if (diff_sec > state.checkpointInterval && state.checkpointInterval > 0)

+    {

+        createCheckpoint(state);

+        lastCheckpoint = bpt_now();

+    }

+}

+

+static void storeSamplerInfo(GapsInternalState &state, Vector &atomsA,

+Vector &atomsP, Vector &chi2)

+{

+    chi2[state.iter] = state.sampler.chi2();

+    atomsA[state.iter] = state.sampler.totalNumAtoms('A');

+    atomsP[state.iter] = state.sampler.totalNumAtoms('P');

+    state.nIterA = gaps::random::poisson(std::max(atomsA[state.iter], 10.f));

+    state.nIterP = gaps::random::poisson(std::max(atomsP[state.iter], 10.f));

+}

+

+static void displayStatus(GapsInternalState &state, const std::string &type,

+unsigned nIterTotal)

+{

+    if ((state.iter + 1) % state.nOutputs == 0 && state.messages)

+    {

+        Rcpp::Rcout << type << state.iter + 1 << " of " << nIterTotal

+            << ", Atoms:" << state.sampler.totalNumAtoms('A') << "("

+            << state.sampler.totalNumAtoms('P') << ") Chi2 = "

+            << state.sampler.chi2() << '\n';

+    }

+}

+

+static void takeSnapshots(GapsInternalState &state)

+{

+    if (state.nSnapshots && !((state.iter+1)%(state.nSample/state.nSnapshots)))

+    {

+        state.snapshotsA.push_back(state.sampler.getNormedMatrix('A'));

+        state.snapshotsP.push_back(state.sampler.getNormedMatrix('P'));

+    }    

+}

+

+static void runBurnPhase(GapsInternalState &state)

+{

+    for (; state.iter < state.nEquil; ++state.iter)

+    {

+        makeCheckpointIfNeeded(state);

+        float temp = ((float)state.iter + 2.f) / ((float)state.nEquil * 2.f);

+        state.sampler.setAnnealingTemp(std::min(1.f,temp));

+        updateSampler(state);

+        displayStatus(state, "Equil: ", state.nEquil);

+        storeSamplerInfo(state, state.nAtomsAEquil, state.nAtomsPEquil,

+            state.chi2VecEquil);

+    }

+}

+

+static void runCoolPhase(GapsInternalState &state)

+{

+    for (; state.iter < state.nEquilCool; ++state.iter)

+    {

+        makeCheckpointIfNeeded(state);

+        updateSampler(state);

+    }

+}

+

+static void runSampPhase(GapsInternalState &state)

 {

-    for (; state.iter < nIterTotal; ++state.iter)

+    for (; state.iter < state.nSample; ++state.iter)

     {

-        bpt::ptime now = bpt::microsec_clock::local_time();

-        bpt::time_duration diff = now - lastCheckpoint;

-        if (diff.total_milliseconds() > state.checkpointInterval * 1000

-        && state.checkpointInterval > 0)

-        {

-            createCheckpoint(state);

-            lastCheckpoint = bpt::microsec_clock::local_time();

-        }

-

-        if (state.phase == GAPS_BURN)

-        {

-            state.sampler.setAnnealingTemp(std::min(1.f,

-                ((float)state.iter + 2.f) / ((float)state.nEquil / 2.f)));

-        }

-        

-        state.nUpdatesA += state.nIterA;

-        state.nUpdatesP += state.nIterP;

-        for (unsigned j = 0; j < state.nIterA; ++j)

-        {

-            state.sampler.update('A');

-        }

-        for (unsigned j = 0; j < state.nIterP; ++j)

-        {

-            state.sampler.update('P');

-        }

-

-        if (state.phase == GAPS_SAMP)

-        {

-            state.sampler.updateStatistics();

-            if (state.nSnapshots && !((state.iter + 1) %

-            (nIterTotal / state.nSnapshots)))

-            {

-                state.snapshotsA.push_back(state.sampler.getNormedMatrix('A'));

-                state.snapshotsP.push_back(state.sampler.getNormedMatrix('P'));

-            }

-        }

-

-        if (state.phase != GAPS_COOL)

-        {

-            float nAtomsA = state.sampler.totalNumAtoms('A');

-            float nAtomsP = state.sampler.totalNumAtoms('P');

-            aAtomVec[state.iter] = nAtomsA;

-            pAtomVec[state.iter] = nAtomsP;

-            state.nIterA = gaps::random::poisson(std::max(nAtomsA, 10.f));

-            state.nIterP = gaps::random::poisson(std::max(nAtomsP, 10.f));

-

-            if ((state.iter + 1) % state.nOutputs == 0 && state.messages)

-            {

-                std::string ph(state.phase == GAPS_BURN ? "Equil: " : "Samp: ");

-                Rcpp::Rcout << ph << state.iter + 1 << " of " << nIterTotal

-                    << ", Atoms:" << aAtomVec[state.iter] << "("

-                    << pAtomVec[state.iter] << ") Chi2 = "

-                    << state.sampler.chi2() << '\n';

-            }

-        }

+        makeCheckpointIfNeeded(state);

+        updateSampler(state);

+        state.sampler.updateStatistics();

+        takeSnapshots(state);

+        displayStatus(state, "Samp: ", state.nSample);

+        storeSamplerInfo(state, state.nAtomsASample, state.nAtomsPSample,

+            state.chi2VecSample);

     }

 }

@@ -114,25 +147,24 @@ Vector &chi2Vec, Vector &aAtomVec, Vector &pAtomVec)

 static Rcpp::List runCogaps(GapsInternalState &state)

 {

     // reset the checkpoint timer

-    lastCheckpoint = bpt::microsec_clock::local_time();

+    lastCheckpoint = bpt_now();

     // cascade down the various phases of the algorithm

     // this allows for starting in the middle of the algorithm

-    Vector trash(1);

     switch (state.phase)

     {

         case GAPS_BURN:

-            runGibbsSampler(state, state.nEquil, state.chi2VecEquil,

-                state.nAtomsAEquil, state.nAtomsPEquil);

+            runBurnPhase(state);

             state.iter = 0;

             state.phase = GAPS_COOL;

+

         case GAPS_COOL:

-            runGibbsSampler(state, state.nEquilCool, trash, trash, trash);

+            runCoolPhase(state);

             state.iter = 0;

             state.phase = GAPS_SAMP;

+

         case GAPS_SAMP:

-            runGibbsSampler(state, state.nSample, state.chi2VecSample,

-                state.nAtomsASample, state.nAtomsPSample);

+            runSampPhase(state);

     }

     // combine chi2 vectors

@@ -165,43 +197,36 @@ static Rcpp::List runCogaps(GapsInternalState &state)

 }

 // [[Rcpp::export]]

-Rcpp::List cogaps_cpp(const Rcpp::NumericMatric &DMatrix,

-const Rcpp::NumericMatrix &SMatrix, const Rcpp::NumericMatrix &fixedPatterns,

-const Rcpp::S4 &params)

+Rcpp::List cogaps_cpp(const Rcpp::NumericMatrix &D,

+const Rcpp::NumericMatrix &S, unsigned nFactor, unsigned nEquil, unsigned nEquilCool,

+unsigned nSample, unsigned nOutputs, unsigned nSnapshots, float alphaA,

+float alphaP, float maxGibbmassA, float maxGibbmassP, int seed, bool messages,

+bool singleCellRNASeq, char whichMatrixFixed, const Rcpp::NumericMatrix &FP,

+unsigned checkpointInterval)

 {

     // set seed

     uint32_t seedUsed = static_cast<uint32_t>(seed);

     if (seed < 0)