[Gmpm-commits] r17 - in pkg: . R man

noreply at r-forge.r-project.org noreply at r-forge.r-project.org
Thu Jun 10 19:35:17 CEST 2010


Author: dalebarr
Date: 2010-06-10 19:35:17 +0200 (Thu, 10 Jun 2010)
New Revision: 17

Removed:
   pkg/man/gmpmCtrl.Rd
Modified:
   pkg/DESCRIPTION
   pkg/NAMESPACE
   pkg/R/basemethods.R
   pkg/R/generics.R
   pkg/R/gmpm.R
   pkg/R/helpers.R
   pkg/man/GMPMSummary-class.Rd
   pkg/man/gmpm.Rd
   pkg/man/kb07.Rd
   pkg/man/mdTest.Rd
Log:
GMPM now takes advantage of multicore processing

Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION	2010-03-09 22:16:48 UTC (rev 16)
+++ pkg/DESCRIPTION	2010-06-10 17:35:17 UTC (rev 17)
@@ -1,12 +1,12 @@
 Package: gmpm
 Type: Package
 Title: Generalized Multilevel Permutation Models
-Version: 0.4-1
-Date: 2010-02-22
-Author: Dale Barr <dalejbarr3 at gmail.com>
-Maintainer: Dale Barr <dalejbarr3 at gmail.com>
+Version: 0.5-1
+Date: 2010-06-10
+Author: Dale Barr <d.barr at psy.gla.ac.uk>
+Maintainer: Dale Barr <d.barr at psy.gla.ac.uk>
 Description: Permutation methods for testing hypotheses on multilevel experimental data
 License: GPL (>=2)
 LazyLoad: yes
 Depends: methods
-Suggests: nnet
+Suggests: nnet, multicore

Modified: pkg/NAMESPACE
===================================================================
--- pkg/NAMESPACE	2010-03-09 22:16:48 UTC (rev 16)
+++ pkg/NAMESPACE	2010-06-10 17:35:17 UTC (rev 17)
@@ -1,7 +1,6 @@
 export(
        "gmpm",
-       "gmpmCreate",
-       "gmpmCtrl"
+       "gmpmCreate"
        )
 
 exportMethods(

Modified: pkg/R/basemethods.R
===================================================================
--- pkg/R/basemethods.R	2010-03-09 22:16:48 UTC (rev 16)
+++ pkg/R/basemethods.R	2010-06-10 17:35:17 UTC (rev 17)
@@ -34,6 +34,7 @@
 
                         ncomp="numeric", # n of runs completed
                         ndigits="numeric", # n of digits to round p value to
+                        nCores="integer", # N of available processing cores
                         "VIRTUAL"), # factor matrix for the model
          
          prototype=prototype(
@@ -48,7 +49,7 @@
                         gmpmRegSum="list", # main regression
                         showReg="logical" # whether to show reg coef?
                         ),
