diff --git a/.Rbuildignore b/.Rbuildignore
new file mode 100644
index 0000000..91114bf
--- /dev/null
+++ b/.Rbuildignore
@@ -0,0 +1,2 @@
+^.*\.Rproj$
+^\.Rproj\.user$
diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..807ea25
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,3 @@
+.Rproj.user
+.Rhistory
+.RData
diff --git a/DESCRIPTION b/DESCRIPTION
new file mode 100644
index 0000000..62a9fd9
--- /dev/null
+++ b/DESCRIPTION
@@ -0,0 +1,10 @@
+Package: descend
+Title: Gene Expression Distribution Deconvolution for Single Cell RNA-seq
+Version: 0.0.0.9000
+Authors@R: person("First", "Last", email = "first.last@example.com", role = c("aut", "cre"))
+Description: What the package does (one paragraph).
+Depends: R (>= 3.3.0)
+License: What license is it under?
+Encoding: UTF-8
+LazyData: true
+RoxygenNote: 6.0.1
diff --git a/NAMESPACE b/NAMESPACE
new file mode 100644
index 0000000..901923a
--- /dev/null
+++ b/NAMESPACE
@@ -0,0 +1,8 @@
+# Generated by roxygen2: do not edit by hand
+
+export(DESCEND)
+export(deconvG)
+export(deconvSingle)
+export(descend)
+export(getEstimates)
+exportClasses(DESCEND)
diff --git a/R/deconvSingle.R b/R/deconvSingle.R
new file mode 100644
index 0000000..2fa7c90
--- /dev/null
+++ b/R/deconvSingle.R
@@ -0,0 +1,410 @@
+#' Deconvolve the true gene expression distribution  of a single gene
+#'
+#' The deconvolution is computed by using the function \code{deconvG}. This function can automatically discretize the underlying distribution and find the proper tuning parameter of the penalty term in G-modelling. Besides, it computs the estimates and standard deviations of five distribution based statistics (active fraction, active intensity, mean, CV and gini coefficient), as well as the estimated coefficients of the covariates on active intensity (Z) and active fraction (Z0).
+#'
+#' @param y a vector of observed counts of a single gene
+#' @param scaling.consts a vector of cell specific scaling constants, either the cell efficiency or the library size
+#' @inheritParams deconvG
+#' @param whether plot the density curve of the deconvolved the distribution or not. The zero inflation part has been smoothed into the density curve for visualization. Default is True.
+#' @param do.LRT.test whether do LRT test on the coefficients and active fraction or not. Default is True
+#' @param verbose verbose the estimation and testing procedures or not. Default is True.
+#' @param control settings see see {\code{\link{deconv.control}}}
+#'
+#' @return a DESCEND object
+#' @export
+
+
+deconvSingle <- function(y, 
+                         scaling.consts = NULL, 
+                         Z = NULL,
+                         Z0 = NULL,
+                         plot.density = T,
+                         do.LRT.test = T,
+                         family = c("Poisson", "Negative Binomial"),
+                         NB.size = 100,
+                         verbose = T, control = list()) {
+  family <- match.arg(family)
+  control <- do.call("DESCEND.control", control)
+
+  if (mean(y == 0) > control$max.sparse)
+    stop("Too sparse data!")
+
+  if (is.null(scaling.consts)) {
+    offset <- rep(0, length(y))
+  } else {
+    if (min(scaling.consts) <= 0)
+      stop("Cell specific scaling constants should be positive!")
+    offset <- log(scaling.consts)
+  }
+
+
+  lam.max <- quantile(y/exp(offset), probs = control$max.quantile)
+  if (control$only.integer) {
+    if (lam.max <= control$n.points)
+      tau <- log(1:ceiling(lam.max))
+    else 
+      tau <- c(log(1:5), log(seq(5, lam.max, length.out = control$n.points - 5 + 1)[-1]))
+  } else {
+    lam.min <- lam.max/(2 * control$n.points - 1)
+    tau <- log(seq(lam.min, lam.max, length.out = control$n.points))
+  }
+
+  ## rescale Z to guarantee that it does not affect the range of theta, will scale back later
+  if (!is.null(Z)) {
+    Z <- as.matrix(Z)
+    Z.means <- log(colMeans(exp(Z)))
+    Z <- t(t(Z) - Z.means)
+  }
+
+  if (!is.null(Z0)) 
+    Z0 <- as.matrix(Z0)
+
+  result <- list()
+  value <- Inf
+  c0.start <- control$c0.start
+
+  if (!control$zeroInflate)
+    p <- control$pDegree
+  else {
+    p <- control$pDegree + 1
+    if (!is.null(Z0))
+      p <- p + ncol(Z0)
+  }
+
+  for (i in 1:control$nStart) {
+    if (verbose)
+      print(paste("Compute penalized MLE, the ", i, "th time", sep = ""))
+    aStart <- rnorm(p, control$aStart, control$start.sd)
+    if (!is.null(Z0))
+      bStart <- rnorm(ncol(Z0), control$bStart, control$start.sd)
+    if (!is.null(Z))
+      gStart <- rnorm(ncol(Z), control$gStart, control$start.sd)
+
+    if (verbose)
+      print(paste("c0 is ", c0.start))
+    temp <- try(deconvG(tau, y, offset = offset,
+                        Z = Z, Z0 = Z0,
+                        family = family, NB.size = NB.size,
+                        zeroInflate = control$zeroInflate,
+                        c0 = c0.start,
+                        pDegree = control$pDegree,
+                        aStart = aStart, bStart = bStart, gStart = gStart))
+
+    if (verbose)
+      print(paste("Relative fake information is:", round(temp$S * 100, 2), "%", sep = ""))
+
+    if (class(temp) == "try-error" || is.na(temp$S))
+      return(NA)
+   
+    iter <- 0
+    while (1) {
+      if (iter == control$max.c0.iter) {
+        temp <- list(value = Inf)
+        if (verbose)
+          warning("Fail to find a proper c0!")
+        break
+      }
+
+      if (temp$S < control$rel.info.range[2] & temp$S > control$rel.info.range[1])
+        break
+      scale.c0 <- control$rel.info.guide/temp$S
+      scale.c0 <- min(scale.c0, 100)
+      scale.c0 <- max(scale.c0, 0.01)
+      c0.start <- max(scale.c0 * c0.start, control$c0.min)
+      if (verbose)
+        print(paste("c0 is ", c0.start))
+      temp <- try(deconvG(tau, y, offset = offset,
+                            Z = Z, Z0 = Z0,
+                            family = family, NB.size = NB.size,
+                            zeroInflate = control$zeroInflate,
+                            c0 = c0.start,
+                            pDegree = control$pDegree,
+                            aStart = aStart, bStart = bStart, gStart = gStart))
+      if (verbose)
+        print(paste("Relative fake information is:", round(temp$S * 100, 2), "%", sep = ""))
+      if (class(temp) == "try-error")
+        return(NA)
+      iter <- iter + 1
+    }
+    if (temp$value < value) {
+      value <- temp$value
+      result <- temp
+    }
+  }
+
+  if (!is.finite(value))
+    return(NA)
+
+  ## collect estimation of statistics
+  theta <- result$stats$mat.dist[, 1]
+  g <- result$stats$mat.dist[, 2]
+  bias.g <- result$stats$mat.dist[, 6]
+
+  mu <- c(sum(theta * g), sum(theta * bias.g),
+          sqrt(t(theta) %*% result$cov.g %*% theta))
+  sd <- sqrt(sum(theta^2 * g) - mu[1]^2)
+  kappa <- sum(theta^2 * g)
+  CV2.devg <- theta^2  / mu[1]^2 - 2 * kappa * theta / mu[1]^3
+  CV <- sd/mu[1]
+  CV.devg <- CV2.devg / (2 * CV)
+  CV <- c(CV, sum(CV.devg * bias.g), sqrt(t(CV.devg) %*% result$cov.g %*% CV.devg))
+  gamma.devg <- sapply(1:length(g), function(i)sum(abs(theta[i] - theta) * g))
+  gamma <- sum(g * gamma.devg) /2
+  gini.devg <- gamma.devg / mu[1] - gamma * theta/mu[1]^2
+  gini <- c(gamma/mu[1], sum(gini.devg * bias.g), sqrt(t(gini.devg) %*% result$cov.g %*% gini.devg))
+  if (control$zeroInflate) {
+    p0 <- c(result$stats$mat.dist[1, c(2, 6, 3)])
+    if (!is.null(Z)) {
+      coefs <- result$stats$mat.coef[1:ncol(Z), , drop = F]
+      idx <- (nrow(result$cov) - ncol(Z) + 1):nrow(result$cov)
+      cov.coef <- result$cov[idx, idx, drop = F]
+      mu.pos.val <- sum(theta[-1] * g[-1]/(1 - g[1]))
+      mu.pos.adj <- exp(-sum(coefs[, 1] * Z.means))
+      mu.pos.dev <- theta[-1]/(1-g[1]) - mu[1]/(1-g[1])^2
+      mu.pos.adj.dev <- -mu.pos.adj * Z.means
+      mu <- c(mu[1] * mu.pos.adj, 
+              mu.pos.adj * mu[2] + mu[1] * sum(mu.pos.adj.dev * coefs[, 2]),
+              sqrt(mu.pos.adj^2 * t(theta) %*% result$cov.g %*% theta+ 
+                   2 * mu.pos.adj * mu[1] * t(theta) %*% 
+                   result$cov.g.gamma %*% mu.pos.adj.dev + 
+                   mu[1]^2 * t(mu.pos.adj.dev) %*% cov.coef %*% mu.pos.adj.dev))
+      mu.pos <- c(mu.pos.val * mu.pos.adj, 
+                  mu.pos.adj * sum(mu.pos.dev * bias.g[-1]) + mu.pos.val * sum(mu.pos.adj.dev * coefs[, 2]),
+                  sqrt(mu.pos.adj^2 * t(mu.pos.dev) %*% result$cov.g[-1, -1] %*% mu.pos.dev+ 
+                       2 * mu.pos.adj * mu.pos.val * t(mu.pos.dev) %*% 
+                       result$cov.g.gamma[-1, , drop = F] %*% mu.pos.adj.dev + 
+                       mu.pos.val^2 * t(mu.pos.adj.dev) %*% cov.coef %*% mu.pos.adj.dev))
+    } else {
+      mu.pos.val <- sum(theta[-1] * g[-1]/(1 - g[1]))
+      mu.pos.dev <- theta[-1]/(1-g[1]) - mu[1]/(1-g[1])^2
+      mu.pos <- c(mu.pos.val, 
+                  sum(mu.pos.dev * bias.g[-1]), 
+                  sqrt(t(mu.pos.dev) %*% result$cov.g[-1, -1] %*% mu.pos.dev))
+
+    }
+    sd.pos <- sqrt(sum(theta[-1]^2 * g[-1]/(1 - g[1])) - mu.pos[1]^2)
+  } else {
+    p0 <- rep(0, 3)
+    mu.pos <- mu
+  }
+
+  one.minus.p0 <- c(1 - p0[1], -p0[2], p0[3])
+
+  estimates <- rbind(one.minus.p0, mu.pos, mu, CV, 
+                      gini)
+  colnames(estimates) <- c("est", "bias", "sd")
+  rownames(estimates) <- c("Active Fraction", "Active Intensity", "Mean", "CV", "Gini")
+  if (!is.null(result$stats$mat.coef))
+    estimates <- rbind(estimates, result$stats$mat.coef)
+
+  if (!control$only.integer) {
+    if (control$zeroInflate) {
+      theta <- theta[-1]
+      g1 <- g
+      if (control$dense.add.0) {
+        g1[2] <- g1[1] + g1[2]
+        g1 <- g1[-1]
+      } else {
+        g1 <- g1[-1]/sum(g1[-1])
+      }
+    } else
+      g1 <- g
+    density.points <- cbind(theta = theta, density = g1/(theta[2] - theta[1]))
+  } else
+    density.points <- cbind(theta, probability = g)
+
+  if (plot.density) {
+    if (!control$only.integer)
+      plot(density.points, pch = 20, type = "l", ylab = "density", 
+           main = "Deconvolved Distribution")
+    else
+      plot(density.points, pch = 20, ylab = "density", 
+           main = "Deconvolved Distribution")
+  }
+
+  if (do.LRT.test) {   
+    LRT.Z.select <- control$LRT.Z.select
+    LRT.Z0.select <- control$LRT.Z0.select
+    LRT.Z.values <- control$LRT.Z.values
+    condition.list <- list(list(control$zeroInflate, Z, Z0, p, offset))
+    test.name <-  c("Full model")
+    if (control$zeroInflate) {
+      if (!is.null(Z0))
+        p1 <- p - 1 - ncol(Z0)
+      else
+        p1 <- p-1
+      condition.list <- c(condition.list, p0 = list(list(!control$zeroInflate, Z, NULL, p1, offset)))
+      names(condition.list)[length(condition.list)] <- "Active Fraction = 1"
+      test.name <- c(test.name, "Active fraction = 1")
+    }
+    if (!is.null(Z)) {
+      if (length(LRT.Z.select) == 1 || length(LRT.Z.select) == ncol(Z) &&
+          (length(LRT.Z.values) == 1 || length(LRT.Z.values) == ncol(Z))) {
+        if (length(LRT.Z.select) == 1)
+          LRT.Z.select <- rep(LRT.Z.select, ncol(Z))
+        if (length(LRT.Z.values) == 1)
+          LRT.Z.values <- rep(LRT.Z.values, ncol(Z))
+        if (ncol(Z) == 1 && LRT.Z.select) {
+          offset1 <- offset + LRT.Z.values * Z
+          condition.list <- c(condition.list, Z = list(list(control$zeroInflate, NULL, Z0, p, offset1)))
+          names(condition.list)[length(condition.list)] <- paste("Effect of Z: gamma = ",
+                                                                 LRT.Z.values[i], sep = "")
+          test.name <- c(test.name, paste("Effect of Z: gamma = ", LRT.Z.values[i], sep = ""))
+        } else
+          for(i in 1:ncol(Z)) {
+            if (LRT.Z.select[i]) {
+              offset1 <- offset + LRT.Z.values[i] * Z[, i]
+              condition.list <- c(condition.list, 
+                                  list(list(control$zeroInflate, Z[, -i, drop = F], Z0, p, offset1)))
+              names(condition.list)[length(condition.list)] <- paste("Effect of Z: gamma", i, 
+                                                                      " = ", LRT.Z.values[i], sep = "")
+              test.name <- c(test.name, paste("Effect of Z: gamma", i, " = ", LRT.Z.values[i], sep = ""))
+            }
+          }
+      } else 
+        warning("length of LRT.Z.select should be either 1 or the same as 
+                ncol(Z), LRT test on Z not performed!")
+    }
+    if (!is.null(Z0) && control$zeroInflate) {
+      if (length(LRT.Z0.select) == 1 || length(LRT.Z0.select) == ncol(Z0)) {
+        if (length(LRT.Z0.select) == 1)
+        LRT.Z0.select <- rep(LRT.Z0.select, ncol(Z0))
+        if (ncol(Z0) == 1 && LRT.Z0.select) {
+          condition.list <- c(condition.list, Z0 = list(list(control$zeroInflate, Z, NULL, p - 1, offset)))
+          names(condition.list)[length(condition.list)] <- "Effect of Z0: beta = 0"
+          test.name <- c(test.name, "Effect of Z0: beta = 0")
+        } else
+          for(i in 1:ncol(Z0)) {
+            if (LRT.Z0.select[i]) {
+              condition.list <- c(condition.list, 
+                                  list(list(control$zeroInflate, Z, Z0[, -i, drop = F], p-1, offset)))
+              names(condition.list)[length(condition.list)] <- paste("Effect of Z0: beta", i, "= 0", sep = "")
+              test.name <- c(test.name, paste("Effect of Z0: beta", i, " = 0", sep = ""))
+            }
+          }
+      } else 
+        warning("length of LRT.Z.select should be either 1 or the same as 
+                ncol(Z), LRT test on Z0 not performed!")
+    }
+ #   print(condition.list)
+
+
+    if (length(condition.list) > 1) {
+      if (verbose)
+        print(paste("Compute LRT tests, in total", length(condition.list) - 1, "tests"))
+      lk <- lapply(1:length(condition.list), function(k) {
+                   if (verbose)
+                     print(paste("Compute LRT test, ", test.name[k], sep = ""))
+                   ll <- condition.list[[k]]
+                   val <- Inf
+                   for (i in 1:control$nStart.lrt) {
+                     aStart <- rnorm(ll[[4]], control$aStart, control$start.sd)
+                     if (!is.null(ll[[3]]))
+                       bStart <- rnorm(ncol(ll[[3]]), control$bStart, control$start.sd)
+                     if (!is.null(ll[[2]]))
+                       gStart <- rnorm(ncol(ll[[2]]), control$gStart, control$start.sd)
+
+                     temp <- try(deconvG(tau, y, offset = ll[[5]],
+                                           Z = ll[[2]], Z0 = ll[[3]],
+                                           family = family, NB.size = NB.size,
+                                           zeroInflate = ll[[1]],
+                                           c0 = 0,
+                                           pDegree = control$pDegree,
+                                           aStart = aStart, bStart = bStart, gStart = gStart,
+                                           only.value = T))
+                     if (class(temp) == "try-error")
+                       temp <- Inf
+                     if (temp < val)
+                       val <- temp
+                   }
+                   df <- max(1, p-ll[[4]])
+                   if (k == 1)
+                     df <- NA
+                   if (verbose)
+                     print(paste("neg loglikehod is ", as.numeric(val), 
+                                 "df is ", df))
+                   return(c(as.numeric(val), max(1, p - ll[[4]])))
+    })
+      if (is.finite(lk[[1]][1])) {
+        pval <- sapply(lk[-1], function(vec) {
+                       if (is.finite(vec[1]))
+                          return(pchisq(2 * (vec[1] - lk[[1]][1]), 
+                                        vec[2], lower.tail = F))
+                        return(NA)
+                      })
+        pval <- as.matrix(pval, ncol = 1)
+        colnames(pval) <- "P.value"
+        rownames(pval) <- names(condition.list)[-1]
+      } else
+        pval <- NA
+    } else
+      pval <- NA
+  } else {
+    pval <- NA
+  }
+
+  estimates <- cbind(estimates, DESCEND.sd = sqrt(estimates[, 2]^2 + estimates[, 3]^2))
+
+  result <- new("DESCEND", 
+                distribution = result$stats$mat.dist,
+                estimates = estimates,
+                density.points = density.points)
+  if (!is.na(pval))
+    result@pval <- pval
+
+  return(result)
+}
+
+#' The control parameters of the function \code{\link{descend}} and \code{\link{deconvSingle}}
+DESCEND.control <- function(n.points = 50,
+                            nStart = 2,
+                            nStart.lrt = 2,
+                            max.quantile = 0.98,
+                            max.sparse = 0.95,
+                            rel.info.range = c(0.0005, 0.01),
+                            rel.info.guide = 0.005,
+                            c0.start = 1,
+                            aStart = 1,
+                            bStart = 0,
+                            gStart = 0, 
+                            start.sd = 0.2, 
+                            penalty.Z0 = T,
+                            pDegree = 5,
+                            max.c0.iter = 5,
+                            c0.min = 1e-5,
+                            LRT.Z.select = T,
+                            LRT.Z0.select = T,
+                            LRT.Z.values = 0,
+                            zeroInflate = T,
+                            dense.add.0 = T,
+                            only.integer = F) {
+
+
+  n.points <- max(n.points, 10)
+
+  list(n.points = n.points, nStart = nStart, nStart.lrt = nStart.lrt, 
+       max.quantile = max.quantile,
+       max.sparse = max.sparse,
+       rel.info.range = rel.info.range, rel.info.guide = rel.info.guide,
+       c0.start = c0.start, 
+       aStart = aStart, bStart = bStart, gStart = gStart, start.sd = start.sd,
+       penalty.Z0 = penalty.Z0, pDegree = pDegree, max.c0.iter = max.c0.iter,
+       c0.min = c0.min,
+       LRT.Z.select = LRT.Z.select, 
+       LRT.Z0.select = LRT.Z0.select,
+       LRT.Z.values = LRT.Z.values,
+       zeroInflate = zeroInflate, 
+       dense.add.0 = dense.add.0,
+       only.integer = only.integer)
+}
+
+
+#' @export DESCEND
+#' @exportClass DESCEND
+
+require(methods)
+
+setClass("DESCEND", representation(distribution = "matrix", 
+                                   estimates = "matrix", 
+                                   pval = "matrix", 
+                                   density.points = "matrix"))
diff --git a/R/descend.R b/R/descend.R
new file mode 100644
index 0000000..f8550bd
--- /dev/null
+++ b/R/descend.R
@@ -0,0 +1,163 @@
+#' Default is to set max.quantile = max.sparse = 0.98. For those that are too sparse, set max.quantile = max.sparse = 0.99 and try again.
+#'
+#' @inheritParams deconvSingle
+#' @param count.matrix the observed UMI count matrix. Each row is a gene and each column is a cell
+#' @param n.cores the number of cores used for parallel computing. Default is 1. Used only when parallel computing is done in a single machine. For using multi-machine cores, need to assign \code{cl} explicitly
+#' @param cl an object of class "cluster". See more details in \code{\link[parallel]{makeCluster}}
+#' 
+#' @return a list of DESCEND objects
+#'
+#' @export
+
+descend <- function(count.matrix,
+                    scaling.consts = NULL,
+                    Z = NULL,
+                    Z0 = NULL,
+                    n.cores = 1, cl = NULL,
+                    do.LRT.test = F, 
+                    family = c("Poisson", "Negative Binomial"),
+                    NB.size = 100,
+                    verbose = T, 
+                    control = list()) {
+
+   control <- do.call("DESCEND.control", control)
+
+  Y <- as.matrix(count.matrix)
+  gene.names <- rownames(Y)
+
+  print(paste("DESCEND starts deconvolving distribution of", length(gene.names), "genes!"))
+
+  if (is.null(cl)) {
+    require(doParallel)
+    registerDoParallel(n.cores)
+
+    results <-  foreach(v = iter(Y, "row"), .noexport = c("Y")) %dopar% {
+      if (verbose)
+        print("Compute for gene")
+      if (mean(v == 0) > control$max.sparse) {
+        control$max.quantile <- 0.99
+        control$max.sparse <- 0.99
+      }
+
+      temp <- try(deconvSingle(v,
+                               scaling.consts = scaling.consts,
+                               Z = Z, Z0 = Z0,
+                               plot.density = F,
+                               do.LRT.test = do.LRT.test,
+                               family = family, 
+                               NB.size = NB.size,
+                               control = control, verbose = F))
+      if (class(temp) == "try-error")
+        return(NA)
+      return(temp)
+    }
+
+  } else if (!is.null(cl)) {
+    Y <- as.list(as.data.frame(t(Y)))
+    clusterExport(cl, c("scaling.consts", "Z", "Z0", "do.LRT.test", "family", "NB.size",
+                        "control"), environment())
+
+    results <- parLapply(cl, Y, function(v) {
+
+                         source("~/Dropbox/sparse_factor_bic/code/g_model/package_code/deconvSingle.R")
+                         source("~/Dropbox/sparse_factor_bic/code/g_model/package_code/g_model.R")
+
+                         if ("logMsg" %in% ls())
+                           logMsg(do.LRT.test)
+
+                         if (mean(v == 0) > control$max.sparse) {
+                           control$max.quantile <- 0.99
+                           control$max.sparse <- 0.99
+                         }
+
+                         temp <- try(deconvSingle(v,
+                                                  scaling.consts = scaling.consts,
+                                                  Z = Z, Z0 = Z0,
+                                                  plot.density = F,
+                                                  do.LRT.test = do.LRT.test,
+                                                  family = family, 
+                                                  NB.size = NB.size,
+                                                  control = control, verbose = F))
+                         if (class(temp) == "try-error")
+                           return(NA)
+                         return(temp)
+                    })
+  } 
+  names(results) <- gene.names
+
+  return(results)
+}
+
+#' Grep the value and standard deviation of the estimated parameters calculated from a list of descend objects 
+#'
+#' @param descend.list a list of descend objects
+#' 
+#' @export
+
+getEstimates <- function(descend.list) {
+  temp <- lapply(descend.list, function(ll) {
+                      if (class(ll) == "DESCEND")
+                        return(ll@estimates[, c(1, 4)])
+                      else 
+                        return(NA)
+                    })
+  idx <- which(!sapply(temp, function(mat)is.na(mat[1])))[1]
+  n.est <- nrow(temp[[idx]])
+  est.names <- rownames(temp[[idx]])
+  report.names <- colnames(temp[[idx]])
+
+  temp <- lapply(temp, function(mat) {
+                      if (is.na(mat[1])) {
+                        mm <- matrix(NA, n.est, 2)
+                        colnames(mm) <- report.names
+                        rownames(mm) <- est.names
+                        return(mm)
+                      }
+                      else
+                        return(mat)
+                    })
+
+  est.list <- lapply(1:length(est.names), function(i) {
+                      result <- t(sapply(temp, function(mat) mat[i, ]))
+                      colnames(result) <- report.names
+                      rownames(result) <- names(descend.list)
+                      return(result)
+                    })
+
+  names(est.list) <- est.names
+  return(est.list)
+}
+
+getPval <- function(descend.list) {
+  temp <- lapply(descend.list, function(ll) {
+                      if (class(ll) == "DESCEND")
+                        return(ll@pval)
+                      else 
+                        return(NA)
+                    })
+
+  idx <- which(!sapply(temp, is.na))[1]
+  n.test <- nrow(temp[[idx]])
+  if (n.test == 0) 
+    return(NA)
+  test.names <- rownames(temp[[idx]])
+  report.names <- colnames(temp[[idx]])
+
+  temp <- lapply(temp, function(mat) {
+                      if (is.na(mat[1])) {
+                        mm <- matrix(NA, n.est, 1)
+                        colnames(mm) <- report.names
+                        rownames(mm) <- test.names
+                      }
+                      else
+                        return(mat)
+                    })
+
+
+  test.result <- t(sapply(temp, function(mat) mat[, 1]))
+
+  rownames(test.result) <- names(descend.list)
+
+  return(test.result)
+
+}
diff --git a/R/g_model.R b/R/g_model.R
new file mode 100644
index 0000000..8a9e920
--- /dev/null
+++ b/R/g_model.R
@@ -0,0 +1,643 @@
+#' Base function of G-modeling deconvolution
+#'
+#' Base function used by \code{\link{deconvSingle}} to deconvolve the underlying distribution. 
+#' We assume X \sim F(T) where F is the noise distribution. We assume that 
+#' \deqn{log(T) = offset + \gammaZ + \epsilon}
+#' \deqn{P(T = 0) = \beta_0 + \beta_1Z0}
+#' The goal is the recover the distribution of exp(log(T) - offset - \gammaZ), which has density g and is discretized at exp(tau) (add 0 when zero inflation is allowed).
+#'
+#' @param tau log of the discrete points of the deconvolved distribution
+#' @param X a vector of observed counts 
+#' @param family family of the noise distribution, support either "Poisson" or "Negative Binomial" with known tuning parameter
+#' @param ignoreZero whether ignore the zero count. If true, then use truncated Poisson / Negative Binomial distribution. Default is False
+#' @param zeroInflate whether add zero inflation part to the deconvolved distribution to reflect transcriptional bursting. Default is True. 
+#' @param Z covariates for active intensity. Default is NULL. 
+#' @param Z0 covariates for active fraction. Used only when zeroInflate is True. Default is NULL.
+#' @param c0 the tuning parameter on the L2 penalty term. Default is 1. c0 will be selected automatically in \code{deconvSingle}
+#' @param NB.size over-dispersion parameter when the family is Negative Binomial: mu = mu + mu^2/size
+#' @param only.value whether not to compute the estimation statistics but only the value of the optimized lieklihood. Used for likelihood ratio test.
+#' @param pDegree the degree of the Spline matrix. Default is 5.
+#' @param aStart,bStart,gStart initial values of the density parameters, the coefficients of Z0 and coefficients of Z
+#' @param ... extra parameters for the \code{\link[stats]{nlm}} function
+#'
+#' @return a list with elements
+#' \item{stats}{a list of two elements. One is the \code{mat.dist}, which is the matrix of the estimated distribution. The other is \code{mat.coef}, which is the matrix of the coefficients of Z and Z0}
+#' \item{mle}{the estimated parameters of the the density function}
+#' \item{mle.g}{the estimated coefficients of Z}
+#' \item{value}{the optimized penalized negative log-likelihood value}
+#' \item{S}{the fake information proportion}
+#' \item{cov}{the covariance of the parameters}
+#' \item{bias}{the bias of the parameters}
+#' \item{cov.g}{the covaraince of the estimated density points}
+#' \item{cov.g.gamma}{the covariance between the estimated density points and the coefficient of Z}
+#' \item{loglik}{the objective function being optimized}
+#' \item{statsFunction}{the function computing the relavant statistics}
+#'
+#' @note This is an extension of the G-modeling package
+#' @export
+
+deconvG <- 
+function(tau, X,  
+         offset, 
+         family = c("Poisson", "Negative Binomial"),
+         ignoreZero = F, ## whether ignore zero count
+         zeroInflate = F, ## whether add zero inflation part to the deconvolved distribution
+         Z = NULL,
+         Z0 = NULL,
+         c0 = 1, 
+         NB.size = 100, ## over-dispersion parameter when family is Negative Binomial
+         only.value = F,
+         pDegree = 5, 
+         aStart = 1,
+         bStart = 0,
+         gStart = 0,
+         ...) 
+{
+
+  family <- match.arg(family)
+  require(stats)
+  if (!is.null(Z))
+    Z <- as.matrix(Z)
+  if (!is.null(Z0)) {
+    Z0 <- as.matrix(Z0)
+    Z0.means <- colMeans(Z0)
+    Z0 <- scale(Z0, scale = F)
+  }
+    
+  if (ignoreZero) {
+    if (!is.null(Z))
+      Z <- Z[X!=0, , drop= F]
+    if (!missing(offset))
+      offset <- offset[X!=0]
+    X <- X[X != 0]
+  }
+
+  N <- length(X)
+  if (missing(offset))
+    offset <- rep(0, N)
+
+  y <- rep(1, N)
+  lam <- exp(tau)
+
+  size <- NB.size
+
+
+##################### Construction step ################################################  
+
+  if (is.null(Z)) {  ## then P is a matrix of N * m
+    Q <- scale(splines::ns(lam, pDegree), center = TRUE, 
+               scale = FALSE)
+    Q <- svd(Q)$u
+    Q <- apply(Q, 2, function(w) w/sqrt(sum(w * w))) ## matrix of m * (pDegree + 1)
+    if (ignoreZero) {
+      if (family == "Poisson")
+        P <- sapply(tau, function(theta) 
+                    dpois(x = X, 
+                          lambda = exp(theta + offset))/(1 - exp(-(exp(theta + offset)))))
+      else if (family == "Negative Binomial")
+        P <- sapply(tau, function(theta) 
+                    dnbinom(x = X, 
+                            mu = exp(theta + offset), size = size)/
+                    (1 - dnbinom(0, mu = exp(theta + offset), size = size)))
+    } else {
+      # P is dimension N * m
+      if (family == "Poisson")
+        P <- sapply(tau, function(theta) dpois(x = X, 
+                                               lambda = exp(theta + offset)))
+      else if (family == "Negative Binomial")
+        P <- sapply(tau, function(theta) 
+                    dnbinom(x = X, 
+                            mu = exp(theta + offset), size = size))
+      if (zeroInflate) {
+        #      print("check")
+        P <- cbind(as.numeric(X == 0), P)
+        lam <- c(0, lam)
+        Q <- scale(splines::ns(lam, pDegree), center = TRUE, 
+                   scale = FALSE)
+        Q <- svd(Q)$u
+        Q <- apply(Q, 2, function(w) w/sqrt(sum(w * w))) ## matrix of (m + 1) * (pDegree)
+
+        if (is.null(Z0))
+          Q <- cbind(c(1, rep(0, length(tau))), Q)
+        else
+          Q.fun <- function(Z.vec) {
+            add.on <- rbind(c(1, Z.vec), matrix(0, length(tau), length(Z.vec)+1))
+            return(cbind(add.on, Q))
+          }
+      }
+    }
+  } else { ## when Z is not NULL
+    Q <- scale(splines::ns(lam, pDegree), center = TRUE, 
+               scale = FALSE)
+    Q <- svd(Q)$u
+    Q <- apply(Q, 2, function(w) w/sqrt(sum(w * w))) ## matrix of m * (pDegree + 1)
+    if (ignoreZero) {
+      P.fun <- function(gamma) 
+        sapply(tau, function(theta) { 
+               if (family == "Poisson")
+                 dpois(x = X, 
+                       lambda = exp(theta + offset + Z %*% gamma))/ 
+               (1 - exp(-(exp(theta + offset + Z %*% gamma))))
+               else if (family == "Negative Binomial")
+                 dnbinom(x = X, 
+                         mu = exp(theta + offset + Z%*% gamma), size = size)/
+               (1 - dnbinom(0, mu = exp(theta + offset + Z %*% gamma), size = size))
+               })
+    } else {
+      P.fun <- function(gamma) {
+        P <- sapply(tau, function(theta) {
+                    if (family == "Poisson")
+                      dpois(x = X, 
+                            lambda = exp(theta + offset + Z %*% gamma))
+                    else if (family == "Negative Binomial")
+                      dnbinom(x = X, 
+                              mu = exp(theta + offset + Z %*% gamma), size = size)
+               })
+        if (zeroInflate) 
+          P <- cbind(as.numeric(X == 0), P)
+        return(P)
+      }
+      if (zeroInflate) {
+        lam <- c(0, lam)
+        Q <- scale(splines::ns(lam, pDegree), center = TRUE, 
+                   scale = FALSE)
+        Q <- svd(Q)$u
+        Q <- apply(Q, 2, function(w) w/sqrt(sum(w * w))) ## matrix of (m+1) * (pDegree)
+
+        if (is.null(Z0))
+          Q <- cbind(c(1, rep(0, length(tau))), Q)
+        else
+          Q.fun <- function(Z.vec) {
+            add.on <- rbind(c(1, Z.vec), matrix(0, length(tau), length(Z.vec)+1))
+            return(cbind(add.on, Q))
+          }
+      }
+    }
+  }
+
+######################### Optimization Step ########################################
+  # Initialization
+  p <- ncol(Q)
+  if (!is.null(Z0))
+    p <- p + 1 + ncol(Z0)
+
+    
+  pGiven <- length(aStart)
+  if (pGiven == 1) {
+    aStart <- rep(aStart, p)
+    if (!is.null(Z0)) {
+      if (length(bStart) == 1)
+        bStart <- rep(bStart, ncol(Z0))
+      aStart[2:(1+ ncol(Z0))] <- bStart
+    }
+  }
+  else {
+    if (pGiven != p) 
+      stop(sprintf("Wrong length (%d) for initial parameter, expecting length (%d)", 
+                   pGiven, p))
+  }
+
+  if (!is.null(Z)) {
+    if (length(gStart) == 1)
+      gStart <- rep(gStart, ncol(Z))
+  }
+
+  if (!is.null(Z0)) 
+    g.fun <- function(Z.vec, a) {
+      Qi <- Q.fun(Z.vec)
+      gi <- exp(Qi %*% a)
+      gi <- as.vector(gi /sum(gi))
+      return(gi)
+    }
+
+
+  ## optimization objective function
+  if (is.null(Z)) {
+    loglik <- function(a) {
+      if (is.null(Z0)) {
+        g <- exp(Q %*% a)
+        g <- as.vector(g/sum(g))
+        f <- as.vector(P %*% g)
+        ## objective value
+        value <- -sum(y * log(f)) + c0 * sum(a^2)^0.5
+        Pt <- P/f
+        W <- g * (t(Pt) - 1)
+        ## derivative
+        qw <- t(Q) %*% W
+      } else {
+        m.temp <- 1 + ncol(Z0)
+        g.temp <- Q %*% a[-(1:m.temp)]
+        G.temp1 <- cbind(rep(1, nrow(Z0)), Z0) %*% a[1:m.temp]
+        G <- matrix(rep(g.temp, nrow(Z0)), nrow = length(g.temp))
+        G[1, ] <- G[1, ] + G.temp1
+        G <- exp(G)
+        G <- t(G) / colSums(G)
+        f <- rowSums(P * G)
+        value <- -sum(y * log(f)) + c0 * sum(a^2)^0.5
+        Pt <- P/f
+        W <- t(G * (Pt - 1)) # dimension is as G: N * (m + 1)
+
+        qw.temp2 <- t(Q) %*% W
+        qw.temp1 <- cbind(rep(1, nrow(Z0)), Z0) * W[1, ]
+        qw <- rbind(t(qw.temp1), qw.temp2)
+
+      }
+      aa <- sqrt(sum(a^2))
+      sDot <- c0 * a/aa
+
+      attr(value, "gradient") <- -(qw %*% y) + sDot
+      value
+    }
+
+    result <- stats::nlm(f = loglik, p = aStart, gradtol = 1e-5,
+                         ...)
+    mle <- result$estimate
+    value <- loglik(mle)
+    mle.a <- result$estimate
+    mle.g <- NULL
+
+  } else {
+    loglik <- function(a.gamma) {
+      a <- a.gamma[1:p]
+      gamma <- a.gamma[(p+1):length(a.gamma)]
+      P <- P.fun(gamma)
+      if (is.null(Z0)) {
+        g <- exp(Q %*% a)
+        g <- as.vector(g/sum(g))
+        f <- as.vector(P %*% g)
+        ## objective value
+        value <- -sum(y * log(f)) + c0 * sum(a^2)^0.5
+        Pt <- P/f
+        W <- g * (t(Pt) - 1)
+        ## derivative
+        qw <- t(Q) %*% W
+      } else {
+        m.temp <- 1 + ncol(Z0)
+        g.temp <- Q %*% a[-(1:m.temp)]
+        G.temp1 <- cbind(rep(1, nrow(Z0)), Z0) %*% a[1:m.temp]
+        G <- matrix(rep(g.temp, nrow(Z0)), nrow = length(g.temp))
+        G[1, ] <- G[1, ] + G.temp1
+        G <- exp(G)
+        G <- t(G) / colSums(G)
+
+        f <- rowSums(P * G)
+        value <- -sum(y * log(f)) + c0 * sum(a^2)^0.5
+        Pt <- P/f
+        W <- t(G * (Pt - 1)) # dimension is as G: N * (m + 1)
+
+        qw.temp2 <- t(Q) %*% W
+        qw.temp1 <- cbind(rep(1, nrow(Z0)), Z0) * W[1, ]
+        qw <- rbind(t(qw.temp1), qw.temp2)
+      }
+
+      aa <- sqrt(sum(a^2))
+      sDot <- c0 * a/aa
+      dot.l.a <- -rowSums(qw) + sDot # length p
+
+      # derivative for gamma
+      eta <- rep(1, N) %*% t(tau) + as.vector(Z %*% gamma) + offset
+      if (family == "Poisson" || family == "Negative Binomial")
+        mu <- exp(eta)
+      if (zeroInflate){
+        P.tilde <- P[, -1]
+        if (is.null(Z0))
+          g.tilde <- g[-1]
+        else
+          G.tilde <- G[, -1]
+      } else {
+        P.tilde <- P
+        g.tilde <- g
+      }
+      if (family == "Poisson")
+        A <- (X - mu) * P.tilde # dimension is N * m
+      if (family == "Negative Binomial")
+        A <- (X - mu) / (1 + mu/size) * P.tilde
+      if (is.null(Z0))
+        agf <- A %*% g.tilde / f
+      else
+        agf <- rowSums(A * G.tilde)/f
+      dot.l.gamma <- -t(Z) %*% agf
+
+      attr(value, "gradient") <- c(dot.l.a, dot.l.gamma)
+
+      value
+    }
+    result <- stats::nlm(f = loglik, p = c(aStart, gStart), 
+                         gradtol = 1e-5, 
+                        ...)
+
+    mle <- result$estimate
+    mle.a <- mle[1:p]
+    mle.g <- mle[-(1:p)]
+    value <- loglik(mle)
+  }
+
+
+################################# inference functions #####################################
+
+  ## used when Z is NULL
+  statsFunction <- function(a) {
+    if (is.null(Z0)) {
+      g <- exp(Q %*% a)
+      g <- as.vector(g/sum(g))
+      f <- as.vector(P %*% g)
+      ## objective value
+
+      value <- -sum(y * log(f)) + c0 * sum(a^2)^0.5
+      Pt <- P/f
+      W <- g * (t(Pt) - 1)
+      ## derivative
+      qw <- t(Q) %*% W
+      qG <- matrix(rep(t(Q) %*% g, ncol(W)), ncol = ncol(W))
+      Wplus <- rowSums(W)
+      qWq <- t(Q) %*% (Wplus * Q)
+    } else {
+      m.temp <- 1 + ncol(Z0)
+      g.temp <- Q %*% a[-(1:m.temp)]
+      G.temp1 <- cbind(rep(1, nrow(Z0)), Z0) %*% a[1:m.temp]
+      G <- matrix(rep(g.temp, nrow(Z0)), nrow = length(g.temp))
+      G[1, ] <- G[1, ] + G.temp1
+      G <- exp(G)
+      G <- t(G) / colSums(G)
+
+      f <- rowSums(P * G)
+
+      value <- sum(y * log(f)) + c0 * sum(a^2)^0.5
+      Pt <- P/f
+      W <- t(G * (Pt - 1)) # dimension is as G: N * (m + 1)
+      qw <- sapply(1:nrow(Z0), function(i) {
+                   Qi <- Q.fun(Z0[i, ])
+                   return(t(Qi) %*% W[, i])
+                         })
+      g <- colMeans(G)
+      qG <- sapply(1:nrow(Z0), function(i) {
+                         Qi <- Q.fun(Z0[i, ])
+                         return(t(Qi) %*% G[i, ])
+                         })
+      qWq <- matrix(0, nrow(qw), nrow(qw))
+      for (i in 1:nrow(Z0)) {
+        Qi <- Q.fun(Z0[i, ])
+        qWq <- qWq + t(Qi) %*% (W[, i] * Qi)
+      }
+    }
+
+
+    yHat <- if (length(y) == 1 && y == 1 || length(y) == length(X)) 
+      y
+    else sum(y) * f
+    # estimate W
+    I1 <- qw %*% (yHat * t(qw))
+
+    aa <- sqrt(sum(a^2))
+    sDot <- c0 * a/aa
+    sDotDot <- (c0/aa) * (diag(length(a)) - outer(a, a)/aa^2)
+    R <- sum(diag(sDotDot))/sum(diag(I1))
+    temp <- qw %*% t(qG)
+    I2 <- solve(I1 + temp + t(temp) - qWq + sDotDot)
+
+    bias <- as.vector(-I2 %*% sDot)
+    Cov <- I2 %*% (I1 %*% t(I2))
+
+    if (is.null(Z0))
+      mat.coef <- NULL
+    else {
+      idx <- 2:(1+ncol(Z0))
+      mat.coef <- cbind(mle[idx], 
+                         bias[idx],
+                         sqrt(diag(Cov[idx, idx, drop = F])))
+      Q0 <- cbind(rbind(c(1, -Z0.means), matrix(0, nrow(Q) - 1, 1 + ncol(Z0))), Q)
+      temp <- as.vector(exp(Q0 %*% a))
+      beta0 <- log(temp[1]/sum(temp[-1]))
+      temp.dev <- temp * Q0
+      beta0.dev <- c(1/temp[1], -rep(1/sum(temp[-1]), length(temp)-1))
+      beta0.dev <- t(temp.dev) %*% beta0.dev
+      beta0 <- c(beta0, t(beta0.dev) %*% bias, 
+                 sqrt(t(beta0.dev) %*% Cov %*% beta0.dev))
+      mat.coef <- rbind(beta0, mat.coef)
+      rownames(mat.coef) <- paste("Z0 effect: beta", 0:ncol(Z0), sep = "")
+      colnames(mat.coef) <- c("est", "bias", "sd")
+    }
+
+
+    if (is.null(Z0)) {
+      Dq <- (diag(g) - outer(g, g)) %*% Q
+    } else {
+      Dq <- matrix(0, ncol(G), p)
+      sapply(1:nrow(G), function(i) {
+             Qi <- Q.fun(Z0[i, ])
+             temp <- (diag(G[i, ]) - outer(G[i, ], G[i, ])) %*% Qi
+             Dq <<- Dq + temp
+                         })
+      Dq <- Dq/nrow(G)
+    }
+    bias.g <- Dq %*% bias
+    Cov.g <- Dq %*% Cov %*% t(Dq)
+    se.g <- diag(Cov.g)^0.5
+    D <- diag(length(lam))
+    D[lower.tri(D)] <- 1
+    accuG <- cumsum(g)
+    Cov.G <- D %*% (Cov.g %*% t(D))
+    se.G <- diag(Cov.G)^0.5
+    mat <- cbind(lam, g, se.g, accuG, se.G, bias.g)
+    colnames(mat) = c("theta", "g", "SE.g", "G", "SE.G", 
+                      "Bias.g")
+    list(S = R, cov = Cov, bias = bias, cov.g = Cov.g, cov.g.gamma= NULL,
+         mat = list(mat.dist = mat, mat.coef = mat.coef))
+  }
+
+  
+
+  ## Used when Z is not NULL
+  statsFunctionZ <- function(a, gamma) {
+    P <- P.fun(gamma)
+    if (is.null(Z0)) {
+      g <- exp(Q %*% a)
+      g <- as.vector(g/sum(g))
+      f <- as.vector(P %*% g)
+      ## objective value
+      value <- -sum(y * log(f)) + c0 * sum(a^2)^0.5
+      Pt <- P/f
+      W <- g * (t(Pt) - 1)
+      ## derivative
+      qw <- t(Q) %*% W
+      qG <- matrix(rep(t(Q) %*% g, ncol(W)), ncol = ncol(W))
+      Wplus <- rowSums(W)
+      qWq <- t(Q) %*% (Wplus * Q)
+    } else {
+      m.temp <- 1 + ncol(Z0)
+      g.temp <- Q %*% a[-(1:m.temp)]
+      G.temp1 <- cbind(rep(1, nrow(Z0)), Z0) %*% a[1:m.temp]
+      G <- matrix(rep(g.temp, nrow(Z0)), nrow = length(g.temp))
+      G[1, ] <- G[1, ] + G.temp1
+      G <- exp(G)
+      G <- t(G) / colSums(G)
+
+      f <- rowSums(P * G)
+      value <- sum(y * log(f)) + c0 * sum(a^2)^0.5
+      Pt <- P/f
+      W <- t(G * (Pt - 1)) # dimension is as G: N * (m + 1)
+      qw <- sapply(1:nrow(Z0), function(i) {
+                   Qi <- Q.fun(Z0[i, ])
+                   return(t(Qi) %*% W[, i])
+                      })
+      g <- colMeans(G)
+      qG <- sapply(1:nrow(Z0), function(i) {
+                         Qi <- Q.fun(Z0[i, ])
+                         return(t(Qi) %*% G[i, ])
+                         })
+      qWq <- matrix(0, nrow(qw), nrow(qw))
+      for (i in 1:nrow(Z0)) {
+        Qi <- Q.fun(Z0[i, ])
+        qWq <- qWq + t(Qi) %*% (W[, i] * Qi)
+      }
+    }
+
+    ## compute Igg
+    eta <- rep(1, N) %*% t(tau) + as.vector(Z %*% gamma) + offset
+    if (family == "Poisson" || family == "Negative Binomial") {
+      mu <- exp(eta)
+    } 
+    if (zeroInflate){
+      P.tilde <- P[, -1]
+      if (is.null(Z0)) {
+        g.tilde <- g[-1]
+        Q.tilde <- Q[-1, ]
+      }
+      else {
+        G.tilde <- G[, -1]
+        g.tilde <- G.tilde[1, ]
+        Q.tilde <- cbind(matrix(0, nrow(Q)-1, 1 + ncol(Z0)), Q[-1, ])
+      }
+    } else {
+      P.tilde <- P
+      g.tilde <- g
+      Q.tilde <- Q
+    }
+    if (family == "Poisson") {
+      A <- (X - mu) * P.tilde # dimension is N * m
+      B <- ((X - mu)^2 - mu) * P.tilde
+    }
+    if (family == "Negative Binomial") {
+      A <- (X - mu) / (1 + mu/size) * P.tilde
+      B <- ((X - mu)^2 / (1 + mu/size)^2 - mu/(1 + mu/size) - 
+            (X - mu) * mu / (1 + mu/size)^2/size) * P.tilde
+    }
+
+
+    if (is.null(Z0)) {
+      agf <- as.vector(A %*% g.tilde / f)
+      bgf <- as.vector(B %*% g.tilde / f)
+    }
+    else {
+      agf <- rowSums(A * G.tilde)/f
+      bgf <- rowSums(B * G.tilde)/f
+    }
+
+
+    dot.li.gamma <- t(agf * Z)
+
+    deriv.comb <- rbind(qw, dot.li.gamma)
+    aa <- sqrt(sum(a^2))
+    sDot <- c0 * a/aa
+    sDotDot <- (c0/aa) * (diag(length(a)) - outer(a, a)/aa^2)
+
+    I1full <- deriv.comb %*% t(deriv.comb)
+    Iaa <- I1full[1:length(a), 1:length(a)]
+    Igg <- t(Z) %*% (as.vector(agf^2 - bgf) * Z)
+    if (is.null(Z0))
+      Iag <- t(Q.tilde) %*% (g.tilde * (t(P.tilde/f) %*% (agf * Z) - t(A/f) %*% Z))
+    else
+      Iag <- t(Q.tilde) %*% (t(G.tilde * P.tilde/f) %*% (agf * Z) - t(G.tilde * A/f) %*% Z)
+    I2full <- rbind(cbind(Iaa, Iag), cbind(t(Iag), Igg))
+    I2full[1:length(a), 1:length(a)] <- I2full[1:length(a), 1:length(a)] + 
+                                        sDotDot
+
+    Ifull.inv <- solve(I2full)
+
+      ## ratio of artificial to genuine information
+    R <- sum(diag(sDotDot))/sum(diag(I1full[1:length(a), 1:length(a)])) 
+
+
+    bias <- as.vector(-Ifull.inv %*% c(sDot, rep(0, ncol(Z))))
+    Cov <- Ifull.inv %*% (I1full %*% t(Ifull.inv))
+    if (is.null(Z0)) {
+      mat.coef <- cbind(mle.g, bias[-(1:length(a))],
+                        sqrt(diag(Cov[-(1:length(a)), -(1:length(a)), drop = F])))
+      rownames(mat.coef) <- c(paste("Z:gamma", 1:ncol(Z), sep = ""))
+      colnames(mat.coef) <- c("est", "bias", "sd")
+    }
+    else {
+      idx1 <- 2:(1+ncol(Z0))
+      idx <- c((length(a)+1):length(bias), idx1)
+      mat.coef <- cbind(mle[idx], 
+                         bias[idx],
+                         sqrt(diag(Cov[idx, idx, drop = F])))
+
+      Q0 <- cbind(rbind(c(1, -Z0.means), matrix(0, nrow(Q) - 1, 1 + ncol(Z0))), Q)
+      temp <- as.vector(exp(Q0 %*% a))
+      beta0 <- log(temp[1]/sum(temp[-1]))
+      temp.dev <- temp * Q0
+      beta0.dev <- c(1/temp[1], -rep(1/sum(temp[-1]), length(temp)-1))
+      beta0.dev <- t(temp.dev) %*% beta0.dev
+      beta0 <- c(beta0, t(beta0.dev) %*% bias[1:length(a)], 
+                 sqrt(t(beta0.dev) %*% Cov[1:length(a), 1:length(a), drop = F] %*% beta0.dev))
+      mat.coef <- rbind(mat.coef[1:ncol(Z), ],
+                        beta0,
+                        mat.coef[-(1:ncol(Z)), ])
+
+      rownames(mat.coef) <- c(paste("Z effect: gamma", 1:ncol(Z), sep = ""),
+                              paste("Z0 effect: beta", 0:ncol(Z0), sep = ""))
+      colnames(mat.coef) <- c("est", "bias", "sd")
+    }
+
+    if (is.null(Z0)) {
+      Dq <- (diag(g) - outer(g, g)) %*% Q
+    } else {
+      Dq <- matrix(0, ncol(G), p)
+      sapply(1:nrow(G), function(i) {
+             Qi <- Q.fun(Z0[i, ])
+             temp <- (diag(G[i, ]) - outer(G[i, ], G[i, ])) %*% Qi
+             Dq <<- Dq + temp
+                      })
+      Dq <- Dq/nrow(G)
+    }
+
+    bias.g <- Dq %*% bias[1:p]
+    Cov.g <- Dq %*% Cov[1:p, 1:p] %*% t(Dq)
+    Cov.g.gamma <- Dq %*% Cov[1:p, -(1:p), drop = F]
+    se.g <- diag(Cov.g)^0.5
+    D <- diag(length(lam))
+    D[lower.tri(D)] <- 1
+    accuG <- cumsum(g)
+    Cov.G <- D %*% (Cov.g %*% t(D))
+    se.G <- diag(Cov.G)^0.5
+    mat <- cbind(lam, g, se.g, accuG, se.G, bias.g)
+    colnames(mat) = c("theta", "g", "SE.g", "G", "SE.G", 
+                      "Bias.g")
+    list(S = R, cov = Cov, bias = bias, cov.g = Cov.g, cov.g.gamma = Cov.g.gamma,
+         mat = list(mat.dist = mat, mat.coef = mat.coef))
+  }
+
+
+  if (only.value)
+    return(value)
+
+  if (is.null(Z))
+    stats <- statsFunction(mle.a)
+  else
+    stats <- statsFunctionZ(mle.a, mle.g)
+  list(stats = stats$mat, 
+       mle = mle.a, 
+       mle.g = mle.g,
+       value = value,
+       S = stats$S, 
+       cov = stats$cov, 
+       bias = stats$bias,
+       cov.g = stats$cov.g, 
+       cov.g.gamma = stats$cov.g.gamma,
+       loglik = loglik, 
+       statsFunction = statsFunction)
+}
+
+
+
diff --git a/descend.Rproj b/descend.Rproj
new file mode 100644
index 0000000..d848a9f
--- /dev/null
+++ b/descend.Rproj
@@ -0,0 +1,16 @@
+Version: 1.0
+
+RestoreWorkspace: No
+SaveWorkspace: No
+AlwaysSaveHistory: Default
+
+EnableCodeIndexing: Yes
+Encoding: UTF-8
+
+AutoAppendNewline: Yes
+StripTrailingWhitespace: Yes
+
+BuildType: Package
+PackageUseDevtools: Yes
+PackageInstallArgs: --no-multiarch --with-keep.source
+PackageRoxygenize: rd,collate,namespace
diff --git a/man/DESCEND.control.Rd b/man/DESCEND.control.Rd
new file mode 100644
index 0000000..87017ea
--- /dev/null
+++ b/man/DESCEND.control.Rd
@@ -0,0 +1,16 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/deconvSingle.R
+\name{DESCEND.control}
+\alias{DESCEND.control}
+\title{The control parameters of the function \code{\link{descend}} and \code{\link{deconvSingle}}}
+\usage{
+DESCEND.control(n.points = 50, nStart = 2, nStart.lrt = 2,
+  max.quantile = 0.98, max.sparse = 0.95, rel.info.range = c(5e-04, 0.01),
+  rel.info.guide = 0.005, c0.start = 1, aStart = 1, bStart = 0,
+  gStart = 0, start.sd = 0.2, penalty.Z0 = T, pDegree = 5,
+  max.c0.iter = 5, c0.min = 1e-05, LRT.Z.select = T, LRT.Z0.select = T,
+  LRT.Z.values = 0, zeroInflate = T, dense.add.0 = T, only.integer = F)
+}
+\description{
+The control parameters of the function \code{\link{descend}} and \code{\link{deconvSingle}}
+}
diff --git a/man/deconvG.Rd b/man/deconvG.Rd
new file mode 100644
index 0000000..a31e355
--- /dev/null
+++ b/man/deconvG.Rd
@@ -0,0 +1,62 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/g_model.R
+\name{deconvG}
+\alias{deconvG}
+\title{Base function of G-modeling deconvolution}
+\usage{
+deconvG(tau, X, offset, family = c("Poisson", "Negative Binomial"),
+  ignoreZero = F, zeroInflate = F, Z = NULL, Z0 = NULL, c0 = 1,
+  NB.size = 100, only.value = F, pDegree = 5, aStart = 1, bStart = 0,
+  gStart = 0, ...)
+}
+\arguments{
+\item{tau}{log of the discrete points of the deconvolved distribution}
+
+\item{X}{a vector of observed counts}
+
+\item{family}{family of the noise distribution, support either "Poisson" or "Negative Binomial" with known tuning parameter}
+
+\item{ignoreZero}{whether ignore the zero count. If true, then use truncated Poisson / Negative Binomial distribution. Default is False}
+
+\item{zeroInflate}{whether add zero inflation part to the deconvolved distribution to reflect transcriptional bursting. Default is True.}
+
+\item{Z}{covariates for active intensity. Default is NULL.}
+
+\item{Z0}{covariates for active fraction. Used only when zeroInflate is True. Default is NULL.}
+
+\item{c0}{the tuning parameter on the L2 penalty term. Default is 1. c0 will be selected automatically in \code{deconvSingle}}
+
+\item{NB.size}{over-dispersion parameter when the family is Negative Binomial: mu = mu + mu^2/size}
+
+\item{only.value}{whether not to compute the estimation statistics but only the value of the optimized lieklihood. Used for likelihood ratio test.}
+
+\item{pDegree}{the degree of the Spline matrix. Default is 5.}
+
+\item{aStart, bStart, gStart}{initial values of the density parameters, the coefficients of Z0 and coefficients of Z}
+
+\item{...}{extra parameters for the \code{\link[stats]{nlm}} function}
+}
+\value{
+a list with elements
+\item{stats}{a list of two elements. One is the \code{mat.dist}, which is the matrix of the estimated distribution. The other is \code{mat.coef}, which is the matrix of the coefficients of Z and Z0}
+\item{mle}{the estimated parameters of the the density function}
+\item{mle.g}{the estimated coefficients of Z}
+\item{value}{the optimized penalized negative log-likelihood value}
+\item{S}{the fake information proportion}
+\item{cov}{the covariance of the parameters}
+\item{bias}{the bias of the parameters}
+\item{cov.g}{the covaraince of the estimated density points}
+\item{cov.g.gamma}{the covariance between the estimated density points and the coefficient of Z}
+\item{loglik}{the objective function being optimized}
+\item{statsFunction}{the function computing the relavant statistics}
+}
+\description{
+Base function used by \code{\link{deconvSingle}} to deconvolve the underlying distribution. 
+We assume X \sim F(T) where F is the noise distribution. We assume that 
+\deqn{log(T) = offset + \gammaZ + \epsilon}
+\deqn{P(T = 0) = \beta_0 + \beta_1Z0}
+The goal is the recover the distribution of exp(log(T) - offset - \gammaZ), which has density g and is discretized at exp(tau) (add 0 when zero inflation is allowed).
+}
+\note{
+This is an extension of the G-modeling package
+}
diff --git a/man/deconvSingle.Rd b/man/deconvSingle.Rd
new file mode 100644
index 0000000..f7e3c67
--- /dev/null
+++ b/man/deconvSingle.Rd
@@ -0,0 +1,37 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/deconvSingle.R
+\name{deconvSingle}
+\alias{deconvSingle}
+\title{Deconvolve the true gene expression distribution  of a single gene}
+\usage{
+deconvSingle(y, scaling.consts = NULL, Z = NULL, Z0 = NULL,
+  plot.density = T, do.LRT.test = T, family = c("Poisson",
+  "Negative Binomial"), NB.size = 100, verbose = T, control = list())
+}
+\arguments{
+\item{y}{a vector of observed counts of a single gene}
+
+\item{scaling.consts}{a vector of cell specific scaling constants, either the cell efficiency or the library size}
+
+\item{Z}{covariates for active intensity. Default is NULL.}
+
+\item{Z0}{covariates for active fraction. Used only when zeroInflate is True. Default is NULL.}
+
+\item{do.LRT.test}{whether do LRT test on the coefficients and active fraction or not. Default is True}
+
+\item{family}{family of the noise distribution, support either "Poisson" or "Negative Binomial" with known tuning parameter}
+
+\item{NB.size}{over-dispersion parameter when the family is Negative Binomial: mu = mu + mu^2/size}
+
+\item{verbose}{verbose the estimation and testing procedures or not. Default is True.}
+
+\item{control}{settings see see {\code{\link{deconv.control}}}}
+
+\item{whether}{plot the density curve of the deconvolved the distribution or not. The zero inflation part has been smoothed into the density curve for visualization. Default is True.}
+}
+\value{
+a DESCEND object
+}
+\description{
+The deconvolution is computed by using the function \code{deconvG}. This function can automatically discretize the underlying distribution and find the proper tuning parameter of the penalty term in G-modelling. Besides, it computs the estimates and standard deviations of five distribution based statistics (active fraction, active intensity, mean, CV and gini coefficient), as well as the estimated coefficients of the covariates on active intensity (Z) and active fraction (Z0).
+}
diff --git a/man/descend.Rd b/man/descend.Rd
new file mode 100644
index 0000000..95148f5
--- /dev/null
+++ b/man/descend.Rd
@@ -0,0 +1,39 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/descend.R
+\name{descend}
+\alias{descend}
+\title{Default is to set max.quantile = max.sparse = 0.98. For those that are too sparse, set max.quantile = max.sparse = 0.99 and try again.}
+\usage{
+descend(count.matrix, scaling.consts = NULL, Z = NULL, Z0 = NULL,
+  n.cores = 1, cl = NULL, do.LRT.test = F, family = c("Poisson",
+  "Negative Binomial"), NB.size = 100, verbose = T, control = list())
+}
+\arguments{
+\item{count.matrix}{the observed UMI count matrix. Each row is a gene and each column is a cell}
+
+\item{scaling.consts}{a vector of cell specific scaling constants, either the cell efficiency or the library size}
+
+\item{Z}{covariates for active intensity. Default is NULL.}
+
+\item{Z0}{covariates for active fraction. Used only when zeroInflate is True. Default is NULL.}
+
+\item{n.cores}{the number of cores used for parallel computing. Default is 1. Used only when parallel computing is done in a single machine. For using multi-machine cores, need to assign \code{cl} explicitly}
+
+\item{cl}{an object of class "cluster". See more details in \code{\link[parallel]{makeCluster}}}
+
+\item{do.LRT.test}{whether do LRT test on the coefficients and active fraction or not. Default is True}
+
+\item{family}{family of the noise distribution, support either "Poisson" or "Negative Binomial" with known tuning parameter}
+
+\item{NB.size}{over-dispersion parameter when the family is Negative Binomial: mu = mu + mu^2/size}
+
+\item{verbose}{verbose the estimation and testing procedures or not. Default is True.}
+
+\item{control}{settings see see {\code{\link{deconv.control}}}}
+}
+\value{
+a list of DESCEND objects
+}
+\description{
+Default is to set max.quantile = max.sparse = 0.98. For those that are too sparse, set max.quantile = max.sparse = 0.99 and try again.
+}
diff --git a/man/getEstimates.Rd b/man/getEstimates.Rd
new file mode 100644
index 0000000..4086759
--- /dev/null
+++ b/man/getEstimates.Rd
@@ -0,0 +1,14 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/descend.R
+\name{getEstimates}
+\alias{getEstimates}
+\title{Grep the value and standard deviation of the estimated parameters calculated from a list of descend objects}
+\usage{
+getEstimates(descend.list)
+}
+\arguments{
+\item{descend.list}{a list of descend objects}
+}
+\description{
+Grep the value and standard deviation of the estimated parameters calculated from a list of descend objects
+}