R/analysis_gradient_flow.R

In HISKP-LQCD/hadron: Analysis Framework for Monte Carlo Simulation Data in Physics

# these should also have references
ref_gf_scales <- list()
ref_gf_scales[["sqrt_t0_Eplaq"]] <- list(val=0.1416, dval=0.0008,
                                         label="t_0^{\\mathrm{plaq}}/a^2", 
                                         obs="tsqEplaq",
                                         obslabel="$\\langle t^2 E_{\\mathrm{plaq}}(t) \\rangle$",
                                         physobslabel="\\langle t_0^2 E_{\\mathrm{plaq}}(t_0) \\rangle")

ref_gf_scales[["sqrt_t0_Esym"]] <- list(val=0.1416, dval=0.0008,
                                        label="t_0^{\\mathrm{sym}}/a^2",
                                        obs="tsqEsym",
                                        obslabel="$\\langle t^2 E_{\\mathrm{sym}}(t) \\rangle$",
                                        physobslabel="\\langle t_0^2 E_{\\mathrm{sym}}(t_0) \\rangle")

ref_gf_scales[["w0_Wsym"]] <- list(val=0.1755, dval=0.0019,
                                   label="(w_0^{\\mathrm{sym}}/a)^2",
                                   obs="Wsym",
                                   obslabel="$\\langle W_{\\mathrm{sym}}(t) \\rangle$",
                                   physobslabel="\\langle W_{\\mathrm{sym}}(w_0^2) \\rangle")


#' analysis_gradient_flow
#'
#' @param path string. path to data files
#' @param outputbasename string. basename of output files
#' @param basename string. basename of input files, for example "gradflow"
#' @param read.data boolean. Indicates whether to read data fresh from
#'   data files or to use `basename.raw.gradflow.Rdata` instead
#' @param pl boolean. If set to `TRUE` plots will be generated
#' @param skip integer. number of measurements to skip
#' @param start integer. start value for time
#' @param scale numeric. scale factor for the MD time, should be set to
#'        the stridelength (in units of trajectories or configurations)
#'        which was used to produce the gradient flow files, such that
#'        the distance between measurements can be interpreted
#'        correctly and the reported autocorrelation times scaled appropriately.
#' @param plotsize numeric. Plot sidelength, this is passed to
#'        \code{tikz.init}.
#' @param dbg boolean. If set to `TRUE` debugging output will be provided.
#'
#' @description
#' function to analyse the gradient flow output files generated by
#' the tmLQCD software, see references.
#'
#' @references
#' K. Jansen and C. Urbach, Comput.Phys.Commun. 180 (2009) 2717-2738