-         prototype(showReg=FALSE)
+         prototype(showReg=TRUE)
          )
 
 setClass(
@@ -309,7 +310,7 @@
 
 setMethod("summary",
     signature(object = "GMPM"),
-          function (object, showReg=FALSE, ...) 
+          function (object, showReg=TRUE, ...) 
           {
 #            print("~~~ in summary (GMPM) ~~~")
             x <- object

Modified: pkg/R/generics.R
===================================================================
--- pkg/R/generics.R	2010-03-09 22:16:48 UTC (rev 16)
+++ pkg/R/generics.R	2010-06-10 17:35:17 UTC (rev 17)
@@ -185,3 +185,7 @@
 setGeneric(name=".collapseMultinomPmx",
            def=function(x,index){
              standardGeneric(".collapseMultinomPmx")})
+
+setGeneric(name=".setOpts",
+           def=function(x,opts) {
+             standardGeneric(".setOpts")})

Modified: pkg/R/gmpm.R
===================================================================
--- pkg/R/gmpm.R	2010-03-09 22:16:48 UTC (rev 16)
+++ pkg/R/gmpm.R	2010-06-10 17:35:17 UTC (rev 17)
@@ -632,15 +632,20 @@
           function(x, ix, maxruns, elapsed) {
             fmtstr1 <- paste("%",nchar(as.character(maxruns)),"d",
                             sep="")
-            fmtstr2 <- paste(fmtstr1, "/", fmtstr1, ":", sep="")
-            stmp1 <- sprintf(fmtstr2, ix, maxruns)
+            ixx <- ifelse( (ix*x at nCores)>x at gmpmControl[["maxruns"]],
+                          x at gmpmControl[["maxruns"]], ix*x at nCores)
+            fmtstr1.2 <- paste("%", nchar(as.character(x at gmpmControl[["maxruns"]])),"d", sep="")
+            fmtstr2 <- paste(fmtstr1, "/", maxruns, " (", fmtstr1.2,
+                             "/", x at gmpmControl[["maxruns"]],
+                             ") ", ":", sep="")
+            stmp1 <- sprintf(fmtstr2, ix, ixx, maxruns)
             #stmp2 <- sprintf("% 1.3f", tmp1)
             if (length(elapsed) > 1) {
               efac <- mean(elapsed)
             } else {
               efac <- elapsed[1]
             }
-            cat(stmp1, sprintf("%1.3fs/run, ", efac))
+            cat(stmp1, sprintf("%1.3fs/sweep, ", efac))
 
             totsecs <- efac*(maxruns-ix)
             cat(sprintf("%02dh:", floor(totsecs / 3600)))
@@ -1347,16 +1352,21 @@
           function(x,gmpmControl)
           {
             #print("~~~ in gmpmFit (GMPM) ~~~")
-            if (missing(gmpmControl)) {
-              gmpmControl <- x at gmpmControl
-            }
+            if (!missing(gmpmControl)) {
+              .setOpts(x, gmpmControl)
+              #gmpmControl <- x at gmpmControl
+            } else {}
+
+            x at nCores <- .calculateCores(x at gmpmControl)
+
             if (is.null(x at coef0)) {
               x at coef0 <- coef(fitOnce(x))
             }
 
-            maxruns <- gmpmControl$maxruns
-            report.interval <- gmpmControl$report.interval
-            outfile <- gmpmControl$outfile
+            maxruns <- x at gmpmControl[["maxruns"]]
+            report.interval <- x at gmpmControl[["report.interval"]]
+            outfile <- x at gmpmControl[["outfile"]]
+            
             elapsed <- rep(NA, maxruns)
 
             if (!is.null(outfile)) {
@@ -1370,31 +1380,64 @@
             if (x at nBetween > 0) {
               pbnames <- names(x at psec)[1:x at nBetween]
             } else {}
+
+            # one "sweep" is defined as a set of parallel runs
+            # on each of the requested processing cores.
+            #
+            # the following code sets up the calls for doing the
+            # fitting over multiple processors.
+            nsweeps <- floor(maxruns / x at nCores)
+            listsize <- rep(x at nCores, nsweeps)
+            if ((maxruns %% x at nCores) > 0) {
+              nsweeps <- nsweeps + 1
+              listsize <- c(listsize, maxruns %% x at nCores)
+            } else {}
+
+            if ("multicore" %in% installed.packages()) {
+              mycall <- "mclapply"
+            } else {
+              mycall <- "lapply"
+            }
             
-            for (i in 1:maxruns) {
+            x.list <- list()
+            for (i in 1:x at nCores) {
+              x.list[[i]] <- x
+            }
+            
+            for (i in 1:nsweeps) {
+              if (i == nsweeps) {
+                x.list <- x.list[1:listsize[i]]
+              } else {}
+              thisrow <- (i-1)*x at nCores+2
+              
               t1 <- proc.time()["elapsed"]
-
+              
               if (x at nBetween > 0) {
                 for (j in 1:x at nBetween) {
-                  tiv <- x at ivBetween[j]
-                  x2 <- permute(x,tiv)
-                  myFit <- fitOnce(x2)
-                  .storeFitResult(x, myFit, pbnames[j], i+1)
+                  x2.list <- do.call(mycall, args=list(X=x.list, FUN=permute, thisiv=x at ivBetween[j]))
+                  myFit.list <- do.call(mycall, args=list(X=x2.list, FUN=fitOnce))
+                  for (mm in 1:listsize[i]) {
+                    .storeFitResult(x, myFit.list[[mm]], pbnames[j], thisrow+mm-1)
+                  }
+                  #myFit <- fitOnce(x2)
                   #x at psec[[pbnames[j]]][i+1,] <- cfs
                 }
               } else {}
 
               if (x at nWithin > 0) {
-                x2 <- permute(x)
-                myFit <- fitOnce(x2)
-                .storeFitResult(x, myFit, "within", i+1)
+                x2.list <- do.call(mycall, args=list(X=x.list, FUN=permute))
+                #x2 <- permute(x)
+                myFit.list <- do.call(mycall, args=list(X=x2.list, FUN=fitOnce))
+                for (mm in 1:listsize[i]) {
+                  .storeFitResult(x, myFit.list[[mm]], "within", thisrow+mm-1)
+                }
                 #x at psec[["within"]][i+1,] <- cfs
                 #.storeFitResult(x, ff, i+1)
               } else {}
 
-              if (!is.null(outfile)) {
-                .writeFit(x, myFit, outfile, TRUE)
-              } else {}
+              #if (!is.null(outfile)) {
+              #  .writeFit(x, myFit, outfile, TRUE)
+              #} else {}
 
               t2 <- proc.time()["elapsed"]
               elapsed[i] <- t2-t1
@@ -1402,12 +1445,12 @@
               if (report.interval > 0) {
                 if ((i == maxruns) ||
                     ((i %% report.interval)==0)) {
-                  .reportProgress(x, i, maxruns, elapsed[1:i])
+                  .reportProgress(x, i, nsweeps, elapsed[1:i])
                 } else {}
               } else {}
             }
 
-            x at ncomp <- i
+            x at ncomp <- sum(listsize[1:i])
             x at ndigits <- ifelse(x at ncomp > 9, ceiling(log(x at ncomp+1,base=10)), 1)
 
             #print("~~~ leaving gmpmFit ~~~")            
@@ -1451,3 +1494,28 @@
               names(psec) <- names(x at psec)
               return(psec)
           })
+
+setMethod(".setOpts",
+          signature(x="GMPM"),
+          function(x, opts) {
+            nameObject <- deparse(substitute(x))            
+            
+            if (!is.null(opts[["nCores"]])) {
+              x at gmpmControl[["nCores"]] <- .calculateCores(opts)
+              x at nCores <- as.integer(x at gmpmControl[["nCores"]])
+            } else {}
+
+            if (!is.null(opts[["maxruns"]])) {
+              x at gmpmControl[["maxruns"]] <- opts[["maxruns"]]
+            } else {}
+
+            if (!is.null(opts[["report.interval"]])) {
+              x at gmpmControl[["report.interval"]] <- opts[["report.interval"]]
+            } else {}
+
+            if (!is.null(opts[["outfile"]])) {
+              x at gmpmControl[["outfile"]] <- opts[["outfile"]]
+            } else {}
+            assign(nameObject, x, envir=parent.frame())            
+            return(invisible(x))
+          })

Modified: pkg/R/helpers.R
===================================================================
--- pkg/R/helpers.R	2010-03-09 22:16:48 UTC (rev 16)
+++ pkg/R/helpers.R	2010-06-10 17:35:17 UTC (rev 17)
@@ -7,7 +7,7 @@
 }
 
 gmpmCreate <- function(formula, family, data=parent.frame(), ivars=c(),
-                gmpmControl=gmpmCtrl(), ...)
+                gmpmControl=list(), ...)
 {
   gmpmFormals <- setdiff(names(formals(gmpmCreate)), "...")
   
@@ -57,7 +57,8 @@
     }
   else {}
   
-  x at gmpmControl <- gmpmControl
+  x at gmpmControl <- .gmpmCtrl()
+  .setOpts(x, gmpmControl)
   
   .parseFormula(x, formula)
 
@@ -102,13 +103,13 @@
   x at coefTerms <- .getCoefTerms(x)  
 
   x at ivix <- .getIVix(x)
- 
+  
   return(x)
 }
   
 
 gmpm <- function(formula, family, data=parent.frame(), ivars=c(),
-                gmpmControl=gmpmCtrl(), ...)
+                gmpmControl=list(), ...)
 {
   mc <- match.call()
   mc[[1]] <- gmpmCreate
@@ -119,12 +120,16 @@
   # fit the object
   return(gmpmEstimate(x))
 }
-  
-gmpmCtrl <- function(maxruns = 999, report.interval=10,
-                        outfile = NULL)
+
+gmpmCtrl <- function(...) {
+  stop("gmpmCtrl function deprecated as of gmpm version 0.5-1.")
+}
+
+.gmpmCtrl <- function(maxruns = 999, report.interval=10,
+                        outfile = NULL, nCores="all")
 {
   return(list(maxruns=maxruns, report.interval=report.interval,
-              outfile=outfile))
+              outfile=outfile, nCores=nCores))
 }
 
 .getIVtypes <- function(x, ivars, munit)
