# R/plotutils.R In hadron: Analysis Framework for Monte Carlo Simulation Data in Physics

#### Documented in compute.plotlimsplot.averxplot.cfitplothlinewitherrorplot.massfitplot.ofitplotwitherrorpointswithslantederror

#' compute.plotlims
#'
#' @description
#' Computes limits for plots
#'
#' @param val Numeric. Value.
#' @param logscale Boolean.
#' @param cumul.dval Numeric. Cumulative error.
#' @param cumul.mdval Numeric. Cumulative error.
#'
#' @return
#' The computed plot limits are returned as a two component
#' numeric vector.
compute.plotlims <- function(val, logscale, cumul.dval, cumul.mdval){
tmp <- val - 0.1*abs(val)
tmpp <- val + 0.1*abs(val)
if(!is.null(cumul.dval)) {
## cumul.mdx is implicitly negative
tmp <- val+2*apply(X=cumul.mdval,MARGIN=1,FUN=min,na.rm=TRUE)
tmpp <- val+2*apply(X=cumul.dval,MARGIN=1,FUN=max,na.rm=TRUE)
}
if(logscale) {
tmp <- tmp[ tmp > 0 ]
tmpp <- tmpp[ tmpp > 0 ]
if( ( all(is.na(tmp)) && all(is.na(tmpp)) ) | ( length(tmp) == 0 & length(tmpp) == 0 ) ){
warning("compute.plotlims: log scale requested but there are no positive data, setting default range\n")
tmp <- 10^(-6)
tmpp <- 10^2
}
}
range(c(as.vector(tmp),as.vector(tmpp)),na.rm=TRUE)
}

## there are two possibilities for one-dimensional vectors: the vector class or the array class
is.vectorial <- function(x) {
( is.vector(x) || length(dim(x))==1 )
}

errorpos <- function(dx,errsum.method="linear") {
stop(sprintf("errorpos: called with unknown errsum.method %s",errsum.method))
}
norm <- 1

norm <- apply(X=dx,MARGIN=1,FUN=function(x){sqrt(sum(x^2))})
}

rval <- NULL
rval <- array(dim=c(nrow(dx),(ncol(dx)+2)))
} else {
rval <- array(dim=c(nrow(dx),(ncol(dx)+1)))
}

rval[,1] <- 0
rval[,2] <- dx[,1]
rval[,2] <- dx[,1]^2
}

## compute error bar deviations from the columns of dx depending on the errsum.method
## so for the "quadrature" method we will have
# norm[j] <- sqrt(sum(dx[j,]^2))
# rval[j,] == c(0, dx[j,1]^2/norm[j], (dx[j,1]^2+dx[j,2]^2)/norm[j], (dx[j,1]^2+dx[j,2]^2+dx[j,3]^2)/norm[j],..,norm[j])
if(ncol(dx)>1){
for( i in 2:ncol(dx) ){
## when dx has only one row, dx[,1:i] is returned as a vector, making apply fail....
if(nrow(dx)>1){
rval[,(i+1)] <- apply( X=dx[,1:i],MARGIN=1,
FUN=function(x){
else return( sum(x) )
} )
} else {
rval[,(i+1)] <- sum(dx[,1:i]^2)
} else {
rval[,(i+1)] <- sum(dx[,1:i])
}
}
}
}
## unlike above, even if dx has only one row, this works fine
rval[,(ncol(dx)+2)] <- apply(X=dx,MARGIN=1,FUN=function(x){ sqrt(sum(x^2)) })
}
rval <- rval/norm
rval
}

