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

the_data <- import_biom("data/ravel/ravel.biom")


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]]))))

fitted_model <- BM_poisson_covariates("LBM", adj = per_taxa_networks[[2]], covariates = list(covar_mat))
fitted_model$estimate()

fitted_model$model_parameters[[which.max(fitted_model$ICL)]]