@@ -194,3 +199,39 @@
   }
   return(sigvec)
 }
+
+.calculateCores <- function(gmpmControl = list()) {
+  nCores <- 1
+
+  if (is.null(gmpmControl[["nCores"]])) {
+    nCores <- 1
+  } else {
+    if ("multicore" %in% installed.packages()) {
+      library(multicore)
+      if (gmpmControl[["nCores"]] == "all") {
+        nCores = multicore:::detectCores()
+      } else {
+        if (gmpmControl[["nCores"]] == "all.but.one") {
+          nCores = multicore:::detectCores() - 1
+          if (nCores <= 0) {
+            warning("only one core available; 'all.but.one' option in gmpmControl was ignored")
+            nCores = 1
+          } else {}
+        } else {
+          if (is.numeric(gmpmControl[["nCores"]])) {
+            if (gmpmControl[["nCores"]] > multicore:::detectCores()) {
+              nCores = multicore:::detectCores()              
+              warning("more processing cores requested (", gmpmControl[["nCores"]], ") than available (", nCores, "); using only the first ", nCores)
+            } else {
+              nCores = as.integer(gmpmControl[["nCores"]])
+            }
+          } else {
+            stop("unrecognized argument for 'nCores' to gmpmControl")
+          }
+        }
+      }
+    } else {}
+  }
+  
+  return(nCores)
+}