#' Plot Command For XY Plots With Error Bars
#'
#' Plot command for XY scatterplots based on plot and points which provides
#' support for multiple, non-symmetric error bars. Error bars are drawn as
#' vertical or horizontal lines originating from the point with narrow,
#' perpendicular lines at the end of the error bar (end caps). When multiple
#' errors are drawn, the width of the perpendicular line increases from the
#' innermost error bar to the outermost one. Different summation methods for
#' the individual errors are supported.
#'
#'
#' @param x vector of x coordinates
#' @param xlim limits for x-axis
#' @param y vector of y coordinates
#' @param ylim limits for y-axis
#' @param col colour of plotted data
#' @param dy one of: \itemize{ \item Vector of errors on y coordinates.
#' \item Array, matrix or data frame if multiple error bars are to be drawn,
#' such that each column refers to one error. The individual errors should be
#' provided as is, because they are summed internally to draw the final error
#' bars.  A given column can also be provided with 0 entries, in which case the
#' error bar will be drawn, but it will have zero length, such that only the
#' end caps for this error will be visible.  }
#' @param dx Same as \code{dy}, but for the x coordinate.
#' @param mdx Support for non-symmetric error bars. Same as \code{dx}, but for
#' errors in the negative x-direction. Errors should be provided as positive
#' numbers, the correct sign will be added internally. If not provided,
#' \code{dx} is used as a symmetric error.
#' @param mdy Same as \code{mdx} but for the y coordinate.
#' @param errsum.method Determines how the invidual errors should be summed for
#' display purposes. Valid argument values are: \itemize{ \item "linear"
#' \itemize{ \item Individual errors are summed linearly, such that the distance
#' from the point to the \eqn{i}'th error bar, \eqn{l_i}, is \deqn{ l_i =
#' \sum_{j=1}^i e_j } Hence, the third error bar, for example, would be located
#' at \deqn{ l_3 = e_1 + e_2 + e_3 } while the second error bar is at \deqn{
#' l_2 = e_1 + e_2 }
#'
#' }
#'
#' and error bars are drawn at the fractional position according to the
#' following formula: \deqn{ l_{max} = \sqrt{ \sum_{j=1}^{max} e_j^2 } } \deqn{
#' l_i = \sum_{j=1}^i e_j^2 / l_{max} }
#'
#' }
#'
#' \item "linear.quadrature" \itemize{ \item Errors are summed as for "linear",
#' but the total error summed in quadrature is also indicated as an end cap of
#' triple line width }
#'
#' }
#' @param rep If set to \code{TRUE}, operate like "replot" in gnuplot. Allows
#' adding points with error bars to the current plot. Switches the underlying
#' plotting routine from \code{\link[base]{plot}} to
#' @param ...  any graphic options passed over to \code{\link[base]{plot}}
#' @return a plot with error bars is drawn on the current device
#' @author Carsten Urbach, \email{urbach@@hiskp.uni-bonn.de} \cr Bartosz
#' Kostrzewa, \email{bartosz.kostrzewa@@desy.de}
#'
#' @return
#' Returns for convenience a list with elements xlim and ylim representing the
#' x- and y-limits chosen by the routine.
#'
#' @examples
#'
#' # Create some random data, set one error to zero.
#' x <- 1:50
#' y <- runif(50, 0, 1)
#' dy <- runif(50, 0.1, 0.2)
#' dy[4] <- 0
#'
#' plotwitherror(x, y, dy)
#'
#' @export plotwitherror
plotwitherror <- function(x, y, dy, ylim = NULL, dx, xlim = NULL, mdx, mdy, errsum.method="linear.quadrature", rep=FALSE, col="black", ...) {
if(!missing(mdy) && missing(dy)){
stop("plotwitherror: if 'mdy' is provided, 'dy' must be too (it can be 0)")
}
if(!missing(mdx) && missing(dx)){
stop("plotwitherror: if 'mdx' is provided, 'dx' must be too (it can be 0)")
}
if( (!missing(dx) && is.null(ncol(dx)) && length(dx)>length(x) ) || ( !missing(mdx) && is.null(ncol(mdx)) && length(mdx)>length(x) ) ||
(!missing(dy) && is.null(ncol(dy)) && length(dy)>length(y) ) || ( !missing(mdy) && is.null(ncol(dy)) && length(dy)>length(y) ) ) {
stop("plotwitherror: if any dx, mdx, dy or mdy is a simple vector, it must be of the same length as the corresponing x/y, multiple errors must ALWAYS be specified in the form of arrays/data frames/matrices")
}

## if no plotting limits are passed, we calculate them automatically
## for this to work properly, we need to know if any of the axes
## are requested in log scale
dots <- list(...)
ylog <- FALSE
xlog <- FALSE
if( 'log' %in% names(dots) ) {
logidx <- which(names(dots) == 'log')
if(any( grepl('y', dots[logidx]) ) ){ ylog <- TRUE }
if(any( grepl('x', dots[logidx]) ) ){ xlog <- TRUE }
}
my.xlim <- c()
my.ylim <- c()

## cumulative errors as computed with errsum.method
cumul.dx <- NULL
cumul.mdx <- NULL
cumul.dy <- NULL
cumul.mdy <- NULL

## compute cumulative error bar positions, first convert vectorial dx, dy into arrays
## see above for description of what errorpos does
if(!missing(dx)){
## if dx is a simple vector, convert into an appropriately sized array
if(is.vectorial(dx)){
dx <- array(data=dx,dim=c(length(dx),1))
}

cumul.dx <- errorpos(dx,errsum.method)

if(missing(mdx)){
cumul.mdx <- -errorpos(dx,errsum.method)
} else {
if(is.vectorial(mdx)){
mdx <- array(data=mdx,dim=c(length(mdx),1))
}
cumul.mdx <- -errorpos(mdx,errsum.method)
}
}

if(!missing(dy)){
if(is.vectorial(dy)){
dy <- array(data=dy,dim=c(length(dy),1))
}
cumul.dy <- errorpos(dy,errsum.method)

if(missing(mdy)){
cumul.mdy <- -errorpos(dy,errsum.method)
} else {
if(is.vectorial(mdy)){
mdy <- array(data=mdy,dim=c(length(mdy),1))
}
cumul.mdy <- -errorpos(mdy,errsum.method)
}
}

if( is.null(xlim) ){
my.xlim <- compute.plotlims(val=x, logscale=xlog, cumul.dval=cumul.dx, cumul.mdval=cumul.mdx)
} else {
my.xlim <- xlim
}

if( is.null(ylim) ){
my.ylim <- compute.plotlims(val=y, logscale=ylog, cumul.dval=cumul.dy, cumul.mdval=cumul.mdy)
} else {
my.ylim <- ylim
}

if(rep) {
points(x, y, col=col, ...)
}
else if(!is.null(my.xlim) && !is.null(my.ylim)) {
plot(x, y, ylim=my.ylim, xlim=my.xlim, col=col, ...)
}
else if(is.null(my.xlim) && !is.null(my.ylim)) {
plot(x, y, ylim=my.ylim, col=col, ...)
}
else if(!is.null(my.xlim) && is.null(my.ylim)) {
plot(x, y, xlim=my.xlim, col=col, ...)
}
else {
plot(x, y, col=col, ...)
}

## We will need some length in inches, so let's define the proportionality factor as
## size of the output in inches divided by the difference of boundary coordinates.
x.to.inches <- par("pin")[1]/diff(par("usr")[1:2])
y.to.inches <- par("pin")[2]/diff(par("usr")[3:4])

if(!is.null(cumul.dy)) {
for(cumul.err in list(cumul.dy,cumul.mdy)){
rng <- 2:ncol(cumul.err)
if(ncol(cumul.err)>2 && errsum.method=="linear.quadrature") rng <- 2:(ncol(cumul.err)-1)
## this loop is necessary because the "length" parameter of "arrows" is not vectorial...
## so to accommodate the generalisations below, we need to draw the error for each point
## individually, it doesn't make it much slower
## the length of the arrowhead lines will depend on the "level" (the more errors, the longer the arrowhead lines)
arwhd.len <- 0.02
for(level in rng){
start <- y+cumul.err[,(level-1)]
end <- y+cumul.err[,level]
proper.val <- is.finite(start) & is.finite(end)
start <- start[proper.val]
end <- end[proper.val]

if(par("ylog")){
## logarithmic scaling means distances are multiplicative (the inner abs() is just supposed to prevent irrelevant warnings)
dy.in.inches <- ifelse(start > 0 & end > 0, abs(log10(abs(end/start)))*y.to.inches, 0)
## An arrow can only be plotted if it is longer than 1/1000 inches.
## We want to plot an errorbar of minimum size anyway, so as to show the point is plotted with an error.
start <- ifelse(dy.in.inches <= 1/1000, y[proper.val] / 10^(1.01/2000/y.to.inches), start)
end <- ifelse(dy.in.inches <= 1/1000, y[proper.val] * 10^(1.01/2000/y.to.inches), end)
}else{
dy.in.inches <- abs(end-start)*y.to.inches
start <- ifelse(dy.in.inches <= 1/1000, y[proper.val] - 1.01/2000/y.to.inches, start)
end <- ifelse(dy.in.inches <= 1/1000, y[proper.val] + 1.01/2000/y.to.inches, end)
}

if(length(col) == 1){
rw <- which(dy.in.inches > 1/1000)
arrows(x[proper.val][rw], start[rw], x[proper.val][rw], end[rw], length=arwhd.len, angle=90, code=2, col=col)
rw <- which(dy.in.inches > 0)
arrows(x[proper.val][rw], start[rw], x[proper.val][rw], end[rw], length=arwhd.len, angle=90, code=3, col=col)
}else{
## this loop is necessary because the "col" parameter of "arrows" is not vectorial...
## it does make it much slower
for(rw in seq(proper.val)){
if(dy.in.inches[rw] > 1/1000)
arrows(x[proper.val][rw], start[rw], x[proper.val][rw], end[rw], length=arwhd.len, angle=90, code=2, col=col[proper.val][rw])
else if(dy.in.inches[rw] > 0)
arrows(x[proper.val][rw], start[rw], x[proper.val][rw], end[rw], length=arwhd.len, angle=90, code=3, col=col[proper.val][rw])
}
}

arwhd.len <- arwhd.len + 0.01
}
## for the linear.quadrature method, show the total error as a line of triple thickness
## without drawing any "arrowstems"
## to be consistent, drawX/Ybars uses inches just like arrows
arwhd.len <- arwhd.len + 0.02
if(length(col) == 1)
drawYbars(x=x,y=y,dy=cumul.err[,ncol(cumul.err)],length=arwhd.len,lwd=3,col=col)
else{
for(rw in seq(y))
drawYbars(x=x[rw],y=y[rw],dy=cumul.err[rw,ncol(cumul.err)],length=arwhd.len,lwd=3,col=col[rw])
}
}
}
}
if(!is.null(cumul.dx)) {
for(cumul.err in list(cumul.dx,cumul.mdx)){
rng <- 2:ncol(cumul.err)
if(ncol(cumul.err)>2 && errsum.method=="linear.quadrature") rng <- 2:(ncol(cumul.err)-1)
arwhd.len <- 0.02
for(level in rng){
start <- x+cumul.err[,(level-1)]
end <- x+cumul.err[,level]
proper.val <- is.finite(start) & is.finite(end)
start <- start[proper.val]
end <- end[proper.val]

if(par("xlog")){
dx.in.inches <- ifelse(start > 0 & end > 0, abs(log10(abs(end/start)))*x.to.inches, 0)
start <- ifelse(dx.in.inches <= 1/1000, x[proper.val] / 10^(1.01/2000/x.to.inches), start)
end <- ifelse(dx.in.inches <= 1/1000, x[proper.val] * 10^(1.01/2000/x.to.inches), end)
}else{
dx.in.inches <- abs(end-start)*x.to.inches
start <- ifelse(dx.in.inches <= 1/1000, x[proper.val] - 1.01/2000/x.to.inches, start)
end <- ifelse(dx.in.inches <= 1/1000, x[proper.val] + 1.01/2000/x.to.inches, end)
}

if(length(col) == 1){
rw <- which(dx.in.inches > 1/1000)
arrows(start[rw], y[proper.val][rw], end[rw], y[proper.val][rw], length=arwhd.len, angle=90, code=2, col=col)
rw <- which(dx.in.inches > 0)
arrows(start[rw], y[proper.val][rw], end[rw], y[proper.val][rw], length=arwhd.len, angle=90, code=3, col=col)
}else{
for(rw in seq(proper.val)){
if(dx.in.inches[rw] > 1/1000)
arrows(start[rw], y[proper.val][rw], end[rw], y[proper.val][rw], length=arwhd.len, angle=90, code=2, col=col[proper.val][rw])
else if(dx.in.inches[rw] > 0)
arrows(start[rw], y[proper.val][rw], end[rw], y[proper.val][rw], length=arwhd.len, angle=90, code=3, col=col[proper.val][rw])
}
}

arwhd.len <- arwhd.len + 0.01
}
## to be consistent, drawX/Ybars uses inches just like arrows
arwhd.len <- arwhd.len + 0.02
if(length(col) == 1)
drawXbars(x=x,y=y,dx=cumul.err[,ncol(cumul.err)],length=arwhd.len,lwd=3,col=col)
else{
for(rw in seq(x))
drawXbars(x=x[rw],y=y[rw],dx=cumul.err[rw,ncol(cumul.err)],length=arwhd.len,lwd=3,col=col[rw])
}
}
}
}

