Origin/reweighting

DESCRIPTION

@@ -15,6 +15,7 @@ Maintainer: Carsten Urbach <[email protected]>
 SystemRequirements: Gnu Scientific Library version >= 1.8
 Description: Toolkit to extract hadronic quantities from Lattice QCD simulations. It contains functionality for IO, plotting, bootstrap and jackknife resampling, fitting, GEVP solving, error and autocorrelation estimation as well as other areas.
 Imports:
+    dplyr,
     abind,
     boot,
     R6,

NAMESPACE

@@ -1,5 +1,7 @@
 exportPattern("^cf.*")
 exportPattern("*.cf")
+exportPattern("^rw.*")
+exportPattern("*.rw")
 exportPattern("^plot.*")
 exportPattern("^analysis_*")
 exportPattern("^print.*")

R/cf.R

@@ -47,7 +47,8 @@ cf_meta <- function (.cf = cf(), nrObs = 1, Time = NA, nrStypes = 1, symmetrised
   return (.cf)
 }
 
-#' Bootstrapped CF mixin constructor
+
+#' Bootstrapped CF mixin constructor 
 #'
 #' @param .cf `cf` object to extend.
 #' @param boot.R Integer, number of bootstrap samples used.
@@ -108,6 +109,7 @@ cf_boot <- function (.cf = cf(), boot.R, boot.l, seed, sim, cf.tsboot, icf.tsboo
   return (.cf)
 }
 
+
 #' Estimates error from jackknife samples
 #'
 #' Currently this uses the mean over the jackknife samples in order to compute
@@ -374,6 +376,78 @@ gen.block.array <- function(n, R, l, endcorr=TRUE) {
   return(list(starts = st, lengths = lens))
 }
 
