Fix PLNnetwork crash issues from glasso

R/PLNfamily-class.R

@@ -72,16 +72,24 @@ PLNfamily <-
         } else {
           nullModel <- NULL
         }
-        for (model in self$models)
-          model$postTreatment(
-            self$responses,
-            self$covariates,
-            self$offsets,
-            self$weights,
-            config_post=config_post,
-            config_optim=config_optim,
-            nullModel = nullModel
-        )
+        tryCatch({
+          for (i in seq_along(self$models)) {
+            model <- self$models[[i]]
+            model$postTreatment(
+              self$responses,
+              self$covariates,
+              self$offsets,
+              self$weights,
+              config_post = config_post,
+              config_optim = config_optim,
+              nullModel = nullModel
+            )
+          }
+        },
+        error = function(e) {
+          warning(paste("Post-treatment failed for model", i, ":", e$message, "\nTruncating model family to models 1 to", i - 1))
+          self$models <- self$models[1:(i - 1)]
+        })
       },
 
       ## %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

R/PLNnetworkfamily-class.R

@@ -119,6 +119,16 @@ Networkfamily <- R6Class(
 
     },
 
+    ## %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    ## Post treatment --------------------
+    #' @description Update fields after optimization
+    #' @param config_post a list for controlling the post-treatments (optional bootstrap, jackknife, R2, etc.).
+    #' @param config_optim a list for controlling the optimization parameters used during post_treatments
+    postTreatment = function(config_post, config_optim) {
+      super$postTreatment(config_post, config_optim)
+      private$params <- self$penalties[seq_along(self$models)]
+    },
+
     ## %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     ## Extractors ------------------------
     #' @description Extract the regularization path of a [`Networkfamily`]
@@ -432,7 +442,13 @@ PLNnetworkfamily <- R6Class(
         myPLN <- PLNnetworkfamily$new(self$penalties, data, control)
         myPLN$optimize(data, control$config_optim)
         nets <- do.call(cbind, lapply(myPLN$models, function(model) {
-          as.matrix(model$latent_network("support"))[upper.tri(diag(private$p))]
+          # If Omega is null, glasso diverged on the first iteration, so the network is completely unstable
+          if (is.null(model$model_par$Omega)) {
+            support <- matrix(0, nrow = private$p, ncol = private$p)
+          } else {
+            support <- as.matrix(model$latent_network("support"))
+          }
+          support[upper.tri(diag(private$p))]
         }))
         nets
       }, mc.cores = getOption("mc.cores", 1L))
@@ -595,7 +611,12 @@ ZIPLNnetworkfamily <- R6Class(
         myPLN$optimize(data, control$config_optim)
 
         nets <- do.call(cbind, lapply(myPLN$models, function(model) {
-          as.matrix(model$latent_network("support"))[upper.tri(diag(private$p))]
+          if (is.null(model$model_par$Omega)) {
+            support <- matrix(0, nrow = private$p, ncol = private$p)
+          } else {
+            support <- as.matrix(model$latent_network("support"))
+          }
+          support[upper.tri(diag(private$p))]
         }))
         nets
       }, mc.cores = getOption("mc.cores", 1L))