return(invisible(list(xlim=my.xlim, ylim=my.ylim)))
}

#' plothlinewitherror
#'
#' @description
#' plot a horizontal line with error band
#'
#' @param m Numeric. Mean value of the line to plot.
#' @param dp Numeric. Error up.
#' @param dm Numeric. Error down.
#' @param col String. Colour.
#' @param x0 Numeric. Left value of the range of the horizontal line.
#' @param x1 Numeric. Right value of the range of the horizontal line.
#'
#' @return
#' No return value, only graphics is generated.
#'
#' @export
plothlinewitherror <- function(m, dp, dm, col=c("red"), x0, x1) {
if(missing(dm)) {
dm <- dp
}
arrows(x0=x0, y0=m, x1=x1, y1=m, col=col, length=0)
arrows(x0=x0, y0=m+dp, x1=x1, y1=m+dp, col=col, length=0, lwd=c(1))
arrows(x0=x0, y0=m-dm, x1=x1, y1=m-dm, col=col, length=0, lwd=c(1))
}

#' plot.massfit
#'
#' @description
#' Generic function to plot an object of type massfit
#'
#' @param x Object of type massfit
#' @param ... Generic graphical parameter to be passed on to \link{plotwitherror}
#' @param xlab String. Label for x-axis
#' @param ylab String. Lable for y-axis
#'
#' @return
#'
#' @export
plot.massfit <- function(x, ..., xlab = "t", ylab = "m") {
plotwitherror(x$t, x$mass, x$dmass, xlab=xlab, ylab=ylab, ...) } #' plot.c1fit #' #' @description #' Generic function to plot an object of type c1fit #' #' @param x Object of type c1fit #' @param ... Generic graphical parameter to be passed on to \link{plotwitherror} #' #' @return #' No return value, only plots are generated. #' #' @export plot.cfit <- function(x, ...) { fit <- x fit.mass <- abs(fit$fitresult$par[fit$matrix.size+1])
if(!is.null(fit$effmass$mll)) {
new_window_if_appropriate()
plot.effmass(m=fit.mass,
ll=data.frame(t=fit$effmass$t, mass=fit$effmass$mll, dmass=fit$effmass$dmll),
lf=data.frame(t=fit$effmass$t, mass=fit$effmass$mlf, dmass=fit$effmass$dmlf),
ff=data.frame(t=fit$effmass$t, mass=fit$effmass$mff, dmass=fit$effmass$dmff))
}
if(!is.null(fit$effmass$m)) {
new_window_if_appropriate()
plot.effmass(m=fit.mass, ll=data.frame(t=fit$effmass$t, mass=fit$effmass$m, dmass=fit$effmass$dm))
}
if(!is.null(fit$uwerrresultmpcac)) { new_window_if_appropriate() plot(fit$uwerrresultmpcac, main=expression(m[PCAC]))
}
if(!is.null(fit$uwerrresultmps)) { new_window_if_appropriate() plot(fit$uwerrresultmps, main=expression(m[PS]))
}
if(!is.null(fit$uwerrresultfps)) { new_window_if_appropriate() plot(fit$uwerrresultfps, main=expression(f[PS]))
}
if(!is.null(fit$uwerrresultmv)) { new_window_if_appropriate() plot(fit$uwerrresultmv, main=expression(m[V]))
}
if(!is.null(fit$boot)) { new_window_if_appropriate() plot(fit$boot, main="Bootstrap analysis for mps")
}
if(!is.null(fit$tsboot)) { new_window_if_appropriate() plot(fit$tsboot, main="TS Boostrap analysis for mps")
}
if(!is.null(fit$mv.boot)) { new_window_if_appropriate() plot(fit$mv.boot, main="Bootstrap analysis for mv")
}
if(!is.null(fit$mv.tsboot)) { new_window_if_appropriate() plot(fit$mv.tsboot, main="TS Bootstrap analysis for mv")
}
if(!is.null(fit$fitdata)) { new_window_if_appropriate() plot(fit$fitdata$Chi, main="Chi per data point", xlab="data point", ylab=expression(chi), type="h") } if(!is.null(fit$dpaopp)) {
new_window_if_appropriate()
plotwitherror(fit$dpaopp$t, fit$dpaopp$mass, fit$dpaopp$dmass,
main=expression(m[PCAC]), xlab="t", ylab=expression(m[PCAC]))
abline(h=fit$fitresult$par[3]*fit$fitresult$par[2]/fit$fitresult$par[1]/2., col="red")
}
if(!is.null(fit$MChist.dpaopp)) { new_window_if_appropriate() plot(seq(1, length(fit$MChist.dpaopp))+fit$skip, fit$MChist.dpaopp, type="l",
main=expression(m[PCAC]), xlab=expression(t[HMC]), ylab=expression(m[PCAC]))
abline(h=fit$fitresult$par[3]*fit$fitresult$par[2]/fit$fitresult$par[1]/2., col="red")
}
}