+#' Computes the samples for reweighted correlation function
+#'
+#' @param cf `cf` object.
+#' @param rw `rw` object.
+#' @param boot.R Integer
+#' @param boot.l Integer
+#' @param seed Integer
+#' @param sim string
+#' @param endcorr boolean
+#' @export
+bootstrap_rw.cf <- function(cf, rw, boot.R=400, boot.l=2, seed=1234, sim="geom", endcorr=TRUE) {
+  stopifnot(inherits(cf, 'cf_orig'))
+  stopifnot(inherits(rw, 'rw_orig'))
+  stopifnot(inherits(rw, 'rw_meta'))
+  stopifnot(inherits(cf, 'cf_indexed'))
+
+
+  ##We should also check that the cf object and the rw object contains the same gauge configurations
+
+  stopifnot(rw$conf.index == cf$conf.index)
+
+  stopifnot( nrow(cf$cf) == length(rw$conf.index) )
+
+  boot.l <- ceiling(boot.l)
+  boot.R <- floor(boot.R)
+
+  stopifnot(boot.l >= 1)
+  stopifnot(boot.l <= nrow(cf$cf))
+  stopifnot(boot.R >= 1)
+
+  ##Construct correlation function for the reweighting samples
+  rw_cf <- cf
+  rw_cf$cf <- replicate(ncol(cf$cf), rw$rw)
+
+  ## we set the seed for reproducability and correlation
+  old_seed <- swap_seed(seed)
+
+  ## now we bootstrap the correlators*reweighting factor
+  rwcf.tsboot <- boot::tsboot(cf$cf*rw_cf$cf, statistic = function(x){ return(apply(x, MARGIN=2L, FUN=mean))},
+                            R = boot.R, l=boot.l, sim=sim, endcorr=endcorr)
+
+
+  restore_seed(old_seed)
+
+  ## we set the seed for reproducability and correlation
+  old_seed <- swap_seed(seed)
+
+  ## now we bootstrap the reweighting factor
+  rw.tsboot <- boot::tsboot(rw_cf$cf, statistic = function(x){ return(apply(x, MARGIN=2L, FUN=mean))},
+                            R = boot.R, l=boot.l, sim=sim, endcorr=endcorr)
+
+
+  rwcf.tsboot$t0<- rwcf.tsboot$t0/rw.tsboot$t0
+  rwcf.tsboot$t <- rwcf.tsboot$t/rw.tsboot$t
+
+  cf <- cf_boot(cf,
+                boot.R = boot.R,
+                boot.l = boot.l,
+                seed = seed,
+                sim = sim,
+                cf.tsboot = rwcf.tsboot)
+
+  class(cf) <- append(class(cf), 'cfrw_boot')
+
+  class(cf) <- setdiff(class(cf), 'cf_orig')
+
+
+  restore_seed(old_seed)
+
+  return(invisible(cf))
+}
+
 bootstrap.cf <- function(cf, boot.R=400, boot.l=2, seed=1234, sim="geom", endcorr=TRUE) {
   stopifnot(inherits(cf, 'cf_orig'))
 
@@ -412,6 +486,7 @@ bootstrap.cf <- function(cf, boot.R=400, boot.l=2, seed=1234, sim="geom", endcor
   return(invisible(cf))
 }
 
+
 jackknife.cf <- function(cf, boot.l = 1) {
   stopifnot(inherits(cf, 'cf_orig'))
 
@@ -469,8 +544,83 @@ jackknife.cf <- function(cf, boot.l = 1) {
                 resampling_method = 'jackknife')
 
   return (invisible(cf))
+
 }
 
+#' Computes the jackknife samples for reweighted correlation function
+#'
+#' @param cf `cf` object.
+#' @param rw `rw` object.
+#' @param boot.l Integer
+#' @export
+jackknife_rw.cf <- function(cf, rw, boot.l = 1) {
+  stopifnot(inherits(cf, 'cf_orig'))
+  stopifnot(inherits(rw, 'rw_orig'))
+  stopifnot(inherits(rw, 'rw_meta'))
+  stopifnot(inherits(cf, 'cf_indexed'))
+
+  stopifnot(rw1$conf.index == rw2$conf.index)
+
+  ##We should also check that the cf object and the rw object contains the same gauge configurations
+
+  stopifnot( nrow(cf$cf) == length(rw$conf.index) )
+
+
+  stopifnot(boot.l >= 1)
+  boot.l <- ceiling(boot.l)
+
+  ##Construct correlation function for the reweighting samples
+  rw_cf <- cf
+  rw_cf$cf <- replicate(ncol(cf$cf), rw$rw)
+
+
+  ## blocking with fixed block length, but overlapping blocks
+  ## number of observations
+  n <- nrow(cf$cf)
+  ## number of overlapping blocks
+  N <- n-boot.l+1
+
+
+  numerator <- apply(cf$cf*rw_cf$cf, 2, mean)
+  denominator <- apply(rw_cf$cf, 2, mean)
+  t0 <- numerator/denominator
+
+  t <- array(NA, dim = c(N, ncol(cf$cf)))
+  for (i in 1:N) {
+    ## The measurements that we are going to leave out.
+    ii <- c(i:(i+boot.l-1))
+    ## jackknife replications of the mean
+    t[i, ] <- 
+
+    numerator <- apply(cf$cf[-ii, ]* rw_cf$cf[ ii, ], 2L, mean)
+    denominator <- apply( rw_cf$cf[ ii, ] , 2L, mean )
+    t[i, ] < numerator/denominator 
+  }
+
+
+  cf <- invalidate.samples.cf(cf)
+
+  cf.tsboot <- list(t = t,
+                    t0 = t0,
+                    R = N,
+                    l = boot.l)
+
+
+  cf <- cf_boot(cf,
+                boot.R = cf.tsboot$R,
+                boot.l = cf.tsboot$l,
+                seed = 0,
+                sim = 'geom',
+                cf.tsboot = cf.tsboot,
+                resampling_method = 'jackknife')
+
+  class(cf) <- append(class(cf), 'cfrw_boot')
+
+  class(cf) <- setdiff(class(cf), 'cf_orig')
+
+  return (invisible(cf))
+}
+# Gamma method analysis on all time-slices in a 'cf' object
 #' uwerr.cf
 #' @description
 #' Gamma method analysis on all time-slices in a 'cf' object
@@ -513,10 +663,31 @@ addConfIndex2cf <- function(cf, conf.index) {
   if(is.null(cf$conf.index)) {
     cf$conf.index <- conf.index
   }
+  class(cf) <- append(class(cf), 'cf_indexed')
   return(cf)
 }
 
-addStat.cf <- function(cf1, cf2) {
+#' Combine correlation function from different replicas
+#' 
+#' @param cf1 `cf` object: correlation function for replicum A
+#' @param cf2 `cf` object: correlation function for replicum B
+#' @param reverse1 `boolean` After the bifurcation point one of
+#'                           the replicas (chain of correlation 
+#'                           functions in simulation time) has  
+#'                           to be reversed.
+#' @param reverse2 `boolean`
+#'
+#' @examples
+#' Suppose we have correlation functions in replicum A from 0 to 500
+#' in steps of 4 and in replicum B from 4 to 500 in steps of 4.
+#' To combined the two replicas we have to use
+#'
+#' addstat.cf(cf_replicumB, cf_replicumA, TRUE, FALSE)
+#' which means
+#' combined=(cf500 from B, cf496 from B,...,cf004 from B, cf000 from A, ..
+#' cf500 from A) 
+#' @export
+addStat.cf <- function(cf1, cf2,reverse1=FALSE, reverse2=FALSE) {
   stopifnot(inherits(cf1, 'cf'))
   stopifnot(inherits(cf2, 'cf'))
 
@@ -530,15 +701,47 @@ addStat.cf <- function(cf1, cf2) {
   stopifnot(inherits(cf1, 'cf_meta'))
   stopifnot(inherits(cf2, 'cf_meta'))
 
+  ##Either both should have an index or none of them
+  stopifnot(inherits(cf1, 'cf_indexed') == inherits(cf1, 'cf_indexed') )
+
   stopifnot(cf1$Time == cf2$Time)
   stopifnot(dim(cf1$cf)[2] == dim(cf2$cf)[2])
   stopifnot(cf1$nrObs == cf2$nrObs )
   stopifnot(cf1$nrStypes == cf2$nrStypes)
 
   cf <- cf1
 
-  cf$cf <- rbind(cf1$cf, cf2$cf)
-  cf$icf <- rbind(cf1$icf, cf2$icf)
+  cf1_temp<- cf1$cf
+  icf1_temp <- cf1$icf
+  if (reverse1 == TRUE){
+    apply(cf1_temp,2,rev)
+    if ( has_icf(cf1)){
+      apply(icf1_temp,2,rev)
+    }
+  }
+  cf2_temp <- cf2$cf
+  icf2_temp <- cf2$icf
+  if (reverse2 == TRUE){
+    apply(cf2_temp,2,rev)
+    if ( has_icf(cf2)){
+      apply(icf2_temp,2,rev)
+    }
+  }
+  if (inherits(cf1, 'cf_indexed')){
+    conflist_temp1 <- cf1$conf.index 
+    if (reverse1 == TRUE){
+      conflist_temp1 <- rev(conflist_temp1) 
+    }
+    conflist_temp2 <- cf2$conf.index
+    if (reverse2 == TRUE){
+      conflist_temp2 <- rev(conflist_temp2)
+    }
+    cf$conf.index <- c(conflist_temp1,conflist_temp2)
+  }
+
+
+  cf$cf <- rbind(cf1_temp, cf2_temp)
+  cf$icf <- rbind(icf1_temp, icf2_temp)
 
   cf <- invalidate.samples.cf(cf)
 

R/readutils.R

@@ -347,7 +347,50 @@ readtextcf <- function(file, T=48, sym=TRUE, path="", skip=1, check.t=0, ind.vec
 
   return (invisible(ret))
 }
+#' @title reading reweighting factors for a list of gauge configuration 
+#'        and random samples from ASCII files
+#' @param file_names_to_read list of filenames for the reweighting factors
+#' @param gauge_conf_list <- a list of integers with the indices of the gauge configs
+#' @param nsamples number of stochastic samples used for computing the reweighting factors
+read.rw <- function( file_names_to_read, gauge_conf_list, nsamples, monomial_id ) 
+{
+  stopifnot(length(gauge_conf_list)==length(file_names_to_read)) 
+  ret <- rw_meta(conf.index=gauge_conf_list)
+  tmp <- readcmidatafiles(files=file_names_to_read,skip=0,verbose=TRUE,colClasses=c("integer","integer","numeric","numeric","numeric","numeric","numeric"))
+  names(tmp)[1] <- "monomialid"
+  names(tmp)[2] <- "stochastic_index"
+  names(tmp)[3] <- "kappa_target"
+  names(tmp)[4] <- "kappa_original"
+  names(tmp)[5] <- "light_quark_mass_target"
+  names(tmp)[6] <- "light_quark_mass_original"
+  names(tmp)[7] <- "reweightingfactor"
+
+# Select the reweighting factor for a particular monomial
+
+  dplyr_avail <- requireNamespace("dplyr")
+  stopifnot(dplyr_avail)
+
+  tmp <- dplyr::filter(tmp,monomialid==monomial_id)
+
+# Number of reweighted determinants for each gauge configuration
+
+  n_rew_factors <- length(tmp$reweightingfactor)/(nsamples*length(gauge_conf_list))
+  stopifnot(n_rew_factors == 1)
+
 
+# Exponentianing and Averaging over the stochastic samples
+
+  tmp2 <- matrix(tmp$reweightingfactor,nrow=nsamples,ncol=length(gauge_conf_list)*n_rew_factors)
+  tmp3 <- apply(exp(-tmp2),2,mean)
+
+# Normalize the largest reweighting factor to be one and storing this factor
+# this is neccessary due to the large value of the reweighting factor 
+# after exponentiating
+  tmp4 <- tmp3/max(tmp3)
+
+  ret <- rw_orig(ret, rw = tmp4, conf.index=gauge_conf_list, max_value = max(tmp3))
+
+}
 #' @title reader for Nissa text format correlation functions
 #' @param file_basenames_to_read Character vector of file names without the
 #'                              smearing combination suffixes (such as 'll', 'ls', 'sl', 'ss')