[Swift-commit] r4037 - in SwiftApps/SwiftR: . mxtests

[noreply at svn.ci.uchicago.edu]

Author: tga
Date: 2011-01-25 16:22:33 -0600 (Tue, 25 Jan 2011)
New Revision: 4037

Added:
   SwiftApps/SwiftR/mxtests/
   SwiftApps/SwiftR/mxtests/BootstrapParallelBigger.R
Log:
Checking in initial version of larger bootstrap parallel script before modifying to work
with swiftR


Added: SwiftApps/SwiftR/mxtests/BootstrapParallelBigger.R
===================================================================
--- SwiftApps/SwiftR/mxtests/BootstrapParallelBigger.R	                        (rev 0)
+++ SwiftApps/SwiftR/mxtests/BootstrapParallelBigger.R	2011-01-25 22:22:33 UTC (rev 4037)
@@ -0,0 +1,118 @@
+#
+#   Copyright 2007-2010 The OpenMx Project
+#
+#   Licensed under the Apache License, Version 2.0 (the "License");
+#   you may not use this file except in compliance with the License.
+#   You may obtain a copy of the License at
+# 
+#        http://www.apache.org/licenses/LICENSE-2.0
+# 
+#   Unless required by applicable law or agreed to in writing, software
+#   distributed under the License is distributed on an "AS IS" BASIS,
+#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+#   See the License for the specific language governing permissions and
+#   limitations under the License.
+
+require(OpenMx)
+require(snowfall)
+
+hostfile <- Sys.getenv(c("PBS_NODEFILE"))
+
+hostnames <- read.table(hostfile, stringsAsFactors=FALSE)[[1]]
+ncpus <- length(hostnames)
+
+sfSetMaxCPUs(ncpus)
+sfInit(parallel=TRUE, cpus=ncpus, socketHosts=hostnames)
+sfLibrary(OpenMx)
+
+set.seed(10)
+
+# parameters for the simulation: lambda = factor loadings,
+# specifics = specific variances
+nVar <- 75
+nObs <- 10000
+nReps <- 1024
+goodStartValues <- TRUE
+if (!is.logical(goodStartValues)) {
+  stop("'goodStartValues' should be logical. Try again.")
+}
+
+lambda <- matrix(1:nVar*2/nVar, nVar, 1)
+specifics <- diag(nVar)
+chl <- chol(lambda %*% t(lambda) + specifics)
+
+# indices for parameters and hessian estimate in results
+pStrt <- 3
+pEnd <- pStrt + 2*nVar - 1
+hStrt <- pEnd + 1
+hEnd <- hStrt + 2*nVar - 1
+
+# dimension names for OpenMx
+dn <- list()
+dn[[1]] <- paste("Var", 1:nVar, sep="")
+dn[[2]] <- dn[[1]]
+
+# function to get a covariance matrix
+randomCov <- function(nObs, nVar, chl, dn) {
+  x <- matrix(rnorm(nObs*nVar), nObs, nVar)
+  x <- x %*% chl
+  thisCov <- cov(x)
+  dimnames(thisCov) <- dn
+  return(thisCov)  
+}
+
+createNewModel <- function(index, prefix, model) {
+	modelname <- paste(prefix, index, sep='')
+	model at data@observed <- randomCov(nObs, nVar, chl, dn)
+	model at name <- modelname
+	return(model)
+}
+
+getStats <- function(model) {
+	retval <- c(model at output$status[[1]],
+		max(abs(model at output$gradient)),
+		model at output$estimate,
+		sqrt(diag(solve(model at output$hessian))))
+	return(retval)
+}
+
+
+# initialize obsCov for MxModel
+obsCov <- randomCov(nObs, nVar, chl, dn)
+
+# results matrix: get results for each simulation
+results <- matrix(0, nReps, hEnd)
+dnr <- c("inform", "maxAbsG", paste("lambda", 1:nVar, sep=""),
+         paste("specifics", 1:nVar, sep=""),
+         paste("hessLambda", 1:nVar, sep=""),
+         paste("hessSpecifics", 1:nVar, sep=""))
+dimnames(results)[[2]] <- dnr
+
+# instantiate MxModel
+template <- mxModel(name="stErrSim",
+                       mxMatrix(name="lambda", type="Full", nrow=nVar, ncol=1,
+                                free=TRUE, values=(1:nVar*2/nVar)*goodStartValues),
+                       mxMatrix(name="specifics", type="Diag", nrow=nVar,
+                                free=TRUE, values=rep(1, nVar)),
+                       mxAlgebra(lambda %*% t(lambda) + specifics,
+                                 name="preCov", dimnames=dn),
+                       mxData(observed=obsCov, type="cov", numObs=nObs),
+                       mxMLObjective(covariance='preCov'),
+                       independent = TRUE)
+
+topModel <- mxModel(name = 'container')
+
+sfExportAll()
+
+submodels <- lapply(1:nReps, createNewModel, 'stErrSim', template)
+
+names(submodels) <- omxExtractNames(submodels)
+topModel at submodels <- submodels
+
+modelResults <- mxRun(topModel, silent=TRUE, suppressWarnings=TRUE)
+
+print(ncpus)
+print(modelResults at output$wallTime)
+
+sfStop()
+