#' plot.ofit
#'
#' @description
#' Generic function to plot an object of type ofit
#'
#' @param x Object of type ofit
#' @param ... Generic graphical parameter to be passed on to \link{plotwitherror}
#'
#' @return
#'
#' @export
plot.ofit <- function(x, ...) {
plot.cfit(x)
}

#' plot.effmass
#'
#' @param x Object of class effmass
#' @param ... Graphical parameters to be passed on.
#' @param ll local-local effective mass object
#' @param lf local-fuzzed effective mass object
#' @param ff fuzzed-fuzzed effective mass object
#'
#' @return
#' No value returned, only plots are generated.
#'
#' @export
plot.effmass <- function (x, ..., ll, lf, ff) {
m <- x

new_window_if_appropriate()
if(!missing(ff)) {
plot.massfit(ff, ylab=expression(m[eff]), xlab="t", ...)
points((ll$t-0.2), ll$mass, pch=1, col="blue")
arrows((ll$t-0.2), ll$mass-ll$dmass, (ll$t-0.2), ll$mass+ll$dmass, length=0.01,angle=90,code=3)
points((lf$t+0.2), lf$mass, pch=2, col="red")
arrows((lf$t+0.2), lf$mass-lf$dmass, (lf$t+0.2), lf$mass+lf$dmass, length=0.01,angle=90,code=3)
lines(ll$t, rep(m, times=length(ll$t)))
}
else if(!missing(lf)) {
plot.massfit(ll, ylab=expression(m[eff]), xlab="t", ...)
points((lf$t-0.2), lf$mass, pch=1, col="blue")
arrows((lf$t-0.2), lf$mass-lf$dmass, (lf$t-0.2), lf$mass+lf$dmass, length=0.01,angle=90,code=3)
lines(ll$t, rep(m, times=length(ll$t)))
}
else {
plot.massfit(ll, ylab=expression(m[eff]), xlab="t", ...)
lines(ll$t, rep(m, times=length(ll$t)))
}
}

