X-Git-Url: https://git.auder.net/?p=morpheus.git;a=blobdiff_plain;f=pkg%2FR%2FcomputeMu.R;h=bc52bb31df714ab1146c279eb1f22d239cf40105;hp=2a62793ef756b8f8af5c55f59b96a2fb8b5598a5;hb=2b3a6af5c55ac121405e3a8da721626ddf46b28b;hpb=d294ece1cf943b74d96b26cc28b08c00cb191264

diff --git a/pkg/R/computeMu.R b/pkg/R/computeMu.R
index 2a62793..bc52bb3 100644
--- a/pkg/R/computeMu.R
+++ b/pkg/R/computeMu.R
@@ -23,70 +23,69 @@
 #' @examples
 #' io = generateSampleIO(10000, 1/2, matrix(c(1,0,0,1),ncol=2), c(0,0), "probit")
 #' μ = computeMu(io$X, io$Y, list(K=2)) #or just X and Y for estimated K
+#'
 #' @export
 computeMu = function(X, Y, optargs=list())
 {
-	if (!is.matrix(X) || !is.numeric(X) || any(is.na(X)))
-		stop("X: real matrix, no NA")
-	n = nrow(X)
-	d = ncol(X)
-	if (!is.numeric(Y) || length(Y)!=n || any(Y!=0 & Y!=1))
-		stop("Y: vector of 0 and 1, size nrow(X), no NA")
-	if (!is.list(optargs))
-		stop("optargs: list")
+  if (!is.matrix(X) || !is.numeric(X) || any(is.na(X)))
+    stop("X: real matrix, no NA")
+  n = nrow(X)
+  d = ncol(X)
+  if (!is.numeric(Y) || length(Y)!=n || any(Y!=0 & Y!=1))
+    stop("Y: vector of 0 and 1, size nrow(X), no NA")
+  if (!is.list(optargs))
+    stop("optargs: list")
 
-	# Step 0: Obtain the empirically estimated moments tensor, estimate also K
-	M = if (is.null(optargs$M)) computeMoments(X,Y) else optargs$M
-	K = optargs$K
-	if (is.null(K))
-	{
-		# TODO: improve this basic heuristic
-		Σ = svd(M[[2]])$d
-		large_ratio <- ( abs(Σ[-d] / Σ[-1]) > 3 )
-		K <- if (any(large_ratio)) max(2, which.min(large_ratio)) else d
-	}
+  # Step 0: Obtain the empirically estimated moments tensor, estimate also K
+  M = if (is.null(optargs$M)) computeMoments(X,Y) else optargs$M
+  K = optargs$K
+  if (is.null(K))
+  {
+    # TODO: improve this basic heuristic
+    Σ = svd(M[[2]])$d
+    large_ratio <- ( abs(Σ[-d] / Σ[-1]) > 3 )
+    K <- if (any(large_ratio)) max(2, which.min(large_ratio)) else d
+  }
 
-	# Step 1: generate a family of d matrices to joint-diagonalize to increase robustness
-	d = ncol(X)
-	fixed_design = FALSE
-	jd_nvects = ifelse(!is.null(optargs$jd_nvects), optargs$jd_nvects, 0)
-	if (jd_nvects == 0)
-	{
-		jd_nvects = d
-		fixed_design = TRUE
-	}
-	M2_t = array(dim=c(d,d,jd_nvects))
-	for (i in seq_len(jd_nvects))
-	{
-		rho = if (fixed_design) c(rep(0,i-1),1,rep(0,d-i)) else normalize( rnorm(d) )
-		M2_t[,,i] = .T_I_I_w(M[[3]],rho)
-	}
+  # Step 1: generate a family of d matrices to joint-diagonalize to increase robustness
+  d = ncol(X)
+  fixed_design = FALSE
+  jd_nvects = ifelse(!is.null(optargs$jd_nvects), optargs$jd_nvects, 0)
+  if (jd_nvects == 0)
+  {
+    jd_nvects = d
+    fixed_design = TRUE
+  }
+  M2_t = array(dim=c(d,d,jd_nvects))
+  for (i in seq_len(jd_nvects))
+  {
+    rho = if (fixed_design) c(rep(0,i-1),1,rep(0,d-i)) else normalize( rnorm(d) )
+    M2_t[,,i] = .T_I_I_w(M[[3]],rho)
+  }
 
-	# Step 2: obtain factors u_i (and their inverse) from the joint diagonalisation of M2_t
-	jd_method = ifelse(!is.null(optargs$jd_method), optargs$jd_method, "uwedge")
-	V =
-		if (jd_nvects > 1) {
-			#NOTE: increasing itermax does not help to converge, thus we suppress warnings
-			suppressWarnings({jd = jointDiag::ajd(M2_t, method=jd_method)})
-#			if (jd_method=="uwedge") jd$B else solve(jd$A)
-			if (jd_method=="uwedge") jd$B else MASS::ginv(jd$A)
-		}
-		else
-			eigen(M2_t[,,1])$vectors
+  # Step 2: obtain factors u_i (and their inverse) from the joint diagonalisation of M2_t
+  jd_method = ifelse(!is.null(optargs$jd_method), optargs$jd_method, "uwedge")
+  V =
+    if (jd_nvects > 1) {
+      # NOTE: increasing itermax does not help to converge, thus we suppress warnings
+      suppressWarnings({jd = jointDiag::ajd(M2_t, method=jd_method)})
+      if (jd_method=="uwedge") jd$B else MASS::ginv(jd$A)
+    }
+    else
+      eigen(M2_t[,,1])$vectors
 
-	# Step 3: obtain final factors from joint diagonalisation of T(I,I,u_i)
-	M2_t = array(dim=c(d,d,K))
-	for (i in seq_len(K))
-		M2_t[,,i] = .T_I_I_w(M[[3]],V[,i])
-	suppressWarnings({jd = jointDiag::ajd(M2_t, method=jd_method)})
-#	U = if (jd_method=="uwedge") solve(jd$B) else jd$A
-	U = if (jd_method=="uwedge") MASS::ginv(jd$B) else jd$A
-	μ = normalize(U[,1:K])
+  # Step 3: obtain final factors from joint diagonalisation of T(I,I,u_i)
+  M2_t = array(dim=c(d,d,K))
+  for (i in seq_len(K))
+    M2_t[,,i] = .T_I_I_w(M[[3]],V[,i])
+  suppressWarnings({jd = jointDiag::ajd(M2_t, method=jd_method)})
+  U = if (jd_method=="uwedge") MASS::ginv(jd$B) else jd$A
+  μ = normalize(U[,1:K])
 
-	# M1 also writes M1 = sum_k coeff_k * μ_k, where coeff_k >= 0
-	# ==> search decomposition of vector M1 onto the (truncated) basis μ (of size dxK)
-	# This is a linear system μ %*% C = M1 with C of size K ==> C = psinv(μ) %*% M1
-	C = MASS::ginv(μ) %*% M[[1]]
-	μ[,C < 0] = - μ[,C < 0]
-	μ
+  # M1 also writes M1 = sum_k coeff_k * μ_k, where coeff_k >= 0
+  # ==> search decomposition of vector M1 onto the (truncated) basis μ (of size dxK)
+  # This is a linear system μ %*% C = M1 with C of size K ==> C = psinv(μ) %*% M1
+  C = MASS::ginv(μ) %*% M[[1]]
+  μ[,C < 0] = - μ[,C < 0]
+  μ
 }