library(blockmodels)
library(phyloseq)
library(biomformat)
library(here)
source("utils.R")

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[[level]]) %*% t(rep(1, nrow(per_taxa_networks[[level]]))))

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)