#' Plots averx data
#'
#' @param x averx object
#' @param ... ignored
#'
#' @return
#' Returns the plotted data in from of a \link{data.frame} with named
#' columns t (the time index), averx the values of average x and
#' daverx the statistical error estimate.
#'
#' @export
plot.averx <- function(x, ...) {
averx <- x
Thalfp1 <- averx$Cf2pt$Time/2+1
#plot(averx$effmass, ylim=c(averx$effmass$opt.res$par[1]/2, 3/2*averx$effmass$opt.res$par[1]), main=c("Pion Effectivemass"), xlab=c("t/a"), ylab=c("a Meff")) new_window_if_appropriate() plot(averx$Cf3pt, xlab=c("t/a"), ylab=c("C3pt"), ylim=c(0, 2*averx$plateau), main=c("3pt ov 2pt")) plothlinewitherror(m=averx$plateau, dp=sd(averx$plateau.tsboot[,1]), dm=sd(averx$plateau.tsboot[,1]),
x0=averx$t1, x1=averx$t2)
new_window_if_appropriate()
plotwitherror(c(0:(averx$Cf2pt$Time/2)),
averx$Cf3pt$cf0/averx$matrixfit$opt.res$par[1]/averx$Cf2pt$cf0[Thalfp1], apply(averx$Cf3pt$cf.tsboot$t/averx$matrixfit$opt.tsboot[1,]/averx$Cf2pt$cf.tsboot$t[,Thalfp1], 2, sd), ylim=c(0,0.6), xlab=c("t/a"), ylab=c("C3pt"), main=c("<x>") ) new_window_if_appropriate() qqplot(averx$plateau.tsboot[,2], rchisq(n=averx$boot.R, df=averx$dof), main=paste("qqplot chisq"))
return(invisible(data.frame(t=c(0:(averx$Cf2pt$Time/2)),
averx=averx$Cf3pt$cf0/averx$matrixfit$opt.res$par[1]/averx$Cf2pt$cf0[Thalfp1], daverx=apply(averx$Cf3pt$cf.tsboot$t/averx$matrixfit$opt.tsboot[1,]/averx$Cf2pt$cf.tsboot$t[,Thalfp1], 2, sd) ) ) ) } #' plot.pionff #' #' @description #' Generic function to plot an object of type pionff #' #' @param x Object of type pionff #' @param ... Generic graphical parameter to be passed on to \link{plotwitherror} #' #' @return #' No return value, only plots are generated. #' #' @export plot.pionff <- function (x, ...) { ff <- x Time <- ff$Cf2ptp0$Time Thalfp1 <- Time/2+1 plot(mul.cf(ff$Cf3ptp0, 1./ff$Cf2ptp0$cf0[Thalfp1]), main=c("1./Z_V"), xlab=c("t/a"), ylab=c("1/Z_V"))
plothlinewitherror(m=1./ff$Cf2ptp0$cf0[Thalfp1]*ff$plateaufitZV$plateau, dp=sd(ff$plateaufitZV$plateau.tsboot[,1]/ff$Cf2ptp0$cf.tsboot$t[,Thalfp1]), dm=sd(ff$plateaufitZV$plateau.tsboot[,1]/ff$Cf2ptp0$cf.tsboot$t[,Thalfp1]),
x0=ff$plateaufitZV$t1, x1=ff$plateaufitZV$t2)