#' @return
#' Nothing is returned.
#'
#' @export
analysis_gradient_flow <- function(path, outputbasename, basename="gradflow",
                                   read.data=TRUE, pl=FALSE, plotsize=4,
                                   skip=0, start=0,
                                   scale=1, dbg=FALSE) {

  # much like for the analysis of online measurements, we store summary information
  # in the list "gradflow_resultsum" with elements named by "outputbasename"
  # such that when the analysis is rerun, the entries are replaced
  resultsfile <- "gradflow.summary.RData"
  gradflow_resultsum <- list()
  if(file.exists(resultsfile)){
    message("Loading gradient flow results database from ", resultsfile, "\n")
    load(resultsfile)
  }

  if(read.data) {
    raw.gradflow <- readgradflow(path=path, skip=skip, basename=basename)
    save(raw.gradflow,file=sprintf("%s.raw.gradflow.Rdata",outputbasename),compress=FALSE)
  }else{
    warning(sprintf("Warning, reading data from %s.raw.gradflow.Rdata, if the number of samples changed, set read.data=TRUE to reread all output files\n",outputbasename))
    load(sprintf("%s.raw.gradflow.Rdata",outputbasename))
  }
  if(dbg==TRUE) print(raw.gradflow)

  t_vec <- unique(raw.gradflow$t)
  Ncol <- ncol(raw.gradflow)
  Nrow <- length(t_vec)
  # allocate some memory, for each observable, have space for the value, the error and the autocorrelation time
  # create a list with NULL rownams and sensible column names for this purpose
  subnames <- c("value","dvalue","ddvalue","tauint","dtauint")
  cnames <- colnames(raw.gradflow[,3:Ncol])
  outer.cnames <- t(outer(cnames,subnames,FUN=function(x,y) { paste(x,y,sep=".") }))
  grad.dimnames <- list()
  grad.dimnames[[1]] <- NULL
  grad.dimnames[[2]] <- c("t",as.vector(outer.cnames) )
  
  gradflow <- as.data.frame(matrix(data=NA,nrow=Nrow,ncol=(length(subnames)*(Ncol-2)+1),dimnames=grad.dimnames))
  for(i_row in 1:length(t_vec)){
    uwerr.gradflow <- apply(X=raw.gradflow[which(raw.gradflow$t==t_vec[i_row]),3:Ncol],MARGIN=2,FUN=uwerrprimary)
    summaryvec <- c(t_vec[i_row])
    for(i_col in 1:length(cnames)) {
      obs <- uwerr.gradflow[[cnames[i_col]]]
      summaryvec <- c(summaryvec, obs$value, obs$dvalue, obs$ddvalue, obs$tauint*scale, obs$dtauint*scale)
    }
    gradflow[i_row,] <- summaryvec
  }
  save(gradflow,file=sprintf("%s.result.gradflow.Rdata",outputbasename),compress=FALSE)
 
  gf_scales <- list()
  gf_latspacs <- list()
  gf_approx_scales <- list()
  gradflow_resultsum[[outputbasename]] <- NULL
  for( i in 1:length(ref_gf_scales) ){
    # extract colum names of gradflow relevant to the observable currently in the loop
    nms <- c("t",
             outer(ref_gf_scales[[i]]$obs, subnames, FUN=function(x,y){ paste(x,y,sep=".") })
             )

    val <- gradflow[,sprintf("%s.value", ref_gf_scales[[i]]$obs)]
    dval <- gradflow[,sprintf("%s.dvalue", ref_gf_scales[[i]]$obs)] 
    gf_scales[[i]] <- c(approx( x=val+dval,y=gradflow$t,xout=0.3)$y, 
                        approx( x=val, y=gradflow$t, xout=0.3 )$y, 
                        approx( x=val-dval, y=gradflow$t, xout=0.3)$y )


    ## determine which discrete value of t is closest to the scale in question
    for( tidx in 1:length(t_vec) ){
      if( t_vec[tidx] >= gf_scales[[i]][2] ){
        gf_approx_scales[[i]] <- gradflow[tidx,nms]
        colnames(gf_approx_scales[[i]]) <-  c("t", subnames)
        break
      }
    }

    ## summary information about the gradient flow scales:
    ## observable name, value, stat. errors in + and - directions
    ## we also add the autocorrelation time and error of the observable
    ## at the value of t closest to our scale
    gradflow_resultsum[[outputbasename]] <- 
      rbind(gradflow_resultsum[[outputbasename]],
            data.frame(obs=ref_gf_scales[[i]]$obs,
                       val=gf_scales[[i]][2],
                       pdval=gf_scales[[i]][2]-gf_scales[[i]][1],
                       mdval=gf_scales[[i]][3]-gf_scales[[i]][2],
                       tauint=gf_approx_scales[[i]]$tauint,
                       dtauint=gf_approx_scales[[i]]$dtauint
                       )
            )

    # sqrt of our bounds to compute the lattice spacing
    sqrt_gf_scale <- sqrt(gf_scales[[i]])

    gf_latspacs[[i]] <- list()
    for( k in 1:length(sqrt_gf_scale) ){
      gf_latspacs[[i]][[k]] <- compute_ratio(dividend=ref_gf_scales[[i]],
                                             divisor=list(val=sqrt_gf_scale[k],
                                                          dval=0))
    }
  }

  save(gradflow_resultsum,
       file=resultsfile)
   
  if(pl) {
    tikzfiles <- tikz.init(basename=sprintf("%s.gradflow",outputbasename),width=plotsize,height=plotsize)
    for( i in 1:length(ref_gf_scales) ){
      scale_obs <- ref_gf_scales[[i]]$obs
      scale_obslabel <- ref_gf_scales[[i]]$obslabel
      scale_label <- ref_gf_scales[[i]]$label

      val <- gradflow[,sprintf("%s.value", scale_obs)]
      dval <- gradflow[,sprintf("%s.dvalue", scale_obs)]

      gf_scale <- gf_scales[[i]]
      gf_latspac <- gf_latspacs[[i]]

      # set up plot
      # xlim is set to 1.25 times the value of t corresponding to the upper
      # limit of the reference scale in lattice units
      plot(x=gradflow$t, 
           y=val,
           type='n',
           xlim=c(0,1.25*gf_scale[1]),
           ylim=c(0.0,0.4),
           xlab="$t/a^2$",
           ylab=scale_obslabel,
           las=1)
      # draw errorband
      poly.col <- rgb(red=1.0,green=0.0,blue=0.0,alpha=0.6)
      poly.x <- c(gradflow$t,rev(gradflow$t))
      poly.y <- c(val+dval,rev(val-dval))
      polygon(x=poly.x,y=poly.y,col=poly.col)
      abline(h=0.3)
      abline(v=gf_scale)
      lines(x=gradflow$t, y=val)
      legend(x="topleft",
             legend=c(sprintf("$%s = %s$",
                              scale_label,
                              tex.catwitherror(x=gf_scale[2],
                                               # for the error, we use the average of the plus and minus errors
                                               dx=0.5*(abs(gf_scale[3]-gf_scale[2]) + abs(gf_scale[2]-gf_scale[1])),
                                               digits=2,
                                               with.dollar=FALSE, with.cdot=FALSE)
                              ),
                    sprintf("$a=%s$\\,fm", 
                            tex.catwitherror(x=gf_latspac[[2]]$val,
                                             dx=sqrt( (0.5* ( 
                                                             abs( gf_latspac[[3]]$val-gf_latspac[[2]]$val) +
                                                             abs( gf_latspac[[1]]$val-gf_latspac[[2]]$val) 
                                                             ) 
                                                       )^2 + gf_latspac[[2]]$dval^2 ),
                                             digits=2,
                                             with.dollar=FALSE, with.cdot=FALSE)
                            )
                    ),
             bty='n')
      
      ### plot MD history of the scale observable at the value of t closest to the scale
      approx_idx <- which(raw.gradflow$t==gf_approx_scales[[i]]$t)
      tseries <- data.frame(y=raw.gradflow[approx_idx,scale_obs],
                            t=start + c( skip :(skip + 
                                                length(raw.gradflow[approx_idx,scale_obs]) - 1 ) )*scale )
      plot_timeseries(dat=tseries,
                      ylab=sprintf("%s$|_{t/a^2 = %.2f}$", 
                                   scale_obslabel,
                                   gf_approx_scales[[i]]$t),
                      titletext="")

      # indicate integrated autocorrelation time in scaled units
      legend(x="topleft",
             bty='n',
             pch=NA,
             legend=sprintf("$\\tau_{\\mathrm{int}}($ %s $) = %s$ traj.",
                            scale_obslabel,
                            tex.catwitherror(x=gf_approx_scales[[i]]$tauint,
                                             dx=gf_approx_scales[[i]]$dtauint,
                                             digits=2,
                                             with.dollar=FALSE, with.cdot=FALSE)
                            )
             )
    }
   
    if( any(cnames == "Qsym") ){
      ################ TOPOLOGICAL CHARGE ####################
      # set up plot
      plot(x=gradflow$t, 
           y=gradflow$Qsym.value,
           type='n',
           ylim=range(c(gradflow$Qsym.value+gradflow$Qsym.dvalue, gradflow$Qsym.value-gradflow$Qsym.dvalue)),
           xlim=c(0.01,max(gradflow$t)),
           xlab="$t/a^2$",
           ylab="$Q(t)$",
           las=1,
           log='x')
      # draw errorband
      poly.col <- rgb(red=1.0,green=0.0,blue=0.0,alpha=0.6)
      poly.x <- c(gradflow$t,rev(gradflow$t))
      poly.y <- c(gradflow$Qsym.value+gradflow$Qsym.dvalue,rev(gradflow$Qsym.value-gradflow$Qsym.dvalue))
      polygon(x=poly.x,y=poly.y,col=poly.col)
      lines(x=gradflow$t, y=gradflow$Qsym.value)
      
      # plot MD history of Q at maximal flow time
      tmax_md_idx <- which(raw.gradflow$t==max(raw.gradflow$t))
      tmax_idx <- which(gradflow$t==max(gradflow$t))
      tseries <- data.frame(y=raw.gradflow[tmax_md_idx,"Qsym"],
                            t=start + c( skip :( skip + length(raw.gradflow[tmax_md_idx,"Qsym"]) - 1 ) )*scale)
      plot_timeseries(dat=tseries,
                      ylab=sprintf("$Q\\left( t/a^2 = %.2f \\right)$",max(raw.gradflow$t)),
                      titletext="")

      # indicate integrated autocorrelation time in scaled units
      legend(x="topleft",
             bty='n',
             pch=NA,
             legend=sprintf("$\\tau_{\\mathrm{int}}($ %s $) = %s$ traj.",
                            sprintf("$Q\\left( t/a^2 = %.2f \\right)$",max(raw.gradflow$t)),
                            tex.catwitherror(x=gradflow[tmax_idx, "Qsym.tauint"],
                                             dx=gradflow[tmax_idx, "Qsym.dtauint"],
                                             digits=1,
                                             with.dollar=FALSE, with.cdot=FALSE)
                            )
             )
    }
    
    tikz.finalize(tikzfiles)
  }
}