@@ -29,7 +29,8 @@ Imports:

     cluster,

     digest,

     gridExtra,

-    methods

+    methods,

+    lintr

 Suggests:

     testthat,

     knitr,

@@ -26,7 +26,6 @@ export(calculate_loglik_from_variables)

 export(calculate_loglik_from_variables.celda_C)

 export(calculate_loglik_from_variables.celda_CG)

 export(calculate_loglik_from_variables.celda_G)

-export(calculate_marginal_likelihood)

 export(calculate_perplexity)

 export(celda)

 export(celda_C)

@@ -405,7 +405,7 @@ celda_heatmap.celda_C = function(celda.mod, counts, ...) {

 #' visualize_model_performance for celda Cell clustering function

 #' @param celda.list A celda_list object returned from celda()

-#' @param method One of "perplexity", "harmonic", or "loglik"

+#' @param method One of "perplexity", "loglik"

 #' @param title Title for the plot

 #' @param log Currently not working for celda_C objects

 #' @import Rmpfr

@@ -423,7 +423,7 @@ visualize_model_performance.celda_C = function(celda.list, method="perplexity",

   # These methods return Rmpfr numbers that are extremely small and can't be 

   # plotted, so log 'em first

-  if (method %in% c("perplexity", "harmonic")) {

+  if (method %in% c("perplexity")) {

     performance.metric = lapply(performance.metric, log)

     performance.metric = methods::new("mpfr", unlist(performance.metric))

     performance.metric = as.numeric(performance.metric)

@@ -669,9 +669,9 @@ celda_heatmap.celda_CG = function(celda.mod, counts, ...) {

 #' L (cell clusters K), plot the performance of each number of cell

 #' clusters K (gene clusters L).

 #' @param celda.list A list of celda_CG objects returned from celda function

-#' @param method One of "perplexity", "harmonic", or "loglik"

+#' @param method One of "perplexity" or "loglik"

 #' @param title Title for the visualize_model_performance

-#' @param log Set log to TRUE to visualize the log(perplexity) of Celda_CG objects. Does not work for "harmonic" metric

+#' @param log Set log to TRUE to visualize the log(perplexity) of Celda_CG objects.

 #' @import Rmpfr

 #' @export

 visualize_model_performance.celda_CG = function(celda.list, method="perplexity",

@@ -696,11 +696,7 @@ visualize_model_performance.celda_CG = function(celda.list, method="perplexity",

       performance.metric = new("mpfr", unlist(performance.metric))

     }

     y.lab = paste0("Log(",method,")")

-  } else if (method == "harmonic") {

-    performance.metric = lapply(performance.metric, log)

-    performance.metric = new("mpfr", unlist(performance.metric))

-    y.lab = paste0("Log(",method,")")

-  }

+  } 

   performance.metric = as.numeric(performance.metric)

@@ -761,7 +757,7 @@ visualize_model_performance.celda_CG = function(celda.list, method="perplexity",

 #' run for each combination of K/L (cell/gene).

 #' 

 #' @param celda.list A list of celda_CG objects returned from celda function

-#' @param method One of "perplexity", "harmonic", or "loglik", passed through to calculate_performance_metric()

+#' @param method One of "perplexity" or "loglik", passed through to calculate_performance_metric()

 #' @param title The plot title

 #' @import Rmpfr 

 #' @export

@@ -783,7 +779,7 @@ render_interactive_kl_plot = function(celda.list,  method="perplexity",

   # The performance metric methods return Rmpfr numbers that are extremely small and can't be 

   # plotted via ggplot2, so log 'em first. 

   # TODO: celda_list getter that calculates these metrics.

-  if (method %in% c("perplexity", "harmonic")) {

+  if (method %in% c("perplexity")) {

     performance.metric = lapply(performance.metric, log)

     performance.metric = new("mpfr", unlist(performance.metric))

     performance.metric = as.numeric(performance.metric)

@@ -803,7 +799,7 @@ render_interactive_kl_plot = function(celda.list,  method="perplexity",

   # TODO: Return plot or nah?

   method.label = method 

-  if (method %in% c("perplexity, harmonic")) {

+  if (method %in% c("perplexity")) {

      method.label = paste("log(", method, ")", sep="")

   }

   k.l.plot = ggplot2::ggplot(figure.df, aes(x=key, y=metric, label=key)) +

@@ -824,7 +820,7 @@ validate_kl_plot_parameters = function(celda.list, method) {

     stop("celda.list argument must be of class 'celda_list'")

  } else if (celda.list$content.type != "celda_CG") {

     stop("celda.list must be a 'celda.list' of 'celda_CG' objects")

- } else if (!(method %in% c("perplexity","harmonic","loglik"))) {

-    stop("Invalid method, 'method' has to be either 'perplexity', 'harmonic', or 'loglik'")

+ } else if (!(method %in% c("perplexity", "loglik"))) {

+    stop("Invalid method, 'method' has to be either 'perplexity' or 'loglik'")

  } 

 }

@@ -481,7 +481,7 @@ celda_heatmap.celda_G = function(celda.mod, counts, ...) {

 # TODO DRYer implementation in concert with celda_C

 #' visualize_model_performance for the celda Gene function

 #' @param celda.list A celda_list object returned from celda()

-#' @param method One of "perplexity", "harmonic", or "loglik"

+#' @param method One of "perplexity" or "loglik"

 #' @param title Title for the plot

 #' @param log Currently not working for celda.G objects

 #' @import Rmpfr

@@ -499,7 +499,7 @@ visualize_model_performance.celda_G = function(celda.list, method="perplexity",

   # These methods return Rmpfr numbers that are extremely small and can't be 

   # plotted, so log 'em first

-  if (method %in% c("perplexity", "harmonic")) {

+  if (method %in% c("perplexity")) {

     performance.metric = lapply(performance.metric, log)

     performance.metric = methods::new("mpfr", unlist(performance.metric))

     performance.metric = as.numeric(performance.metric)

@@ -10,7 +10,7 @@

 #' @param K The K parameter for the desired model in the results list

 #' @param L The L parameter for the desired model in the results list

 #' @param chain The desired chain for the specified model

-#' @param best Method for choosing best chain automatically. Options are c("perplexity", "harmonic", "loglik"). See documentation for chooseBestModel for details. Overrides chain parameter if provided.

+#' @param best Method for choosing best chain automatically. Options are c("perplexity", "loglik"). See documentation for chooseBestModel for details. Overrides chain parameter if provided.

 #' @return A celda model object matching the provided parameters (of class "celda_C", "celda_G", "celda_CG" accordingly), or NA if one is not found.

 #' @export

 getModel = function(celda.list, K=NULL, L=NULL, chain=1, best=NULL) {

@@ -90,7 +90,7 @@ validate_get_model_params = function(celda.list, K, L, chain, best) {

 #' Determine the best chain among a set of celda_* objects with

 #' otherwise uniform K/L choices.

 #' @param celda.mods A list of celda class objects (celda_C, celda_CG, celda_G)

-#' @param method How to choose the best chain. Choices are c("perplexity", "harmonic", "loglik"). Defaults to perplexity. "perplexity" calculates each chain's perplexity as the inverse of the geometric mean, per the original LDA description. "harmonic" calculates each chain's marginal likelihood as the harmonic mean of each iteration of Gibbs sampling's log likelihoods. "loglik" chooses the chain which reached the maximal log likelihood during Gibbs sampling.

+#' @param method How to choose the best chain. Choices are c("perplexity", "loglik"). Defaults to perplexity. "perplexity" calculates each chain's perplexity as the inverse of the geometric mean, per the original LDA description. "loglik" chooses the chain which reached the maximal log likelihood during Gibbs sampling.

 chooseBestChain = function(celda.mods, method="perplexity") {

   # We want to get the *most negative* perplexity, as opposed to the *least* negative

   # for the other metrics...

@@ -101,10 +101,6 @@ chooseBestChain = function(celda.mods, method="perplexity") {

     return(celda.mods[[best]])

   } 

-  else if (method == "harmonic"){

-    metrics = lapply(celda.mods, function(mod) { calculate_perplexity(mod$completeLogLik) })

-    metrics = new("mpfr", unlist(metrics))

-   }

   else if (method == "loglik"){

     metrics = lapply(celda.mods, function(mod) { max(mod$completeLogLik) })

     metrics = unlist(metrics)

@@ -131,4 +127,4 @@ search_res_list = function(celda_list, K=NULL, L=NULL) {

     stop("K/L parameter(s) requested did not appear for any model in the celda_list. Did you modify the run.params?")

   }

   return(requested.chain)

-}

\ No newline at end of file

+}

@@ -1,19 +1,3 @@

-#' Calculate the marginal likelihood from a single celda model

-#' 

-#' Marginal likelihood is estimated as the harmonic mean of the 

-#' (non-log) likelihood over all iterations of Gibbs sampling.

-#' 

-#' @param completeLogLik The complete Gibbs sampling history of log-likelihoods for a single celda chain

-#' @return The estimated marginal likelihood as an mpfr number

-#' @export

-calculate_marginal_likelihood = function(completeLogLik) {

-  mpfr_log_lik = Rmpfr::mpfr(completeLogLik, 512)

-  complete_likelihood = exp(mpfr_log_lik)

-  marginal_likelihood = (Rmpfr::mean((1/complete_likelihood)))^-1

-  return(marginal_likelihood)

-}

-

-

 #' Calculate the perplexity from a single celda chain

 #' 

 #' Perplexity is defined as the inverse of the geometric mean of the log-likelihoods over all 

@@ -34,10 +18,8 @@ calculate_perplexity = function(completeLogLik, log = FALSE) {

 # Convenience function to calculate performance metrics by specifying a method. 

 calculate_performance_metric = function(log.likelihoods, method="perplexity", log = FALSE) {

-    if (method == "perplexity") {

+  if (method == "perplexity") {

     metric = calculate_perplexity(log.likelihoods, log)

-  } else if (method == "harmonic") {

-    metric = calculate_marginal_likelihood(log.likelihoods)

   } else if (method == "loglik") {

      metric = max(log.likelihoods)

   } else stop("Invalid method specified")

deleted file mode 100644

@@ -1,18 +0,0 @@

-% Generated by roxygen2: do not edit by hand

-% Please edit documentation in R/model_performance.R

-\name{calculate_marginal_likelihood}

-\alias{calculate_marginal_likelihood}

-\title{Calculate the marginal likelihood from a single celda model}

-\usage{

-calculate_marginal_likelihood(completeLogLik)

-}

-\arguments{

-\item{completeLogLik}{The complete Gibbs sampling history of log-likelihoods for a single celda chain}

-}

-\value{

-The estimated marginal likelihood as an mpfr number

-}

-\description{

-Marginal likelihood is estimated as the harmonic mean of the 

-(non-log) likelihood over all iterations of Gibbs sampling.

-}

@@ -10,7 +10,7 @@ chooseBestChain(celda.mods, method = "perplexity")

 \arguments{

 \item{celda.mods}{A list of celda class objects (celda_C, celda_CG, celda_G)}

-\item{method}{How to choose the best chain. Choices are c("perplexity", "harmonic", "loglik"). Defaults to perplexity. "perplexity" calculates each chain's perplexity as the inverse of the geometric mean, per the original LDA description. "harmonic" calculates each chain's marginal likelihood as the harmonic mean of each iteration of Gibbs sampling's log likelihoods. "loglik" chooses the chain which reached the maximal log likelihood during Gibbs sampling.}

+\item{method}{How to choose the best chain. Choices are c("perplexity", "loglik"). Defaults to perplexity. "perplexity" calculates each chain's perplexity as the inverse of the geometric mean, per the original LDA description. "loglik" chooses the chain which reached the maximal log likelihood during Gibbs sampling.}

 }

 \description{

 Determine the best chain among a set of celda_* objects with

@@ -17,7 +17,7 @@ getModel(celda.list, K = NULL, L = NULL, chain = 1, best = NULL)

 \item{chain}{The desired chain for the specified model}

-\item{best}{Method for choosing best chain automatically. Options are c("perplexity", "harmonic", "loglik"). See documentation for chooseBestModel for details. Overrides chain parameter if provided.}

+\item{best}{Method for choosing best chain automatically. Options are c("perplexity", "loglik"). See documentation for chooseBestModel for details. Overrides chain parameter if provided.}

 }

 \value{

 A celda model object matching the provided parameters (of class "celda_C", "celda_G", "celda_CG" accordingly), or NA if one is not found.

@@ -10,7 +10,7 @@ render_interactive_kl_plot(celda.list, method = "perplexity",

 \arguments{

 \item{celda.list}{A list of celda_CG objects returned from celda function}

-\item{method}{One of "perplexity", "harmonic", or "loglik", passed through to calculate_performance_metric()}

+\item{method}{One of "perplexity" or "loglik", passed through to calculate_performance_metric()}

 \item{title}{The plot title}

 }

@@ -7,6 +7,8 @@

 simulateCells(model, ...)

 }

 \arguments{

+\item{model}{The celda generative model to use (one of celda_C, celda_G, celda_CG)}

+

 \item{S}{Total number of samples (celda_C, celda_CG)}

 \item{C}{The number of cells (celda_G)}

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

 \arguments{

 \item{celda.list}{A celda_list object returned from celda()}

-\item{method}{One of "perplexity", "harmonic", or "loglik"}

+\item{method}{One of "perplexity", "loglik"}

 \item{title}{Title for the plot}

@@ -11,11 +11,11 @@

 \arguments{

 \item{celda.list}{A list of celda_CG objects returned from celda function}

-\item{method}{One of "perplexity", "harmonic", or "loglik"}

+\item{method}{One of "perplexity" or "loglik"}

 \item{title}{Title for the visualize_model_performance}

-\item{log}{Set log to TRUE to visualize the log(perplexity) of Celda_CG objects. Does not work for "harmonic" metric}

+\item{log}{Set log to TRUE to visualize the log(perplexity) of Celda_CG objects.}

 }

 \description{

 Plot the performance of a list of celda_CG models returned 

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

 \arguments{

 \item{celda.list}{A celda_list object returned from celda()}

-\item{method}{One of "perplexity", "harmonic", or "loglik"}

+\item{method}{One of "perplexity" or "loglik"}

 \item{title}{Title for the plot}