new_window_if_appropriate()
plot(mul.cf(ff$Cf3ptratio, ff$Cf2ptratio$cf0[Thalfp1]), ylim=c(0,1.3), main=c("F(q^2)"), xlab=c("t/a"), ylab=c("F(q^2)")) plothlinewitherror(m=ff$plateaufitFF$plateau*ff$Cf2ptratio$cf0[Thalfp1], dp=sd(ff$plateaufitFF$plateau.tsboot[,1]*ff$Cf2ptratio$cf.tsboot$t[,Thalfp1]), dm=sd(ff$plateaufitFF$plateau.tsboot[,1]*ff$Cf2ptratio$cf.tsboot$t[,Thalfp1]), x0=ff$plateaufitFF$t1, x1=ff$plateaufitFF$t2) } #' Plot Command For Class Ouputdata #' #' Generic plot routine for class \dQuote{ouputdata}. Currently it plots the #' plaquette history and the history of \eqn{\Delta H}{Delta H} #' #' #' @param x object of class \dQuote{outputdata} obtained from a read with #' \code{readoutputdata} #' @param skip number of trajectories to be skipped in analysis for plaquette #' and \eqn{\exp(-\Delta H)}{exp(-Delta H)}. #' @param ... additional arguments passed to the generic plot function. #' @return list containing the \dQuote{data}, an object of class \dQuote{uwerr} #' called \dQuote{plaq.res} containing the statisical analysis for the #' plaquette and a second object of type \dQuote{uwerr} called \dQuote{dH.res} #' with the statisical analysis for \eqn{\exp(-\Delta }{exp(-Delta H)}\eqn{ #' H)}{exp(-Delta H)}. #' @author Carsten Urbach, \email{curbach@@gmx.de} #' @seealso \code{\link{readoutputdata}}, \code{\link{uwerr}} #' @keywords methods hplot #' @return #' The plotted data is return in form of a \link{list} with named elements data #' containing the input data, plaq.res an object returned by \link{uwerrprimary} #' for the plaquette data dn dH.res an object returned by \link{uwerrprimary} #' for \eqn{\Delta H}{Delta H}. #' #' @export #' @examples #' #' plaq <- readoutputdata(paste0(system.file(package="hadron"), "/extdata/output.data")) #' plaq.plot <- plot(plaq, skip=100) #' summary(plaq.plot$plaq.res)
#'
plot.outputdata <- function (x, skip = 0, ...) {
data <- x
plaq.res <- uwerrprimary( data$V2[skip:length(data$V2)])
dH.res <- uwerrprimary( exp(-data$V3[skip:length(data$V3)]))
plot(data$V1, data$V2, type="l",
main="Plaquette", xlab=expression(t[HMC]), ylab="P", col="red", ...)
abline(h=plaq.res$value, col="blue") abline(h=plaq.res$value+plaq.res$dvalue) abline(h=plaq.res$value-plaq.res$dvalue) new_window_if_appropriate() plot(data$V1, data\$V3, type="l",
main=expression(paste(Delta, "H")), xlab=expression(t[HMC]), ylab=expression(paste(Delta, "H")), ...)
return(invisible(list(data=data, plaq.res=plaq.res, dH.res = dH.res)))
}

