Merge pull request #58 from numericalEFT/zyli

Add Parameter `thermal_ratio` to adjust the thermalization steps for the MCMC method.

Author: peter0627ustc

src/main.jl

@@ -7,6 +7,7 @@
         block=16, 
         measure::Union{Nothing,Function}=nothing,
         measurefreq::Int=1,
+        nburnin::Int=100,
         inplace::Bool=false,
         adapt=true,
         gamma=1.0, 
@@ -35,6 +36,7 @@ Calculate the integrals, collect statistics, and return a `Result` struct contai
     - For `solver = :vegas` or `:vegasmc`, the function signature should be `measure(var, obs, relative_weights, config)`. Here, `obs` is a vector of observable values for each component of the integrand and `relative_weights` are the weights calculated from the integrand multiplied by the probability of the corresponding variables. 
     - For `solver = :mcmc`, the signature should be `measure(idx, var, obs, relative_weight, config)`, where `obs` is the observable vector and `relative_weight` is the weight calculated from the `idx`-th integrand multiplied by the probability of the variables.
 - `measurefreq`: How often the measurement function is called (default: `1`).
+- `nburnin` : Tha thermalization steps for MCMC method
 - `inplace`: Whether to use the inplace version of the integrand. Default is `false`, which is more convenient for integrand with a few return values but may cause type instability. Only useful for the :vegas and :vegasmc solver.
 - `adapt`: Whether to adapt the grid and the reweight factor (default: `true`).
 - `gamma`: Learning rate of the reweight factor after each iteration (default: `1.0`).
@@ -80,6 +82,7 @@ function integrate(integrand::Function;
     ignore::Int=adapt ? 1 : 0, #ignore the first `ignore` iterations in average
     measure::Union{Nothing,Function}=nothing,
     measurefreq::Int=1,
+    nburnin::Int = 100,
     inplace::Bool=false, # whether to use the inplace version of the integrand
     parallel::Symbol=:nothread, # :thread or :nothread
     print=-1, printio=stdout, timer=[],
@@ -151,13 +154,13 @@ function integrate(integrand::Function;
                 Threads.@threads for _ in 1:block/MCUtility.mpi_nprocs()
                     _block!(configs, obsSum, obsSquaredSum, summedConfig, solver, progress,
                         integrand, nevalperblock, print, timer, debug,
-                        measure, measurefreq, inplace, parallel)
+                        measure, measurefreq, nburnin, inplace, parallel)
                 end
             else
                 for _ in 1:block/MCUtility.mpi_nprocs()
                     _block!(configs, obsSum, obsSquaredSum, summedConfig, solver, progress,
                         integrand, nevalperblock, print, timer, debug,
-                        measure, measurefreq, inplace, parallel)
+                        measure, measurefreq, nburnin, inplace, parallel)
                 end
             end
         end
@@ -233,7 +236,7 @@ end
 function _block!(configs, obsSum, obsSquaredSum, summedConfig,
     solver, progress,
     integrand::Function, nevalperblock, print, timer, debug::Bool,
-    measure::Union{Nothing,Function}, measurefreq, inplace, parallel)
+    measure::Union{Nothing,Function}, measurefreq, nburnin, inplace, parallel)
 
     rank = MCUtility.threadid(parallel)
     # println(rank)
@@ -249,7 +252,7 @@ function _block!(configs, obsSum, obsSquaredSum, summedConfig,
             measure=measure, measurefreq=measurefreq, inplace=inplace)
     elseif solver == :mcmc
         MCMC.montecarlo(config_n, integrand, nevalperblock, print, timer, debug;
-            measure=measure, measurefreq=measurefreq)
+            measure=measure, measurefreq=measurefreq, nburnin = nburnin)
     else
         error("Solver $solver is not supported!")
     end

src/mcmc/montecarlo.jl

@@ -73,7 +73,8 @@ function montecarlo(config::Configuration{N,V,P,O,T}, integrand::Function, neval
     verbose=0, timer=[], debug=false;
     measurefreq::Int=1,
     measure::Union{Nothing,Function}=nothing,
-    idx::Int=1 # the integral to start with
+    idx::Int=1, # the integral to start with
+    nburnin::Int=100
 ) where {N,V,P,O,T}
 
     @assert measurefreq > 0
@@ -130,7 +131,7 @@ function montecarlo(config::Configuration{N,V,P,O,T}, integrand::Function, neval
     # end
     startTime = time()
 
-    for i = 1:neval
+    for i = 1:(neval+nburnin) 
         # config.neval += 1
         config.visited[state.curr] += 1
         _update = rand(config.rng, updates) # randomly select an update
@@ -140,7 +141,7 @@ function montecarlo(config::Configuration{N,V,P,O,T}, integrand::Function, neval
         if debug && (isfinite(state.probability) == false)
             @warn("integrand probability = $(state.probability) is not finite at step $(config.neval)")
         end
-        if i % measurefreq == 0 && i >= neval / 100
+        if i % measurefreq == 0 && i >= nburnin
 
             ######## accumulate variable #################
             if state.curr != config.norm