Fix integration bug (thx Claude)

diff --git a/pkg/DESCRIPTION b/pkg/DESCRIPTION
index 055c203..12fef64 100644 (file)
--- a/pkg/DESCRIPTION
+++ b/pkg/DESCRIPTION
@@ -5,9 +5,9 @@ Description: Mixture of logistic regressions parameters (H)estimation with
     (U)spectral methods. The main methods take d-dimensional inputs and a vector
     of binary outputs, and return parameters according to the GLMs mixture model
     (General Linear Model). For more details see chapter 3 in the PhD thesis of
-               Mor-Absa Loum: <https://theses.fr/s156435>, available here
+               Mor-Absa Loum: <https://theses.fr/2018SACLS299>, available here
                <https://theses.hal.science/tel-01877796/document>.
-Version: 1.0-5
+Version: 1.0-6
 Authors@R: c(person(given = "Benjamin",
                     family = "Auder",
                     role = c("aut", "cre"),
@@ -22,7 +22,8 @@ Imports:
     MASS,
     jointDiag,
     methods,
-    pracma
+    pracma,
+    utils
 Suggests:
     devtools,
     flexmix,
@@ -31,7 +32,6 @@ Suggests:
     roxygen2
 License: MIT + file LICENSE
 Date: 2026-04-22
-RoxygenNote: 7.3.3
 URL: https://github.com/yagu0/morpheus
 Collate:
     'utils.R'

diff --git a/pkg/NAMESPACE b/pkg/NAMESPACE
index 75c3ec6..5159028 100644 (file)
--- a/pkg/NAMESPACE
+++ b/pkg/NAMESPACE
@@ -10,9 +10,26 @@ export(optimParams)
 export(plotBox)
 export(plotCoefs)
 export(plotHist)
-importFrom(graphics,barplot,boxplot,hist,matplot,par)
+export(pvalue)
+importFrom(graphics,
+  barplot,
+  boxplot,
+  hist,
+  matplot,
+  par
+)
 importFrom(jointDiag,ajd)
 importFrom(methods,new)
 importFrom(pracma,integral)
-importFrom(stats,integrate,pnorm,rbinom,rmultinom,rnorm,runif,sd,pchisq,var)
+importFrom(stats,
+  integrate,
+  pchisq,
+  pnorm,
+  rbinom,
+  rmultinom,
+  rnorm,
+  runif,
+  sd,
+  var
+)
 useDynLib(morpheus, .registration = TRUE)

diff --git a/pkg/R/optimParams.R b/pkg/R/optimParams.R
index c050e63..cd90258 100644 (file)
--- a/pkg/R/optimParams.R
+++ b/pkg/R/optimParams.R
@@ -143,7 +143,6 @@ setRefClass(
     {
       "Compute the weights matrix from a parameters list"
 
-      require(MASS)
       dd <- d + d^2 + d^3
       M <- Moments(θ)
       Omega <- matrix( .C("Compute_Omega",
@@ -321,23 +320,26 @@ setRefClass(
   # NOTE: weird "integral divergent" error on inputs:
   # link="probit"; order=2; λ=c(531.8099,586.8893,523.5816); b=c(-118.512674,-3.488020,2.109969)
   # Switch to pracma package for that (but it seems slow...)
-  sapply( seq_along(λ), function(k) {
-    res <- NULL
-    tryCatch({
+  vapply( seq_along(λ), function(k) {
+    integrand <- function(z)
+      .deriv[[link]][[order]](λ[k]*z+b[k]) * exp(-z^2/2) / sqrt(2*pi)
+    # NOTE: the value must be *returned* by tryCatch (an assignment inside the
+    # error handler would only modify a local variable of the handler).
+    tryCatch(
       # Fast code, may fail:
-      res <- stats::integrate(
-        function(z) .deriv[[link]][[order]](λ[k]*z+b[k]) * exp(-z^2/2) / sqrt(2*pi),
-        lower=-Inf, upper=Inf )$value
-    }, error = function(e) {
-      # Robust slow code, no fails observed:
-      sink("/dev/null") #pracma package has some useless printed outputs...
-      res <- pracma::integral(
-        function(z) .deriv[[link]][[order]](λ[k]*z+b[k]) * exp(-z^2/2) / sqrt(2*pi),
-        xmin=-Inf, xmax=Inf, method="Kronrod")
-      sink()
-    })
-    res
-  })
+      stats::integrate(integrand, lower=-Inf, upper=Inf)$value,
+      error = function(e) {
+        # Robust slow code, no fails observed.
+        # pracma package has some useless printed outputs: silence them
+        # without touching the user's sink() stack.
+        out <- NULL
+        utils::capture.output({
+          out <- pracma::integral(integrand, xmin=-Inf, xmax=Inf, method="Kronrod")
+        })
+        out
+      }
+    )
+  }, numeric(1) )
 }
 
 # Derivatives list: g^(k)(x) for links 'logit' and 'probit'

diff --git a/pkg/R/plot.R b/pkg/R/plot.R
index f76017d..5d34ce1 100644 (file)
--- a/pkg/R/plot.R
+++ b/pkg/R/plot.R
@@ -123,9 +123,6 @@ plotBox <- function(mr, x, y, ...)
 #' @export
 plotCoefs <- function(mr, params, ...)
 {
-  oldpar <- par(no.readonly = TRUE)
-  on.exit(par(oldpar))
-
   d <- nrow(mr[[1]])
   K <- ncol(mr[[1]])
 
@@ -142,7 +139,11 @@ plotCoefs <- function(mr, params, ...)
     }
   }
 
-  par(cex.axis=1.5, cex.lab=1.5, mar=c(4.7,5,1,1))
+  oldpar <- par(cex.axis=1.5, cex.lab=1.5, mar=c(4.7,5,1,1))
+  # Save/restore only the parameters modified here: restoring the full par()
+  # would also reset the current panel (mfg) and break par(mfrow=...) layouts.
+  on.exit(par(oldpar))
+
   params <- as.double(params)
   o <- order(params)
   avg_param <- as.double(params_hat)

diff --git a/pkg/inst/doc/JSS_article.pdf b/pkg/inst/doc/JSS_article.pdf
index 0b930dc..a9d1cc4 100644 (file)
--- a/pkg/inst/doc/JSS_article.pdf
+++ b/pkg/inst/doc/JSS_article.pdf
@@ -1 +1 @@
-#$# git-fat 51b69cedf9b3a4397c353799b5e03a0b6254219a               404201
+#$# git-fat 925bd9ff026245c6c9f5ff805ca5eb29cca3ad70               268250