## draw Y (X) error bars at coordinates y+dy (x+dx) (where dy (dx) can be negative)
## like for arrows, the bar width is specified in inches
## and additional parameters (like lwd and col) can be passed
## to segments
drawYbars <- function(x,y,dy,length=0.01,...) {
x.inch <- grconvertX(x=x,from="user",to="inches")
xp.inch <- x.inch+length
xm.inch <- x.inch-length
xp <- grconvertX(x=xp.inch,from="inches",to="user")
xm <- grconvertX(x=xm.inch,from="inches",to="user")
segments(xp, y+dy,
xm, y+dy, ...)
}

drawXbars <- function(x,y,dx,length=0.01,...) {
y.inch <- grconvertY(y=y,from="user",to="inches")
yp.inch <- y.inch+length
ym.inch <- y.inch-length
yp <- grconvertY(y=yp.inch,from="inches",to="user")
ym <- grconvertY(y=ym.inch,from="inches",to="user")
segments(x+dx, ym,
x+dx, yp, ...)
}

new_window_if_appropriate <- function () {
is_X11 <- grepl(pattern = "X11", x = names(dev.cur()), ignore.case = TRUE)
is_null <- grepl(pattern = "null", x = names(dev.cur()), ignore.case = TRUE)

if (interactive() && (is_X11 || is_null)) {
message('Opening a new X11 window.\n')
dev.new()
}
}

