Nothing
#' @title Container for raw correlation functions
#'
#' @description
#' This function \code{raw_cf()} creates containers for raw correlation functions
#' of class \code{raw_cf}. This class is particularly designed to deal with
#' complex and matrix-valued correlation functions emerging in statistical
#' mechanics and quantum field theory simulations.
#' Arithmetic operations are defined for this class and utility functions
#' such as \code{is.raw_cf} and \code{is_empty.raw_cf}.
#'
#' @family raw_cf constructors
#'
#' @return
#' An object of S3 class `raw_cf`.
#'
#' @export
raw_cf <- function () {
cf <- list()
class(cf) <- append(class(cf), 'raw_cf')
return (cf)
}
#' @title \code{raw_cf} metadata mixin constructor
#'
#' @param cf initial \link{raw_cf} object
#' @param nrObs Integer, number of different observables assembled in the data field of this container.
#' @param Time Integer, full time extent.
#' @param nts Integer, number of time separations actually stored in the data field.
#' @param nrStypes Integer, number of smearing types.
#' @param dim Integer vector of "internal" dimensions for matrix-valued correlation functions.
#' For a scalar correlation, this should be specified as \code{c(1,1)}.
#'
#' @family raw_cf constructors
#'
#' @return
#' An object of S3 class `raw_cf` with metadat mixing added.
#'
#' @export
raw_cf_meta <- function (cf = raw_cf(), nrObs = 1, Time = NA, nrStypes = 1, dim=c(1,1), nts = Time ) {
stopifnot(inherits(cf, 'raw_cf'))
cf$nrObs <- nrObs
cf$Time <- Time
cf$nts <- nts
cf$nrStypes <- nrStypes
cf$dim <- dim
class(cf) <- append(class(cf), 'raw_cf_meta')
return (cf)
}
#' @title Original data mixin constructor for `raw_cf`
#'
#' @param cf `raw_cf` object to extend.
#' @param data Numeric or complex array, original data for all observables and measurements.
#' This should have dimensions c(Nmeas,cf$Time*cf$nrObs*cf$nrStypes,cf$dim).
#' Having the internal dimensions innermost is not as efficient, but it allows
#' different transformations to be applied to different observables in the same
#' container more easily.
#'
#' @family raw_cf constructors
#'
#' @return
#' An object of S3 class `raw_cf` with original data mixin added.
#'
#' @export
raw_cf_data <- function (cf, data) {
stopifnot(inherits(cf, 'raw_cf'))
stopifnot(inherits(cf, 'raw_cf_meta'))
dims <- dim(data)
if( any( cf$dim != dims[3:length(dims)] ) ){
stop("Dimensions of 'data' and 'raw_cf$dim' container meta-data must match")
}
# raw correlation functions are expected to be complex-valued
if( mode(data) != "complex" ){
cf$data <- as.complex(data)
} else {
cf$data <- data
}
class(cf) <- append(class(cf), 'raw_cf_data')
return (cf)
}
#' @title Extract a particular internal component of a 'raw_cf' into a 'cf'
#' @param x 'raw_cf' container with 'raw_cf_data' and 'raw_cf_meta'
#' @param component Integer vector of the same length as the internal dimension
#' of the 'raw_cf' specifying which component should be extracted.
#' @return 'cf' object
#' @export
raw_cf_to_cf <- function(x, component){
stopifnot(inherits(x, 'raw_cf_meta'))
stopifnot(inherits(x, 'raw_cf_data'))
if( length(component) != length(x$dim) ){
stop("The component specifier must be of the same length as the vector of internal dimensions of the 'raw_cf'!")
}
if( x$Time != x$nts ){
stop("The 'nts' dimension must be the same as the 'Time' dimension to convert to 'cf'")
}
idcs <- idx_matrix.raw_cf(x, component)
dims <- dim(x$data)
subset <- array(x$data[idcs], dim=c(dims[1],dims[2]))
cf <- cf_meta(Time = x$Time,
nrObs = x$nrObs,
nrStypes = x$nrStypes,
symmetrised = FALSE)
cf <- cf_orig(cf,
cf = Re(subset),
icf = Im(subset))
return(cf)
}
#' @title Gamma method analysis on all time-slices in a 'raw_cf' object
#'
#' @param cf Correlation function container of class 'raw_cf'
#' @return The return value is a list with elements
#' \describe{
#' \item{\code{value}}{central value}
#' \item{\code{dvalue}}{statistical error}
#' \item{\code{ddvalue}}{error of the statistical error}
#' \item{\code{tauint}}{auto-correlation time estimate}
#' \item{\code{dtauint}}{error of auto-correlation time estimate}
#' }
#' Each of these is in turn an array of dimension \code{ c( cf$nts, cf$dim ) } and
#' hance lacks the first dimension index compared for \code{cf$data}.
#'
#' @export
uwerr.raw_cf <- function(cf){
stopifnot(inherits(cf, 'raw_cf_data'))
dims <- dim(cf$data)
# prepare some output tensors which are reduced over the first
# dimension (the measurements)
na_array <- array(as.numeric(NA), dim=dims[2:length(dims)])
value <- list(real = na_array, imag = na_array)
dvalue <- list(real = na_array, imag = na_array)
ddvalue <- list(real = na_array, imag = na_array)
tauint <- list(real = na_array, imag = na_array)
dtauint <- list(real = na_array, imag = na_array)
# generate an index matrix for this reduced tensor
args <- list()
for( d in dim(na_array) ){
args[[length(args)+1]] <- 1:d
}
idx_matrix_uwerr <- as.matrix(do.call(expand.grid, args))
# the measurements are the fastest running dimension so we can easily
# perform the Gamma method analysis by stepping through the index matrix
# get the index matrix for the correlator tensor
idcs <- idx_matrix.raw_cf(cf)
stepping <- dims[1]
nsteps <- as.integer( nrow(idcs) / stepping )
for( step in 0:(nsteps-1) ){
istart <- step*stepping + 1
iend <- istart + stepping - 1
# this is really annoying. When we subset a single row of the index matrix
# it's turned into a column vector
out_idcs <- t(idx_matrix_uwerr[step+1,])
tseries <- cf$data[ idcs[istart:iend,] ]
uw_real <- try(uwerrprimary(data = Re(tseries)), silent=TRUE)
uw_imag <- try(uwerrprimary(data = Im(tseries)), silent=TRUE)
if(!inherits(uw_real,'try-error')){
value$real[out_idcs] <- uw_real$value
dvalue$real[out_idcs] <- uw_real$dvalue
ddvalue$real[out_idcs] <- uw_real$ddvalue
tauint$real[out_idcs] <- uw_real$tauint
dtauint$real[out_idcs] <- uw_real$dtauint
}
if(!inherits(uw_imag, 'try-error')){
value$imag[out_idcs] <- uw_imag$value
dvalue$imag[out_idcs] <- uw_imag$dvalue
ddvalue$imag[out_idcs] <- uw_imag$ddvalue
tauint$imag[out_idcs] <- uw_imag$tauint
dtauint$imag[out_idcs] <- uw_imag$dtauint
}
}
return(list(value = value,
dvalue = dvalue,
ddvalue = ddvalue,
tauint = tauint,
dtauint = dtauint,
idx_matrix = idx_matrix_uwerr))
}
#' @title Block average correlation function data
#' @description Block `block_length` sequential measurements of the correlation
#' function together. This occurs, for example, when multiple
#' stochastic noise vectors are used per measurement or multiple
#' source locations. Alternatively, it can also be used to
#' account for auto-correlations in the data. If the total number
#' of measurements is not divisible by `block_length`, the last
#' measurements are discarded.
#'
#' @param cf `raw_cf` object
#' @param block_length Integer, number of successive measurements to average over.
#' @return cf `raw_cf` object with the data member reduced in its first dimension
#' by a factor of `block_length` and restricted (at the end) to the number
#' of measurements divisible by `block_length`.
#'
#' @export
block.raw_cf <- function(cf, block_length){
stopifnot(inherits(cf, 'raw_cf'))
stopifnot(inherits(cf, 'raw_cf_data'))
if( block_length == 1 ){
return(cf)
}
idcs <- idx_matrix.raw_cf(cf)
dims <- dim(cf$data)
nblocks <- as.integer( dims[1] / block_length )
nelem <- as.integer( nrow(idcs) / dims[1] )
# array for blocked data as well as an index matrix for this array
out <- array(as.complex(NA),
dim = c(nblocks, cf$nts*cf$nrObs*cf$nrStypes, cf$dim))
args <- list()
for( d in dim(out) ){
args[[length(args)+1]] <- 1:d
}
idcs_out <- as.matrix(do.call(expand.grid, args))
iout <- 1
for( elem in 0:(nelem-1) ){
for( block in 0:(nblocks-1) ){
istart <- elem*dims[1] + block*block_length + 1
iend <- istart + block_length - 1
# select block of measurements
idcs_in <- idcs[istart:iend,]
# place block average into element of 'out' array
# be mindful of conversion to column vector
idx_out <- t(idcs_out[iout,])
out[idx_out] <- mean( cf$data[ idcs_in ] )
iout <- iout+1
}
}
# replace the data member by the blocked data
cf$data <- out
return(cf)
}
#' @title Extend statistics of an existing `raw_cf` container
#'
#' @param cf1 `raw_cf` container with or without 'data' and 'meta' mixins
#' @param cf2 `raw_cf` container with or without 'data' and 'meta' mixins
#'
#' @details
#' When either of `cf1` or `cf2` does not contain any data,
#' the other object is returned. (allows empty `raw_cf` to be extended).
#' If the dimensions (except for the measurements) of the `data` fields of the
#' two containers match, they are concatenated along the measurement dimension.
#'
#' @return
#' An object of S3 class `raw_cf` identical to the input object but with extended statistics.
#'
#' @export
addStat.raw_cf <- function(cf1, cf2) {
stopifnot(inherits(cf1, 'raw_cf'))
stopifnot(inherits(cf2, 'raw_cf'))
if (is_empty.raw_cf(cf1)) {
return (invisible(cf2))
}
if (is_empty.raw_cf(cf2)) {
return (invisible(cf1))
}
stopifnot(inherits(cf1, 'raw_cf_meta'))
stopifnot(inherits(cf2, 'raw_cf_meta'))
stopifnot(cf1$Time == cf2$Time)
stopifnot(cf1$nrObs == cf2$nrObs )
stopifnot(cf1$nrStypes == cf2$nrStypes)
stopifnot(inherits(cf1, 'raw_cf_data'))
stopifnot(inherits(cf2, 'raw_cf_data'))
stopifnot( all(dim(cf1$data)[2:length(dim(cf2$data))] == dim(cf2$data)[2:length(dim(cf2$data))] ) )
cf <- cf1
cf$data <- abind::abind(cf1$data, cf2$data, along=1)
return (invisible(cf))
}
#' @title add two `raw_cf` objects
#' @param cf1 first 'raw_cf' container with data and meta-data
#' @param cf2 second 'raw_cf' container with data and meta-data
#' @param a Numeric or complex, scaling factor applied to \code{cf1}.
#' @param b Numeric or complex, scaling factor applied to \code{cf2}.
#' @return \code{a*cf1$data + b*cf2$data}
#'
#' @export
add.raw_cf <- function(cf1, cf2, a=1.0, b=1.0) {
stopifnot(inherits(cf1, 'raw_cf'))
stopifnot(inherits(cf2, 'raw_cf'))
stopifnot(inherits(cf1, 'raw_cf_data'))
stopifnot(inherits(cf2, 'raw_cf_data'))
stopifnot(inherits(cf1, 'raw_cf_meta'))
stopifnot(inherits(cf2, 'raw_cf_meta'))
stopifnot(all(dim(cf1$data) == dim(cf2$data)))
stopifnot(cf1$Time == cf2$Time)
stopifnot(cf1$nts == cf2$nts)
cf <- cf1
cf$data <- a*cf1$data + b*cf2$data
return(cf)
}
#' @title add two `raw_cf` objects
#' @param cf1 first 'raw_cf' container to be added
#' @param cf2 second 'raw_cf' container to be added
#' @return `raw_cf` object with \code{cf$data == cf1$data + cf2$data}
#'
#' @export
'+.raw_cf' <- function (cf1, cf2) {
add.raw_cf(cf1, cf2, a = 1.0, b = 1.0)
}
#' @title add two `raw_cf` objects
#' @param cf1 first 'raw_cf' container to be subtracted
#' @param cf2 second 'raw_cf' container to be subtracted
#' @return `raw_cf` object with \code{cf$data == cf1$data - cf2$data}
#'
#' @export
'-.raw_cf' <- function(cf1, cf2) {
add.raw_cf(cf1, cf2, a = 1.0, b = -1.0)
}
#' @title divide two `raw_cf` objects
#' @param cf1 'raw_cf' container with data and meta-data to be the dividend
#' @param cf2 'raw_cf' container with data and meta-data to be the divisor
#' @return `raw_cf` object with \code{cf$data == cf1$data / cf2$data}
#'
#' @export
'/.raw_cf' <- function(cf1, cf2) {
stopifnot(inherits(cf1, 'raw_cf_meta'))
stopifnot(inherits(cf2, 'raw_cf_meta'))
stopifnot(inherits(cf1, 'raw_cf_data'))
stopifnot(inherits(cf2, 'raw_cf_data'))
stopifnot(all(dim(cf1$data) == dim(cf2$data)))
stopifnot(cf1$Time == cf2$Time)
stopifnot(cf1$nts == cf2$nts)
cf <- cf1
cf$data <- cf1$data / cf2$data
return (cf)
}
#' @title multiply two `raw_cf` objects
#' @param cf1 first 'raw_cf' container with data and meta-data to be multiplied
#' @param cf2 second 'raw_cf' container with data and meta-data to be multiplied
#' @return `raw_cf` object with \code{cf$data == cf1$data * cf2$data}
#'
#' @export
'*.raw_cf' <- function(cf1, cf2) {
stopifnot(inherits(cf1, 'raw_cf_meta'))
stopifnot(inherits(cf2, 'raw_cf_meta'))
stopifnot(inherits(cf1, 'raw_cf_data'))
stopifnot(inherits(cf2, 'raw_cf_data'))
stopifnot(all(dim(cf1$data) == dim(cf2$data)))
stopifnot(cf1$Time == cf2$Time)
stopifnot(cf1$nts == cf2$nts)
cf <- cf1
cf$data <- cf1$data * cf2$data
return(cf)
}
#' @title Take the complex conjugate of a `raw_cf` object
#' @param cf `raw_cf` cotnainer with data
#' @return `raw_cf`
conj_raw_cf <- function(cf){
stopifnot(inherits(cf, 'raw_cf_data'))
cf$data <- Conj(cf$data)
return(cf)
}
#' @title scale `raw_cf` data
#' @param cf 'raw_cf' container with data to be scaled by the factor \code{a}
#' @param a Numeric or complex scaling factor, although it could also be
#' an array of dimensions compatible with \code{cf$data}
#' @return `raw_cf` object with \code{res$data == a*cf$data}
#'
#' @export
mul.raw_cf <- function(cf, a=1.) {
stopifnot(inherits(cf, 'raw_cf_data'))
stopifnot(is.numeric(a) | is.complex(a))
cf$data <- a*cf$data
return (cf)
}
#' @title check if an object is of class `raw_cf`
#' @param x object to be checked
#'
#' @return
#' Returns TRUE if `x` is an object of class `raw_cf`, FALSE otherwise.
#'
#' @export
is.raw_cf <- function(x){
inherits(x, "raw_cf")
}
#' @title check if an obect is of class `raw_cf` and empty otherwise
#' @param x object to be checked
#'
#' @return
#' Returns TRUE if `x` is an empty object of class `raw_cf`, FALSE otherwise.
#'
#' @export
is_empty.raw_cf <- function(x){
.raw_cf <- x
setequal(class(.raw_cf), class(raw_cf())) & is.null(names(.raw_cf))
}
#' @title Concatenate `raw_cf` correlation function objects
#'
#' @param ... Zero or multiple objects of type `raw_cf`.
#'
#' @return
#' Returns an object of S3 class `raw_cf`, the concatenation of the
#' input objects.
#'
#' @export
c.raw_cf <- function (...) {
rval <- Reduce(concat.raw_cf, list(...), raw_cf())
return (rval)
}
#' @title Concatenate two `raw_cf` correlation function objects
#' @description The data of the \code{left} and \code{right} objects
#' is concatenated along the second array dimension
#' such that the output contains the tensor slices of \code{right}
#' after the slices of \code{left}
#' @param left `raw_cf` object to be concatenated with \code{right}
#' @param right `raw_cf` object to be concatenated with \code{left}
#'
#' @return
#' Returns an object of S3 class `raw_cf`, the concatenation of the
#' two input objects.
#'
#' @export
concat.raw_cf <- function (left, right) {
stopifnot(inherits(left, 'raw_cf'))
stopifnot(inherits(right, 'raw_cf'))
# In case that one of them does not contain data, the other one is the
# result. This satisfies the neutral element axiom of a monoid.
if (is_empty.raw_cf(left)) {
return (right)
} else if (is_empty.raw_cf(right)) {
return (left)
} else if ( is_empty.raw_cf(left) & is_empty.raw_cf(right) ) {
stop("concat.raw_cf: both `raw_cf` objects empty, cannot concatenate!")
}
stopifnot(inherits(left, 'raw_cf_meta'))
stopifnot(inherits(right, 'raw_cf_meta'))
# At this point both `raw_cf` objects given here have original data, therefore we
# need to concatenate them.
# A few checks for compatability.
stopifnot(left$Time == right$Time)
stopifnot(left$nts == right$nts)
stopifnot(nrow(left$data) == nrow(right$data))
stopifnot(left$nrStypes == right$nrStypes)
stopifnot(all(left$dim == right$dim))
rval <- raw_cf_meta(nrObs = left$nrObs + right$nrObs,
Time = left$Time,
nts = left$nts,
dim = left$im,
nrStypes = left$nrStypes)
rval <- raw_cf_data(cf = rval,
data = abind::abind(left$data, right$data, along=2))
return (rval)
}
#' @title extract data from 'raw_cf' in format convenient to plot
#' @description When dealing with with tensorial `raw_cf` objects
#' pre-processing and reshaping is always required to
#' prepare the data for plotting (or similar). This function
#' conveniently prepares a named list of prepared data.
#' The list elements are themselves lists which contain
#' `val` and `dval` members with the central value and error
#' of the element in question. These are in turn
#' arrays of dimension \code{ c( cf$nts, cf$dim ) } and thus
#' lack the first index compared to \code{cf$data}.
#' @param cf `raw_cf` object with meta-data and data.
#' @param reim String, one of 'real', 'imag' or 'both'. Specifies
#' whether the real and/or imaginary parts should be extracted.
#' @param tauint Boolean, specifies if the tensor of auto-correlation times
#' and corresponding errors should be extracted.
#' @param relerr Boolean, specifies if the return value should also include
#' estimates of the relative error and its error.
#' @return List of up to six named elements (depending on what was passed for
#' \code{reim}, \code{tauint}, \code{relerr}) containing the central
#' values and errors of the real and/or imaginary part of `cf$data`
#' as well as the corresponding arrays of auto-correlation times and
#' relative errors. The list elements come in the order
#' \code{real}, \code{imag}, \code{relerr_real}, \code{relerr_imag},
#' \code{tauint_real}, \code{tauint_imag} if \code{reim} is `both` and
#' \code{tauint} and \code{relerr} are \code{TRUE}. The \code{val} and
#' \code{dval} members of these list elements are arrays of dimension
#' \code{ c( cf$nts, cf$dim ) } and thus lack the first index compared
#' to \code{cf$data}.
#'
#' @export
get_plotdata_raw_cf <- function(cf,
reim,
tauint,
relerr)
{
stopifnot(inherits(cf, 'raw_cf_data'))
stopifnot(inherits(cf, 'raw_cf_meta'))
data <- uwerr.raw_cf(cf)
if( !(reim %in% c('real','imag','both') ) ){
stop("'reim' must be one of 'real', 'imag' or 'both'")
}
data <- uwerr.raw_cf(cf)
# construct a return value
plotdata <- list()
# loop over the internal dimensions
for( idim in 0:(prod(cf$dim)-1) ){
istart <- idim*cf$nrStypes*cf$nrObs*cf$nts + 1
iend <- istart + cf$nrStypes*cf$nrObs*cf$nts - 1
# subset of the index matrix for the internal dimension given by the
iselect <- data$idx_matrix[istart:iend,]
lidx <- length(plotdata)+1
plotdata[[lidx]] <- list()
if( reim %in% c('real','both') ){
# ensure that we append to the end in the right order
plotdata[[lidx]] <- append(plotdata[[lidx]],
list(real =
list(val = data$value$real[iselect],
dval = data$dvalue$real[iselect])
)
)
}
if( reim %in% c('imag','both') ){
plotdata[[lidx]] <- append(plotdata[[lidx]],
list(imag =
list(val = data$value$imag[iselect],
dval = data$dvalue$imag[iselect])
)
)
}
if( relerr ){
if( reim %in% c('real','both') ){
plotdata[[lidx]] <- append(plotdata[[lidx]],
list(relerr_real =
list(val = abs(data$dvalue$real[iselect] / data$value$real[iselect]),
dval = abs(data$ddvalue$real[iselect] / data$value$real[iselect]))
)
)
}
if( reim %in% c('imag','both') ){
plotdata[[lidx]] <- append(plotdata[[lidx]],
list(relerr_imag =
list(val = abs(data$dvalue$imag[iselect] / data$value$imag[iselect]),
dval = abs(data$ddvalue$imag[iselect] / data$value$imag[iselect]))
)
)
}
}
if( tauint ){
if( reim %in% c('real','both') ){
plotdata[[lidx]] <- append(plotdata[[lidx]],
list(tauint_real =
list(val = data$tauint$real[iselect],
dval = data$dtauint$real[iselect])
)
)
}
if( reim %in% c('imag', 'both') ){
plotdata[[lidx]] <- append(plotdata[[lidx]],
list(tauint_imag =
list(val = data$tauint$imag[iselect],
dval = data$dtauint$imag[iselect])
)
)
}
}
}
class(plotdata) <- append(class(plotdata), "plotdata_raw_cf")
return(plotdata)
}
#' @title plot all correlators in `raw_cf` object
#' @param x Object of class `raw_cf` with data and meta-data.
#' @param reim Character vector, may contain 'real', 'imag' or 'both'. Determines
#' whether the real and/or imaginary parts of the correlation funtions
#' should be plotted.
#' @param reim_same Boolean, determines whether the real and imaginary parts, if both
#' are to be plotted, will be plotted in the same plot.
#' @param ... Further parameters passed to \link{plotwitherror}.
#' @return Invisibly returns the plotdata, see \link{get_plotdata_raw_cf}.
#'
#' @export
plot.raw_cf <- function(x,
...,
reim = 'real',
reim_same = FALSE)
{
cf <- x
stopifnot(inherits(cf, 'raw_cf_data'))
stopifnot(inherits(cf, 'raw_cf_meta'))
if( !(reim %in% c('real','imag','both') ) ){
stop("'reim' must be one of 'real', 'imag' or 'both'")
}
if( reim_same & !(reim == 'both') ){
stop("'reim_same' can only be true if 'reim' is 'both'")
}
plotdata <- get_plotdata_raw_cf(cf, reim = reim, tauint = FALSE, relerr = FALSE)
for( lidx in 1:length(plotdata) ){
if( reim %in% c('real', 'both') & reim_same == FALSE){
plotwitherror(x = 0:(cf$nts-1),
y = plotdata[[lidx]]$real$val,
dy = plotdata[[lidx]]$real$dval,
...)
}
if( reim %in% c('imag', 'both') & reim_same == FALSE){
plotwitherror(x = 0:(cf$nts-1),
y = plotdata[[lidx]]$imag$val,
dy = plotdata[[lidx]]$imag$dval,
...)
}
if( reim_same ){
plotwitherror(x = rep(0:(cf$nts-1), times=2),
y = c(plotdata[[lidx]]$real$val, plotdata[[lidx]]$imag$val),
dy = c(plotdata[[lidx]]$real$dval, plotdata[[lidx]]$imag$dval),
...)
}
}
return(invisible(plotdata))
}
#' @title create convenient overview plots for a `raw_cf` object
#' @param cf 'raw_cf' container with data and meta-data
#' @param grid Optional, integer vector which satisfies
#' \code{prod(grid) == prod(cf$dim)}. This is passed to \code{par} via
#' \code{par(mfrow=grid)} to produce a grid of plots as defined by the
#' components of \code{grid}.
#' @param reim Vector of strings, one of 'real', 'imag' or 'both'. Specified whether
#' the real or imaginary parts (or both) should be plotted.
#' @param reim_same Boolean, whether real and imaginary parts should be plotted
#' on the same plot. If \code{TRUE}, then \code{reim} must
#' be 'both'. If this is given, the imaginary part as well as its
#' relative error and per-time-slice integrated autocorreation times
# are plotted in red.
#' @param relerr Boolean, whether a plot of the relative error per time slice
#' should be added.
#' @param tauint Boolean, whether a plot of the integrated auto-correlation time
#' on each time slice should be added.
#' @param value_logplot Boolean, whether the plot of the correlator should be
#' on a logarithmic vertical axis. (does not affect \code{tauint}
#' and \code{relerr}).
#' @param value_factor Numeric, either of length '1' or as long as the number of
#' correlation functions in \code{cf}. The data will be scaled
#' by this factor before plotting.
#' @param title Character vector, will be passed as the \code{main} argument to
#' \link{plotwitherror} which in turn passes it to \link{plot}. Can
#' be either of length '1' or \code{prod(cf$dim)}
#'
#' @return
#' No return value, only plots are generated.
#'
#' @export
overview_plot_raw_cf <- function(cf,
grid,
reim = 'real',
reim_same = FALSE,
relerr = FALSE,
tauint = FALSE,
value_logplot = TRUE,
value_factor = c(1),
title = '')
{
stopifnot(inherits(cf, 'raw_cf_data'))
stopifnot(inherits(cf, 'raw_cf_meta'))
if( !(reim %in% c('real','imag','both') ) ){
stop("'reim' must be one of 'real', 'imag' or 'both'")
}
if( reim_same & !(reim == 'both') ){
stop("'reim_same' can only be true if 'reim' is 'both'")
}
if(!missing(grid)){
if( prod(grid) != prod(cf$dim) ){
stop("'prod(grid)' must be equal to 'prod(cf$dim)'")
}
par_save <- par()
on.exit(par(par_save))
par(mfrow=grid)
}
# produce a vector of plot symbols such that we have one
# symbol per observable, all smearing types will be plotted
# with the same symbol and we wrap around when we run out of symbols
pch <- expand.grid(rep(1:cf$nts, times=cf$nrStypes), c(0:6, 15:18), KEEP.OUT.ATTRS=FALSE )
pch <- rep(pch[,2], length.out=cf$nts*cf$nrObs*cf$nrStypes)
ylabs <- list()
ylabs[["real"]] <- "Re[ C(t) ]"
ylabs[["imag"]] <- "Im[ C(t) ]"
ylabs[["real_imag"]] <- "Re:Im[ C(t) ]"
ylabs[["tauint_real"]] <- "tau_int{ Re[ C(t) ] }"
ylabs[["tauint_imag"]] <- "tau_int{ Im[ C(t) ] }"
ylabs[["tauint_real_imag"]] <- "tau_int{ Re:Im[ C(t) ] }"
ylabs[["relerr_real"]] <- "Re[ dC(t) ] / Re[ C(t) ]"
ylabs[["relerr_imag"]] <- "Im[ dC(t) ] / Im[ C(t) ]"
ylabs[["relerr_real_imag"]] <- "Re:Im[ dC(t) ] / Re:Im[ C(t) ]"
tmax <- cf$nts-1
plotdata <- get_plotdata_raw_cf(cf, reim=reim, relerr=relerr, tauint=tauint)
if( length(value_factor) != 1 && length(value_factor) != prod(cf$dim) ){
message <- sprintf(paste("If 'value_factor' is of length > 1,",
"it must be of the same length at the number of correlators",
"to be plotted [prod(cf$dim)] !.",
"In this case, length(value_factor)=%d, prod(cf$dim)=%d"),
length(value_factor), prod(cf$dim) )
stop(message)
}
# if only a single factor has been provided, replicate it as many times as there
# are correlators to be ploted
if( length(value_factor) == 1 ){
value_factor <- rep(value_factor, times = prod(cf$dim) )
}
if( length(title) != 1 && length(title) != prod(cf$dim) ){
message <- sprintf(paste("If 'title' is of length > 1,",
"it must be of the same length at the number of correlators",
"to be plotted [prod(cf$dim)] !.",
"In this case, length(title)=%d, prod(cf$dim)=%d"),
length(title), prod(cf$dim) )
stop(message)
}
if( length(title) == 1 ){
title <- rep(title, times=prod(cf$dim))
}
ts <- rep(0:tmax, times=cf$nrStypes*cf$nrObs)
if( reim_same ){
ts <- rep(ts,times=2)
}
step <- 1
if( reim_same ) step <- 2
for( lidx in 1:length(plotdata) ){
onames <- names(plotdata[[lidx]])
for( oidx in seq(1,length(onames),step) ){
args <- list()
args$x <- ts
args$xlab <- "t/a"
args$ylab <- ylabs[[onames[oidx]]]
args$main <- title[lidx]
args$pch <- pch
if( (onames[oidx] == "real" | onames[oidx] == "imag") & value_logplot ){
args$log <- 'y'
}
if( reim_same ){
args$y <- c(plotdata[[lidx]][[ onames[oidx] ]]$val, plotdata[[lidx]][[ onames[oidx+1] ]]$val)
args$dy <- c(plotdata[[lidx]][[ onames[oidx] ]]$dval, plotdata[[lidx]][[ onames[oidx+1] ]]$dval)
args$col <- c(rep("black", cf$nrStypes*cf$nrObs*cf$nts),
rep("red", cf$nrStypes*cf$nrObs*cf$nts))
args$ylab <- ylabs[[sprintf("%s_imag",onames[oidx])]]
} else {
args$y <- plotdata[[lidx]][[ onames[oidx] ]]$val
args$dy <- plotdata[[lidx]][[ onames[oidx] ]]$dval
}
# apply the scaling factor to teh expectation value of the observable
if( (onames[oidx] == "real" | onames[oidx] == "imag") ){
args$y <- value_factor[lidx]*args$y
args$dy <- value_factor[lidx]*args$dy
}
# when plotting relative errors, it is not useful to plot errors larger than
# 100%
if( grepl("relerr",onames[oidx]) ){
if( max(args$y, na.rm=TRUE) > 1.0 ){
args$ylim <- c(0,1)
}
}
do.call(plotwitherror, args)
}
}
}
#' @title shift a \code{raw_cf} correlation function by 'places' time-slices
#' @param cf \code{raw_cf} container
#' @param places Integer (possibly a vector), number of time slices that the correlation function
#' should be shifted by. Can be positive or negative. This can either
#' be a single value such that a shift by this many time slices will be
#' applied to every measurement or it can be a vector of values of the
#' same length as the number of measurements in \code{cf}. In that case,
#' a different shift will be applied to each measurement. This is useful
#' if it is important to preserve the absolute time coordinates of a
#' correlation function until some time-dependent transformations
#' have been applied.
#' @details
#' The correlation funtion \eqn{C(t)} is shifted in time to produce:
#' \deqn{C'(t) = C(t+places)}
#' using periodic boundary conditions in time.
#'
#' @return
#' Returns an object of class `raw_cf`, shifted compared to the input object.
#'
#' @export
shift.raw_cf <- function(cf, places) {
stopifnot(inherits(cf, 'raw_cf_meta'))
stopifnot(inherits(cf, 'raw_cf_data'))
if( (length(places) != 1) & (length(places) != dim(cf$data)[1]) ){
stop("'places' should be either of length '1' or of a length equalling the number of measurements in 'cf'")
}
if(length(places) == 1){
if(places == 0){
return(cf)
}
}
dims <- dim(cf$data)
if( length(places) == 1 ){
step <- dims[1]
} else {
step <- 1
}
# if places is of length '1', we do all measurements in one step (see index matrix construction below)
# else we do each measurement separately
for( imeas in seq(from = 1, to = dims[1], by = step) ){
# make sure that we don't do any work when places[
if( places[imeas] == 0 ) next
for( oidx in 0:(cf$nrObs-1) ){
for( sidx in 0:(cf$nrStypes-1) ){
# the 'time' indices of the observables and smearing types
istart <- cf$Time*cf$nrStypes*oidx + cf$Time*sidx + 1
iend <- istart + cf$Time - 1
# construct an argument list for do.call below
args <- list()
args[[1]] <- imeas:(imeas+step-1)
args[[2]] <- istart:iend
for( d in cf$dim ){
args[[length(args)+1]] <- 1:d
}
# construct the tensor index set for the output
out_dof <- as.matrix(do.call(expand.grid, args))
if( places[imeas] < 0 ){
ishift <- c( (iend - abs(places[imeas]) + 1):iend,
(istart:(iend-abs(places[imeas]))) )
} else {
ishift <- c( (istart+places[imeas]):iend,
istart:(istart+places[imeas]-1))
}
args[[2]] <- ishift
# construct the tensor index set for the input
in_dof <- as.matrix(do.call(expand.grid, args))
# shift the correlator tensor
cf$data[out_dof] <- cf$data[in_dof]
}
}
}
return(cf)
}
#' @title Construct the tensor index set for the entire raw correlator
#' @param cf 'raw_cf' container with data and meta-data
#' @param component Integer vector. Optional argument to obtain a subset of the
#' index matrix to access a particular element of the interior
#' dimensions. Must of the the same length as cf$dim.
#'
#' @return
#' An object of type matrix is returned containing the tensor index set.
#'
#' @export
idx_matrix.raw_cf <- function(cf, component){
stopifnot(inherits(cf, 'raw_cf_meta'))
stopifnot(inherits(cf, 'raw_cf_data'))
dims <- dim(cf$data)
args <- list()
for( d in dims ){
args[[length(args)+1]] <- 1:d
}
if(!missing(component)){
if( length(component) != length(cf$dim) ){
stop("'component' has to be of length length(cf$dim)!")
}
cidx <- 1
for( didx in 3:length(dims) ){
args[[didx]] <- component[cidx]
cidx <- cidx + 1
}
}
as.matrix(do.call(expand.grid, args))
}
#' @title Construct tensor index set for the internal degrees of freedom
#' @param cf `raw_cf` container
#'
#' @return
#' Returns a matrix containing the above mentioned index set.
#'
#' @export
int_idx_matrix.raw_cf <- function(cf){
stopifnot(inherits(cf, 'raw_cf_meta'))
args <- list()
for( d in cf$dim ){
args[[length(args)+1]] <- 1:d
}
as.matrix(do.call(expand.grid, args))
}
#' @title Print summary of data contained in `raw_cf` container
#' @param object `raw_cf` container with data and meta-data
#' @param ... ignored
#' @param statistics Boolean, return central value and error
#' for all components of the 'raw_cf'. This can
#' be slow so the default is \code{FALSE}.
#'
#' @return
#' The summary is returned invisibly in form of a data frame.
#'
#' @export
summary.raw_cf <- function(object, ..., statistics = FALSE) {
cf <- object
stopifnot(inherits(cf, 'raw_cf_meta'))
stopifnot(inherits(cf, 'raw_cf_data'))
cat("T = ", cf$Time, "\n")
cat("observations = ", dim(cf$data)[1], "\n")
cat("Nr Stypes = ", cf$nrStypes, "\n")
cat("Nr Obs = ", cf$nrObs, "\n")
cat("dim = ", cf$dim, "\n")
if( statistics ){
uw <- uwerr.raw_cf(cf)
idcs <- uw$idx_matrix
out <- NULL
for( i in 1:nrow(idcs) ){
# annoying: subsetting a single row of the index matrix turns the
# result into a column vector, need to transpose
i_out <- t(idcs[i,])
out <- rbind(out,
data.frame(t = idcs[i,1],
value_real = uw$value$real[ i_out ],
value_imag = uw$value$imag[ i_out ],
dvalue_real = uw$dvalue$real[ i_out ],
dvalue_imag = uw$dvalue$imag[ i_out ])
)
}
rownames(out) <- NULL
return(invisible(out))
}
}
#' @title Print summary of data contained in `raw_cf` container
#' @param x `raw_cf` container with data and meta-data
#' @param ... ignored
#'
#' @return
#' See \link{summary.raw_cf}.
#'
#' @export
print.raw_cf <- function(x, ...) {
cf <- x
summary(cf, ...)
}
#' @title Store a 'raw_cf' correlator in an associative array together with a description
#' The object \code{cf} will be stored as an element of \code{cmap} under key \code{out_key}
#' in the member \code{obj} of \code{cmap}. The data frame passed via \code{desc} will be
#' appended as a row to \code{cmap[[out_key]]$map}. If \code{out_key} does not exist
#' as a key in \code{cmap}, a new element will be created. If it already exists,
#' \code{addStat.raw_cf} is called to add statistics to the existing \code{raw_cf}. Requires
#' the 'hash' package.
#' @return Since objects of class \code{hash} are passed and modified by reference, there
#' is no explicit return value. Instead, the passed \code{cmap} is modified.
#' @param cmap Object of class \code{hash} to act as storage for 'raw_cf' correlators.
#' @param cf Object of class \code{raw_cf} to be stored in \code{cmap}.
#' @param out_key String, key associated with \code{cf} object to be stored in \code{cmap}.
#' @param desc Single row data frame containing some descriptive parameters for \code{cf}.
store_correl <- function(cmap, cf, out_key, desc)
{
hash_avail <- requireNamespace("hash")
if( !hash_avail ){
stop("The 'hash' package is required to use this function!\n")
}
stopifnot( "hash" %in% class(cmap) )
stopifnot( "raw_cf" %in% class(cf) )
stopifnot( "raw_cf_meta" %in% class(cf) )
stopifnot( "raw_cf_data" %in% class(cf) )
if( hash::has.key(out_key, cmap) ){
cmap[[out_key]]$obj <- addStat.raw_cf(cmap[[out_key]]$obj, cf)
cmap[[out_key]]$map <- rbind(cmap[[out_key]]$map, desc)
} else {
cmap[[out_key]] <- list()
cmap[[out_key]]$obj <- cf
cmap[[out_key]]$map <- desc
}
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.