From: Benjamin Auder <benjamin.auder@somewhere>
Date: Mon, 28 Oct 2019 21:38:47 +0000 (+0100)
Subject: Fix default weights, fix reporting scripts
X-Git-Url: https://git.auder.net/?p=morpheus.git;a=commitdiff_plain;h=0ad4c8de650e9f27ec3754c9cb9b2a03db5aff24

Fix default weights, fix reporting scripts
---

diff --git a/pkg/R/optimParams.R b/pkg/R/optimParams.R
index 85c21e7..06d1684 100644
--- a/pkg/R/optimParams.R
+++ b/pkg/R/optimParams.R
@@ -59,7 +59,7 @@ optimParams = function(K, link=c("logit","probit"), optargs=list())
 
   weights <- optargs$weights
   if (is.null(weights))
-    weights <- rep(1, K)
+    weights <- rep(1, 3)
 
 	# Build and return optimization algorithm object
 	methods::new("OptimParams", "li"=link, "M1"=as.double(M[[1]]),
diff --git a/reports/accuracy.R b/reports/accuracy.R
index e9c9d1b..91d9c61 100644
--- a/reports/accuracy.R
+++ b/reports/accuracy.R
@@ -11,11 +11,10 @@ optimBeta <- function(N, n, K, p, beta, b, link, weights, ncores)
 				M <- computeMoments(fargs$X, fargs$Y)
 				fargs$optargs$M <- M
 				mu <- computeMu(fargs$X, fargs$Y, fargs$optargs)
-				res2 <- NULL
+        op <- optimParams(K,fargs$optargs$link,fargs$optargs)
+        x_init <- list(p=rep(1/K,K-1), beta=mu, b=rep(0,K))
 				tryCatch({
-					op <- optimParams(K,fargs$optargs$link,fargs$optargs)
-					x_init <- list(p=rep(1/K,K-1), beta=mu, b=rep(0,K))
-					res2 <- do.call(rbind, op$run(x_init))
+          res2 <- do.call(rbind, op$run(x_init))
 				}, error = function(e) {
 					res2 <- NA
 				})
@@ -50,7 +49,6 @@ optimBeta <- function(N, n, K, p, beta, b, link, weights, ncores)
 			fargs$X = io$X
 			fargs$Y = io$Y
 			fargs$optargs$K = ncol(fargs$beta)
-			fargs$optargs$M = computeMoments(io$X,io$Y)
 			fargs
 		}, N=N, ncores=ncores, verbose=TRUE)
 	p <- c(p, 1-sum(p))
@@ -119,7 +117,7 @@ if (d == 2) {
 }
 
 mr <- optimBeta(N, n, K, p, beta, b, link, weights, ncores)
