BM correct sampling

blockmodels_correct_sampling_effort.R

@@ -1,16 +1,45 @@
 library(blockmodels)
 library(phyloseq)
 library(biomformat)
+library(here)
 source("utils.R")
 
-the_data <- import_biom("data/ravel/ravel.biom")
+options(parallelly.availableCores.custom = function() {
+    ncores <- min(parallelly::availableCores(), 64)
+    max(1L, ncores)
+})
 
+message(paste("Number of cores available:", parallelly::availableCores()))
+args <- commandArgs(trailingOnly = TRUE)
+if (length(args) == 0) {
+    author <- "mach"
+} else {
+    author <- args[1]
+}
+
+message("For author: ", author)
+the_data <- import_biom(here("data", author, paste0(author, ".biom")))
+
+epoch <- as.integer(Sys.time())
+level <- 2
+message("At level: ", level, " starting at: ", epoch)
+
+save_path <- here("results", "seq-effort", paste0(author, "_", level, "_", epoch, ".Rds"))
 
 per_taxa_networks <- collapse_otu_at_taxo(the_data)
 
-covar_mat <- t(colSums(per_taxa_networks[[2]]) %*% t(rep(1, nrow(per_taxa_networks[[2]]))))
+covar_mat <- t(colSums(per_taxa_networks[[level]]) %*% t(rep(1, nrow(per_taxa_networks[[level]]))))
 
-fitted_model <- BM_poisson_covariates("LBM", adj = per_taxa_networks[[2]], covariates = list(covar_mat))
+lcovar_mat <- log(covar_mat)
+
+fitted_model <- BM_poisson_covariates("LBM",
+    adj = per_taxa_networks[[level]],
+    covariates = list(lcovar_mat),
+    autosave = save_path,
+    ncores = parallelly::availableCores()
+)
 fitted_model$estimate()
 
 fitted_model$model_parameters[[which.max(fitted_model$ICL)]]
+message("Saving the model results.")
+saveRDS(fitted_model, save_path)