#' pointswithslantederror
#'
#' plots data points with slanted error bars
#'
#' @description
#' This function plots points with x- and y-errors visualised as a
#' slanted errorbar. The length of the error bar represents x- and y-errors
#' added in quadrature. The slope of the error bar is positive of negative
#' depending on whether the correlation betwenn x and y is positive or
#' negative, respectively.
#'
#' @param x numeric vector. x-values
#' @param y numeric vector. y-values
#' @param dx numeric vector. x-standard errors
#' @param dy numeric vector. y-standard errors
#' @param cor numeric vector. Correlation coefficients between x- and y-
#' errors.
#' @param col the color of the points
#' @param bcol the color of the slanted error bars
#' @param ... further graphical parameters to be passed on to points
#'
#' @examples
#' x <- c(1:5)
#' y <- x^2
#' dx <- c(0.1, 0.2, 0.2, 0.1, 0.05)
#' dy <- c(0.05, 0.2, 0.1, 0.2, 0.1)
#' cor <- c(1, -1, -1, 1, 1)
#' plot(NA, xlim=range(x), ylim=range(y), xlab="y", ylab="y")
#' pointswithslantederror(x=x, y=y, dx=dx, dy=dy, cor=cor)
#'
#' @export
pointswithslantederror <- function(x, y, dx, dy, cor, col="black", bcol="black", ...) {
points(x=x, y=y, col=col, ...)
l <- sqrt(dx^2+dy^2)
fac <- rep(1, times=length(cor))
fac[which(cor < 0)] <- -1
arrows(code=0,
x0=x-fac*dx, y0=y-dy,
x1=x+fac*dx, y1=y+dy,
col=bcol)
}