-mr_params <- list("N"=N, "n"=n, "K"=K, "d"=d, "link"=link,
+mr_params <- list("N"=N, "nc"=ncores, "n"=n, "K"=K, "d"=d, "link"=link,
 	"p"=c(p,1-sum(p)), "beta"=beta, "b"=b, "weights"=weights)
 
 save("mr", "mr_params", file=paste("res_",n,"_",d,"_",link,"_",strw,".RData",sep=""))
diff --git a/reports/multistart.R b/reports/multistart.R
index 9971cd2..5fa80ad 100644
--- a/reports/multistart.R
+++ b/reports/multistart.R
@@ -1,12 +1,76 @@
 library(morpheus)
 
+testMultistart <- function(N, n, K, p, beta, b, link, nstart, ncores)
+{
+  ms <- multiRun(
+    list(n=n,p=p,beta=beta,b=b,optargs=list(K=K,link=link,nstart=nstart)),
+    list(
+      function(fargs) {
+        # 1 start
+        library(morpheus)
+        K <- fargs$optargs$K
+        op <- optimParams(K, fargs$optargs$link, fargs$optargs)
+        x_init <- list(p=rep(1/K,K-1), beta=fargs$mu, b=rep(0,K))
+				res <- NULL
+				tryCatch({
+          res <- do.call(rbind, op$run(x_init))
+				}, error = function(e) {
+					res <- NA
+				})
+				res
+      },
+      function(fargs) {
+        # B starts
+        library(morpheus)
+        K <- fargs$optargs$K
+        op <- optimParams(K, fargs$optargs$link, fargs$optargs)
+        best_val <- Inf
+        best_par <- list()
+        for (i in 1:fargs$optargs$nstart)
+        {
+          #x_init <- list(p=rep(1/K,K-1), beta=i*fargs$mu, b=rep(0,K))
+          M <- matrix(rnorm(d*K), nrow=d, ncol=K)
+          M <- t(t(M) / sqrt(colSums(M^2)))
+          x_init <- list(p=rep(1/K,K-1), beta=M, b=rep(0,K))
+          tryCatch({
+            par <- op$run(x_init)
+          }, error = function(e) {
+            par <- NA
+          })
+          if (!is.na(par[0]))
+          {
+            val <- op$f( op$linArgs(par) )
+            if (val < best_val)
+            {
+              best_par <- par
+              best_val <- val
+            }
+          }
+        }
+        # Bet that at least one run succeded:
+        do.call(rbind,best_par)
+      }
+    ),
+    prepareArgs = function(fargs, index) {
+      library(morpheus)
+      io = generateSampleIO(fargs$n, fargs$p, fargs$beta, fargs$b, fargs$optargs$link)
+      fargs$optargs$M <- computeMoments(io$X, io$Y)
+      mu <- computeMu(io$X, io$Y, fargs$optargs)
+      fargs$mu <- mu
+      fargs
+    }, N=N, ncores=ncores, verbose=TRUE)
+  for (i in 1:2)
+    ms[[i]] <- alignMatrices(ms[[i]], ref=rbind(p,beta,b), ls_mode="exact")
+  ms
+}
+
 #model = binomial
 K <- 2
 p <- .5
 b <- c(-.2, .5)
 # Default values:
 link = "logit"
-N <- 100
+N <- 10
 d <- 2
 n <- 1e4
 ncores <- 1
@@ -39,47 +103,8 @@ betas <- list(
 	matrix( c(1,2,-1,0,3,4,-1,-3,0,2, 2,-3,0,1,0,-1,-4,3,2,0), ncol=K ) ) #d=10
 beta <- betas[[ ifelse( d==2, 1, ifelse(d==5,2,3) ) ]]
 
-ms <- multiRun(
-	list(n=n,p=p,beta=beta,b=b,optargs=list(K=K,link=link,nstart=nstart)), list(
-		function(fargs) {
-			# 1 start
-			library(morpheus)
-			K <- fargs$optargs$K
-			op <- optimParams(K, fargs$optargs$link, fargs$optargs)
-			x_init <- list(p=rep(1/K,K-1), beta=fargs$mu, b=rep(0,K))
-			do.call(rbind,op$run(x_init))
-		},
-		function(fargs) {
-			# B starts
-			library(morpheus)
-			K <- fargs$optargs$K
-			op <- optimParams(K, fargs$optargs$link, fargs$optargs)
-			best_val <- Inf
-			best_par <- list()
-			for (i in 1:fargs$optargs$nstart)
-			{
-				x_init <- list(p=rep(1/K,K-1), beta=i*fargs$mu, b=rep(0,K))
-				par <- op$run(x_init)
-				val <- op$f( op$linArgs(par) )
-				if (val < best_val)
-				{
-					best_par <- par
-					best_val <- val
-				}
-			}
-			do.call(rbind,best_par)
-		}),
-		prepareArgs = function(fargs) {
-			library(morpheus)
-			io = generateSampleIO(fargs$n, fargs$p, fargs$beta, fargs$b, fargs$optargs$link)
-			fargs$optargs$M <- computeMoments(io$X, io$Y)
-			mu <- computeMu(io$X, io$Y, fargs$optargs)
-			fargs$mu <- mu
-		}, N=N, ncores=ncores, verbose=TRUE)
-for (i in 1:2)
-	ms[[i]] <- alignMatrices(ms[[i]], ref=rbind(p,beta,b), ls_mode="exact")
-
-ms_params <- list("N"=N, "nc"=ncores, "n"=n, "K"=K, "link"=link,
-	"p"=p, "beta"=beta, "b"=b, "nstart"=nstart)
+ms <- testMultistart(N, n, K, p, beta, b, link, nstart, ncores)
+ms_params <- list("N"=N, "nc"=ncores, "n"=n, "K"=K, "d"=d, "link"=link,
+	"p"=c(p,1-sum(p)), "beta"=beta, "b"=b, "nstart"=nstart)
 
-save(ms, ms_params, file="multistart.RData")
+save("ms", "ms_params", file="multistart.RData")
diff --git a/reports/test.R b/reports/test.R
deleted file mode 100644
index 2ce9a44..0000000
--- a/reports/test.R
+++ /dev/null
@@ -1,58 +0,0 @@
-library(morpheus)
-morph <- function(fargs) {
-	K <- fargs$optargs$K
-	M <- computeMoments(fargs$X, fargs$Y)
-	fargs$optargs$M <- M
-	mu <- computeMu(fargs$X, fargs$Y, fargs$optargs)
-	res2 <- NULL
-	tryCatch({
-		op <- optimParams(K,link,fargs$optargs)
-		x_init <- list(p=rep(1/K,K-1), beta=mu, b=rep(0,K))
-		res2 <- do.call(rbind, op$run(x_init))
-	}, error = function(e) {
-		res2 <- NA
-	})
-	res2
-}
-
-#model = binomial; default values:
-link = "probit"
-N <- 10
-d <- 2
-n <- 1e4
-ncores <- 1
-
-if (d == 2) {
-	K <- 2
-	p <- .5
-	b <- c(-.2, .5)
-	beta <- matrix( c(1,-2, 3,1), ncol=K )
-} else if (d == 5) {
-	K <- 2
-	p <- .5
-	b <- c(-.2, .5)
-	beta <- matrix( c(1,2,-1,0,3, 2,-3,0,1,0), ncol=K )
-} else if (d == 10) {
-	K <- 3
-	p <- c(.3, .3)
-	b <- c(-.2, 0, .5)
-	beta <- matrix( c(1,2,-1,0,3,4,-1,-3,0,2, 2,-3,0,1,0,-1,-4,3,2,0, -1,1,3,-1,0,0,2,0,1,-2), ncol=K )
-} else if (d == 20) {
-	K <- 3
-	p <- c(.3, .3)
-	b <- c(-.2, 0, .5)
-	beta <- matrix( c(1,2,-1,0,3,4,-1,-3,0,2,2,-3,0,1,0,-1,-4,3,2,0, -1,1,3,-1,0,0,2,0,1,-2,1,2,-1,0,3,4,-1,-3,0,2, 2,-3,0,1,0,-1,-4,3,2,0,1,1,2,2,-2,-2,3,1,0,0), ncol=K )
-}
-
-fargs = list(n=n, p=p, beta=beta, b=b)
-fargs$optargs = list(link=link)
-
-io = generateSampleIO(fargs$n, fargs$p, fargs$beta, fargs$b, fargs$optargs$link)
-fargs$X = io$X
-fargs$Y = io$Y
-fargs$optargs$K = ncol(fargs$beta)
-fargs$optargs$M = computeMoments(io$X,io$Y)
-
-res2 <- morph(fargs)
-
-save("res2", file="test.RData")
diff --git a/reports/test.sh b/reports/test.sh
deleted file mode 100644
index b617a09..0000000
--- a/reports/test.sh
+++ /dev/null
@@ -1,16 +0,0 @@
-#!/bin/bash
-
-# arg --vanilla maybe possible on cluster
-for d in 2 5; do
-	for link in "logit" "probit"; do
-		R --slave --args N=10 n=1e3 nc=3 d=$d link=$link <accuracy.R >out$d$link 2>&1
-	done
-done
-
-#for d in 2 5; do
-#	for n in 5000 10000 100000 500000 1000000; do
-#		for link in "logit" "probit"; do
-#			R --slave --args N=1000 n=$n nc=64 d=$d link=$link <accuracy.R >out_$n$link$d 2>&1
-#		done
-#	done
-#done