added functionality to extrac dedicated gamma matrix combinations

R/readutils.R

@@ -5,38 +5,10 @@ readcmicor <- function(filename) {
   return(invisible(data))
 }
 
-readcmidatafiles <- function(path=".", basename="onlinemeas", skip=1, ind.vec=c(1,2,3,4,5),
-                             last.digits=4, verbose=FALSE, sym.vec) {
+readcmidatafiles <- function(path="./", basename="onlinemeas", skip=1,
+                             last.digits=4, verbose=FALSE) {
   ## read one file to determine parameters
   ofiles <- dir(path=path, pattern=paste(basename, "*", sep=""))
-  tmp <- read.table(paste(path, ofiles[1], sep=""), skip=skip)
-  nrObs <- max(tmp[,ind.vec[1]])
-  if(nrObs != length(unique(tmp[,ind.vec[1]]))) {
-    stop("data inconsistent, nrObs not equal to input data length\n")
-  }
-  nrStypes <- length(unique(tmp[,ind.vec[2]]))
-  Time <-  2*max(tmp[,ind.vec[3]])
-  Thalf <- max(tmp[,ind.vec[3]])+1
-  if(Thalf != length(unique(tmp[,ind.vec[3]]))) {
-    stop("data inconsistent, T not equal to input data length\n")
-  }
-  if(verbose) cat("nrObs=",nrObs, "nrStypes=",nrStypes, "T=", Time, "#files=", length(ofiles), "\n")
-
-  ## memory for correlators
-  cf <- array(0, dim=c(length(ofiles), nrObs*Thalf*nrStypes))
-
-  ## symmetrise or anti-symmetrise for given observable?
-  isym.vec <- rep(+1, times= nrObs*Thalf*nrStypes)
-  if(!missing(sym.vec)) {
-    if(length(sym.vec) != nrObs) {
-      stop("sym.vec was given, but had not the correct length")
-    }
-    for(i in c(1:nrObs)) {
-      if(!sym.vec) {
-        isym.vec[((i-1)*Thalf*nrStypes+1):((i)*Thalf*nrStypes)] <- -1
-      }
-    }
-  }
 
   ## sort input files using the last 4 characters of the filename
   e <- nchar(ofiles[1])
@@ -47,18 +19,55 @@ readcmidatafiles <- function(path=".", basename="onlinemeas", skip=1, ind.vec=c(
 
   ## read single files
   j <- 1
+  cmicor <- data.frame()
   for(i in ii) {
     if(verbose) {
       cat("Reading from file", ofiles[i], "\n")
     }
-    tmp <- read.table(paste(path, ofiles[i], sep=""), skip=skip)
-    ## take care of the zeroth and last times, which don't need to be averaged
-    tmp[(tmp[,ind.vec[3]]!=0 & (tmp[,ind.vec[3]]!=(Thalf-1))),] <-
-      tmp[(tmp[,ind.vec[3]]!=0 & (tmp[,ind.vec[3]]!=(Thalf-1))),]/2
-    cf[j,] <- (tmp[,ind.vec[4]] + isym.vec*tmp[,ind.vec[5]])
+    cmicor <- rbind(cmicor, read.table(paste(path, ofiles[i], sep=""), skip=skip))
     j <- j+1
   }
-  return(invisible(list(cf=cf, gaugeno=gaugeno[ii], Time=Time, nrStypes=nrStypes, nrObs=nrObs)))
+  attr(cmicor, "class") <- c("cmicor", "data.frame")  
+  return(invisible(cmicor))
+}
+
+extract.obs <- function(cmicor, vec.obs=c(1), ind.vec=c(1,2,3,4,5), sym.vec, verbose=FALSE) {
+  if(missing(cmicor)) {
+    stop("data missing in extract.obs\n")
+  }
+
+  nrObs <- max(cmicor[,ind.vec[1]])
+  if(nrObs != length(unique(cmicor[,ind.vec[1]]))) {
+    stop("extract.obs: data inconsistent, nrObs not equal to input data length\n")
+  }
+  nrStypes <- length(unique(cmicor[,ind.vec[2]]))
+  Time <-  2*max(cmicor[,ind.vec[3]])
+  Thalf <- max(cmicor[,ind.vec[3]])+1
+  if(Thalf != length(unique(cmicor[,ind.vec[3]]))) {
+    stop("extract.obs: data inconsistent, T not equal to input data length\n")
+  }
+  if(verbose) cat("nrObs=",nrObs, "nrStypes=",nrStypes, "T=", Time, "\n")
+
+  nrObs <- length(vec.obs)
+  data <- cmicor[cmicor[,ind.vec[1]] %in% vec.obs,]
+  data[(data[,ind.vec[3]]!=0 & (data[,ind.vec[3]]!=(Thalf-1))),ind.vec[c(4,5)]] <-
+      data[(data[,ind.vec[3]]!=0 & (data[,ind.vec[3]]!=(Thalf-1))),ind.vec[c(4,5)]]/2
+  ## symmetrise or anti-symmetrise for given observable?
+  isym.vec <- rep(+1, times= nrObs*Thalf*nrStypes)
+  if(!missing(sym.vec)) {
+    if(length(sym.vec) != nrObs) {
+      stop("sym.vec was given, but had not the correct length")
+    }
+    for(i in c(1:nrObs)) {
+      if(!sym.vec) {
+        isym.vec[((i-1)*Thalf*nrStypes+1):((i)*Thalf*nrStypes)] <- -1
+      }
+    }
+  }
+  cf <- array((data[,ind.vec[4]] + isym.vec*data[,ind.vec[5]]),
+              dim=c(nrObs*Thalf*nrStypes, length(data[,1])/(nrObs*Thalf*nrStypes)))
+
+  return(invisible(list(cf=cf, Time=Time, nrStypes=nrStypes, nrObs=nrObs)))
 }
 
 readhlcor <- function(filename) {

man/extract.obs.Rd

@@ -0,0 +1,60 @@
+\name{extract.obs}
+\title{Extract One or More Gamma Combinations from am CMI Correlator}
+\description{
+  Extracts one or more gamma matrix combinations (observables) from a
+  correlator stored in cmi format
+}
+\usage{
+extract.obs(cmicor, vec.obs=c(1), ind.vec=c(1,2,3,4,5), sym.vec, verbose=FALSE)
+}
+\arguments{
+  \item{vec.obs}{
+    vector containing the numbers of observables to be extracted.
+  }
+  \item{ind.vec}{
+    Index vector indexing the column numbers in cmicor to be used. The
+    first must be the observable index, the second the smearing type
+    index, the third the time, the fourth C(+t) and the fifth C(-t).
+  }
+  \item{verbose}{
+    Increases verbosity of the function.
+  }
+  \item{sym.vec}{
+    a vector of bools of length equal to the number of observables
+    indicating whether C(+t) and C(-t) should be added or subtracted. If
+    not given C(+t) and C(-t) will be averaged for all observables.
+  }
+}
+\value{
+  returns a list containing
+  \item{cf}{
+    array containing the correlation function with dimension number of
+    files times (nrObs*nrStypes*(T/2+1)). C(t) and C(-t) are averaged
+    according to \code{sym.vec}.
+  }
+  \item{Time}{
+    The time extend of the correlation functions.
+  }
+  \item{nrStypes}{
+    The number of smearing combinations.
+  }
+  \item{nrObs}{
+    The number of observables.
+  }
+}
+\details{
+  C(t) and C(-t) are averaged as indicated by \code{sym.vec}. This
+  object can be further processed by e.g. \code{block.ts}.
+}
+\references{
+}
+\seealso{
+  \code{\link{readcmicor}}, \code{\link{readcmidatafiles}},
+  \code{\link{block.ts}}
+}
+\examples{
+\dontrun{cmicor <- readcmidatafiles("outprc", skip=1)}
+\dontrun{cf <- extract.obs(cmicor, vec.obs=c(1,3))
+}
+\author{Carsten Urbach, \email{curbach@gmx.de}}
+\keyword{ts}

man/readcmidatafiles.Rd

@@ -5,8 +5,8 @@
   reads data from single files in Chris Michael format
 }
 \usage{
-readcmidatafiles(path=".", basename="onlinemeas", skip=1, ind.vec=c(1,2,3,4,5),
-                 last.digits=4, verbose=FALSE, sym.vec) 
+readcmidatafiles(path=".", basename="onlinemeas", skip=1,
+                 last.digits=4, verbose=FALSE) 
 }
 \arguments{
   \item{path}{
@@ -38,24 +38,8 @@ readcmidatafiles(path=".", basename="onlinemeas", skip=1, ind.vec=c(1,2,3,4,5),
   }
 }
 \value{
-  returns a list with
-  \item{cf}{
-    array containing the correlation function with dimension number of
-    files times (nrObs*nrStypes*(T/2+1)). C(t) and C(-t) are averaged
-    according to \code{sym.vec}.
-  }
-  \item{gaugeno}{
-    The ordered list of gauge configuration indices.
-  }
-  \item{Time}{
-    The time extend of the correlation functions.
-  }
-  \item{nrStypes}{
-    The number of smearing combinations.
-  }
-  \item{nrObs}{
-    The number of observables.
-  }
+  returns an object of class \code{cmicor}, which is an \code{rbind} of
+  all \code{data.frame}s read from the single files.
 }
 \details{
   This function reads data from single data files. The data in the files is assumed
@@ -66,7 +50,7 @@ readcmidatafiles(path=".", basename="onlinemeas", skip=1, ind.vec=c(1,2,3,4,5),
   time; 5) \eqn{c_2}{c2} correlator value at time value backward in time;
 }
 \seealso{
-  \code{\link{readcmicor}}
+  \code{\link{readcmicor}}, \code{\link{extract.obs}}
 }
 \examples{
 library(hadron)