From: Benjamin Auder <benjamin.auder@somewhere>
Date: Fri, 21 Apr 2017 13:33:07 +0000 (+0200)
Subject: fix for m==1
X-Git-Url: https://git.auder.net/variants/Chakart/doc/scripts/%3C?a=commitdiff_plain;h=ea5860f1b4fc91f06e371a0b26915198474a849d;p=valse.git

fix for m==1
---

diff --git a/pkg/DESCRIPTION b/pkg/DESCRIPTION
index e18fddb..3b33e25 100644
--- a/pkg/DESCRIPTION
+++ b/pkg/DESCRIPTION
@@ -39,3 +39,4 @@ Collate:
     'EMGLLF.R'
     'generateXY.R'
     'A_NAMESPACE.R'
+    'util.R'
diff --git a/pkg/R/EMGLLF.R b/pkg/R/EMGLLF.R
index 6ee7ba7..03f0a75 100644
--- a/pkg/R/EMGLLF.R
+++ b/pkg/R/EMGLLF.R
@@ -47,15 +47,20 @@ EMGLLF <- function(phiInit, rhoInit, piInit, gamInit, mini, maxi, gamma, lambda,
 .EMGLLF_R <- function(phiInit, rhoInit, piInit, gamInit, mini, maxi, gamma, lambda, 
   X, Y, eps)
 {
-  # Matrix dimensions: NOTE: phiInit *must* be an array (even if p==1)
-  n <- dim(Y)[1]
-  p <- dim(phiInit)[1]
-  m <- dim(phiInit)[2]
-  k <- dim(phiInit)[3]
+  # Matrix dimensions
+  n <- nrow(X)
+  p <- ncol(X)
+  m <- ncol(Y)
+  k <- length(piInit)
+
+  # Adjustments required when p==1 or m==1 (var.sel. or output dim 1)
+  if (p==1 || m==1)
+    phiInit <- array(phiInit, dim=c(p,m,k))
+  if (m==1)
+    rhoInit <- array(rhoInit, dim=c(m,m,k))
 
   # Outputs
-  phi <- array(NA, dim = c(p, m, k))
-  phi[1:p, , ] <- phiInit
+  phi <- phiInit
   rho <- rhoInit
   pi <- piInit
   llh <- -Inf
@@ -155,7 +160,7 @@ EMGLLF <- function(phiInit, rhoInit, piInit, gamInit, mini, maxi, gamma, lambda,
     ## E step
 
     # Precompute det(rho[,,r]) for r in 1...k
-    detRho <- sapply(1:k, function(r) det(rho[, , r]))
+    detRho <- sapply(1:k, function(r) gdet(rho[, , r]))
     sumLogLLH <- 0
     for (i in 1:n)
     {
diff --git a/pkg/R/EMGrank.R b/pkg/R/EMGrank.R
index db0b8f2..b85a0fa 100644
--- a/pkg/R/EMGrank.R
+++ b/pkg/R/EMGrank.R
@@ -46,10 +46,10 @@ matricize <- function(X)
 .EMGrank_R <- function(Pi, Rho, mini, maxi, X, Y, tau, rank)
 {
   # matrix dimensions
-  n <- dim(X)[1]
-  p <- dim(X)[2]
-  m <- dim(Rho)[2]
-  k <- dim(Rho)[3]
+  n <- nrow(X)
+  p <- ncol(X)
+  m <- ncol(Y)
+  k <- length(Pi)
 
   # init outputs
   phi <- array(0, dim = c(p, m, k))
@@ -92,7 +92,7 @@ matricize <- function(X)
       for (r in seq_len(k))
       {
         dotProduct <- tcrossprod(Y[i, ] %*% Rho[, , r] - X[i, ] %*% phi[, , r])
-        logGamIR <- log(Pi[r]) + log(det(Rho[, , r])) - 0.5 * dotProduct
+        logGamIR <- log(Pi[r]) + log(gdet(Rho[, , r])) - 0.5 * dotProduct
         # Z[i] = index of max (gam[i,])
         if (logGamIR > maxLogGamIR)
         {
diff --git a/pkg/R/computeGridLambda.R b/pkg/R/computeGridLambda.R
index c2e9c8c..cf762ec 100644
--- a/pkg/R/computeGridLambda.R
+++ b/pkg/R/computeGridLambda.R
@@ -20,9 +20,9 @@ computeGridLambda <- function(phiInit, rhoInit, piInit, gamInit, X, Y, gamma, mi
   maxi, tau, fast)
 {
   n <- nrow(X)
-  p <- dim(phiInit)[1]
-  m <- dim(phiInit)[2]
-  k <- dim(phiInit)[3]
+  p <- ncol(X)
+  m <- ncol(Y)
+  k <- length(piInit)
 
   list_EMG <- EMGLLF(phiInit, rhoInit, piInit, gamInit, mini, maxi, gamma, lambda = 0, 
     X, Y, tau, fast)
diff --git a/pkg/R/constructionModelesLassoMLE.R b/pkg/R/constructionModelesLassoMLE.R
index 90d0a2a..1275ca3 100644
--- a/pkg/R/constructionModelesLassoMLE.R
+++ b/pkg/R/constructionModelesLassoMLE.R
@@ -40,10 +40,10 @@ constructionModelesLassoMLE <- function(phiInit, rhoInit, piInit, gamInit, mini,
     if (verbose) 
       print(paste("Computations for lambda=", lambda))
 
-    n <- dim(X)[1]
-    p <- dim(phiInit)[1]
-    m <- dim(phiInit)[2]
-    k <- dim(phiInit)[3]
+    n <- nrow(X)
+    p <- ncol(X)
+    m <- ncol(Y)
+    k <- length(piInit)
     sel.lambda <- S[[lambda]]$selected
     # col.sel = which(colSums(sel.lambda)!=0) #if boolean matrix
     col.sel <- which(sapply(sel.lambda, length) > 0)  #if list of selected vars
@@ -51,8 +51,8 @@ constructionModelesLassoMLE <- function(phiInit, rhoInit, piInit, gamInit, mini,
       return(NULL)
 
     # lambda == 0 because we compute the EMV: no penalization here
-    res <- EMGLLF(array(phiInit[col.sel, , ],dim=c(length(col.sel),m,k)), rhoInit,
-      piInit, gamInit, mini, maxi, gamma, 0, as.matrix(X[, col.sel]), Y, eps, fast)
+    res <- EMGLLF(array(phiInit,dim=c(p,m,k))[col.sel, , ], rhoInit, piInit, gamInit,
+      mini, maxi, gamma, 0, as.matrix(X[, col.sel]), Y, eps, fast)
 
     # Eval dimension from the result + selected
     phiLambda2 <- res$phi
@@ -71,7 +71,7 @@ constructionModelesLassoMLE <- function(phiInit, rhoInit, piInit, gamInit, mini,
       {
         delta <- (Y %*% rhoLambda[, , r] - (X[, col.sel] %*% t(phiLambda[col.sel, , r])))
       } else delta <- (Y %*% rhoLambda[, , r] - (X[, col.sel] %*% phiLambda[col.sel, , r]))
-      densite <- densite + piLambda[r] * det(rhoLambda[, , r])/(sqrt(2 * base::pi))^m * 
+      densite <- densite + piLambda[r] * gdet(rhoLambda[, , r])/(sqrt(2 * base::pi))^m * 
         exp(-diag(tcrossprod(delta))/2)
     }
     llhLambda <- c(sum(log(densite)), (dimension + m + 1) * k - 1)
diff --git a/pkg/R/constructionModelesLassoRank.R b/pkg/R/constructionModelesLassoRank.R
index 85685e9..dc88f67 100644
--- a/pkg/R/constructionModelesLassoRank.R
+++ b/pkg/R/constructionModelesLassoRank.R
@@ -21,9 +21,9 @@
 constructionModelesLassoRank <- function(S, k, mini, maxi, X, Y, eps, rank.min, rank.max, 
   ncores, fast, verbose)
 {
-  n <- dim(X)[1]
-  p <- dim(X)[2]
-  m <- dim(Y)[2]
+  n <- nrow(X)
+  p <- ncol(X)
+  m <- ncol(Y)
   L <- length(S)
 
   # Possible interesting ranks
diff --git a/pkg/R/initSmallEM.R b/pkg/R/initSmallEM.R
index 056d7e7..44b4b06 100644
--- a/pkg/R/initSmallEM.R
+++ b/pkg/R/initSmallEM.R
@@ -10,9 +10,9 @@
 #' @importFrom stats cutree dist hclust runif
 initSmallEM <- function(k, X, Y, fast)
 {
-  n <- nrow(Y)
-  m <- ncol(Y)
+  n <- nrow(X)
   p <- ncol(X)
+  m <- ncol(Y)
   nIte <- 20
   Zinit1 <- array(0, dim = c(n, nIte))
   betaInit1 <- array(0, dim = c(p, m, k, nIte))
@@ -55,7 +55,7 @@ initSmallEM <- function(k, X, Y, fast)
         dotProduct <- tcrossprod(Y[i, ] %*% rhoInit1[, , r, repet]
           - X[i, ] %*% phiInit1[, , r, repet])
         Gam[i, r] <- piInit1[repet, r] * 
-          det(rhoInit1[, , r, repet]) * exp(-0.5 * dotProduct)
+          gdet(rhoInit1[, , r, repet]) * exp(-0.5 * dotProduct)
       }
       sumGamI <- sum(Gam[i, ])
       gamInit1[i, , repet] <- Gam[i, ]/sumGamI
diff --git a/pkg/R/main.R b/pkg/R/main.R
index fecf519..e741d65 100644
--- a/pkg/R/main.R
+++ b/pkg/R/main.R
@@ -31,9 +31,9 @@ valse <- function(X, Y, procedure = "LassoMLE", selecMod = "DDSE", gamma = 1, mi
   ncores_inner = 1, thresh = 1e-08, size_coll_mod = 10, fast = TRUE, verbose = FALSE, 
   plot = TRUE)
 {
-  p <- dim(X)[2]
-  m <- dim(Y)[2]
-  n <- dim(X)[1]
+  n <- nrow(X)
+  p <- ncol(X)
+  m <- ncol(Y)
 
   if (verbose) 
     print("main loop: over all k and all lambda")
@@ -138,7 +138,7 @@ valse <- function(X, Y, procedure = "LassoMLE", selecMod = "DDSE", gamma = 1, mi
     for (r in 1:length(modelSel$pi))
     {
       sqNorm2 <- sum((Y[i, ] %*% modelSel$rho[, , r] - X[i, ] %*% modelSel$phi[, , r])^2)
-      Gam[i, r] <- modelSel$pi[r] * exp(-0.5 * sqNorm2) * det(modelSel$rho[, , r])
+      Gam[i, r] <- modelSel$pi[r] * exp(-0.5 * sqNorm2) * gdet(modelSel$rho[, , r])
     }
   }
   Gam <- Gam/rowSums(Gam)
diff --git a/pkg/R/selectVariables.R b/pkg/R/selectVariables.R
index bfe4042..d863a4b 100644
--- a/pkg/R/selectVariables.R
+++ b/pkg/R/selectVariables.R
@@ -37,15 +37,22 @@ selectVariables <- function(phiInit, rhoInit, piInit, gamInit, mini, maxi, gamma
     params <- EMGLLF(phiInit, rhoInit, piInit, gamInit, mini, maxi, gamma, lambda, 
       X, Y, eps, fast)
 
-    p <- dim(phiInit)[1]
-    m <- dim(phiInit)[2]
+    p <- ncol(X)
+    m <- ncol(Y)
 
     # selectedVariables: list where element j contains vector of selected variables
     # in [1,m]
     selectedVariables <- lapply(1:p, function(j) {
       # from boolean matrix mxk of selected variables obtain the corresponding boolean
       # m-vector, and finally return the corresponding indices
-      seq_len(m)[apply(abs(params$phi[j, , ]) > thresh, 1, any)]
+      if (m>1) {
+        seq_len(m)[apply(abs(params$phi[j, , ]) > thresh, 1, any)]
+      } else {
+        if (any(params$phi[j, , ] > thresh))
+          1
+        else
+          numeric(0)
+      }
     })
 
     list(selected = selectedVariables, Rho = params$rho, Pi = params$pi)
diff --git a/pkg/R/util.R b/pkg/R/util.R
new file mode 100644
index 0000000..f8b01cc
--- /dev/null
+++ b/pkg/R/util.R
@@ -0,0 +1,7 @@
+# ...
+gdet <- function(M)
+{
+	if (is.matrix(M))
+		return (det(M))
+	return (M[1]) #numeric, double
+}