Modified: pkg/man/GMPMSummary-class.Rd
===================================================================
--- pkg/man/GMPMSummary-class.Rd	2010-03-09 22:16:48 UTC (rev 16)
+++ pkg/man/GMPMSummary-class.Rd	2010-06-10 17:35:17 UTC (rev 17)
@@ -27,7 +27,7 @@
     \sQuote{showReg} below).}
     
     \item{\code{showReg}:}{Object of class \code{"logical"} whether or
-      not the regression summary is to be calculated and displayed.}
+      not the regression summary is to be calculated and displayed (default=TRUE).}
     
   }
 }

Modified: pkg/man/gmpm.Rd
===================================================================
--- pkg/man/gmpm.Rd	2010-03-09 22:16:48 UTC (rev 16)
+++ pkg/man/gmpm.Rd	2010-06-10 17:35:17 UTC (rev 17)
@@ -12,12 +12,11 @@
 
 \usage{
 
-gmpmCreate(formula, family, data=parent.frame(), ivars,
-gmpmControl=gmpmCtrl(), ...)
+gmpmCreate(formula, family, data=parent.frame(), ivars, gmpmControl=list(), ...)
 
 gmpmEstimate(x, gmpmControl)
 
-gmpm(formula, family, data = parent.frame(), ivars, gmpmControl = gmpmCtrl(), ...)
+gmpm(formula, family, data = parent.frame(), ivars, gmpmControl=list(), ...)
 
 }
 
@@ -44,17 +43,27 @@
   \sQuote{covariates} and will not be subject to randomization.}
 
   \item{gmpmControl}{A gmpmControl object used to control the fitting
-    function.  This is a list one or more of the following options:
+    function.  This is a list one or more of the following options.  Any
+    option not specified in the list will be given a default value.
 
     \itemize{
-      \item{maxruns}{The number of Monte Carlo runs to be performed.}
+      \item{maxruns}{The number of Monte Carlo runs to be performed
+      (default = 999).}
 
       \item{report.interval}{The run interval at which to intermittently
 	report, during the fitting process, the number of runs completed and
-	the estimated time remaining.}
-      }
-    }
+	the estimated time remaining (default = 10).}
 
+      \item{nCores}{The number of processing cores to use (default =
+	"all").  Using more processors will generally result in faster
+	performance.  This variable can either be the character
+	expression "all" (use all available cores), the character
+	expression "all.but.one" (leave one core out so that the system
+	is not bogged down) or an integer.  The multicore package must
+	be installed, otherwise the number of cores will default to
+	one.}
+    } }
+
     \item{x}{A GMPM object}
     
   \item{\dots}{Other arguments to be passed on to the underlying function

Deleted: pkg/man/gmpmCtrl.Rd
===================================================================
--- pkg/man/gmpmCtrl.Rd	2010-03-09 22:16:48 UTC (rev 16)
+++ pkg/man/gmpmCtrl.Rd	2010-06-10 17:35:17 UTC (rev 17)
@@ -1,19 +0,0 @@
-\name{gmpmCtrl}
-\alias{gmpmCtrl}
-\title{Return a gmpmCtrl Object with Default Settings}
-\description{
-This is a convenience function that returns a gmpmControl object for
-controlling Monte Carlo runs, with default settings.  (Will probably be
-deprecated in future versions of package \pkg{gmpm}.)
-}
-\usage{
-gmpmCtrl(maxruns = 999, report.interval = 10, outfile = NULL)
-}
-\arguments{
-  \item{maxruns}{total number of Monte Carlo runs}
-  \item{report.interval}{interval at which to report back to user}
-  \item{outfile}{(for future development)}
-}
-\value{
-  A list, with elements as indicated above.
-}

Modified: pkg/man/kb07.Rd
===================================================================
--- pkg/man/kb07.Rd	2010-03-09 22:16:48 UTC (rev 16)
+++ pkg/man/kb07.Rd	2010-06-10 17:35:17 UTC (rev 17)
@@ -54,7 +54,7 @@
 # you will need to increase the number of runs in the command below
 # (to, say, 999) to get sensible results
 
-kb07.binom.gmpm <- gmpmEstimate(kb07.binom.gmpm, gmpmCtrl(maxruns=19))
+kb07.binom.gmpm <- gmpmEstimate(kb07.binom.gmpm, list(maxruns=19))
 
 summary(kb07.binom.gmpm)
 
@@ -67,7 +67,7 @@
 
 # again, you'll need to increase the number of runs to get something
 # sensible
-kb07.mnom.gmpm <- gmpmEstimate(kb07.mnom.gmpm, gmpmCtrl(maxruns=19))
+kb07.mnom.gmpm <- gmpmEstimate(kb07.mnom.gmpm, list(maxruns=19))
 
 summary(kb07.mnom.gmpm)
 

Modified: pkg/man/mdTest.Rd
===================================================================
--- pkg/man/mdTest.Rd	2010-03-09 22:16:48 UTC (rev 16)
+++ pkg/man/mdTest.Rd	2010-06-10 17:35:17 UTC (rev 17)
@@ -38,7 +38,8 @@
 summary(aov(Y ~ A*B + Error(SubjID), df1))
 
 # randomization analysis; should increase maxruns to something like 999
-df1.gmpm <- gmpm(Y ~ A*B | SubjID, gaussian, df1, c("A","B"), gmpmCtrl(maxruns=99))
+df1.gmpm <- gmpm(Y ~ A*B | SubjID, gaussian, df1, c("A","B"),
+                 gmpmControl=list(maxruns=99))
 
 # retrieve the permutation matrix
 pmx <- getPermMx(df1.gmpm)[["within"]]



More information about the Gmpm-commits mailing list