Update mcmc with thermal_ratio::Float64 for thermalization

.github/workflows/CI.yml

@@ -1,89 +1,82 @@
 name: CI
 on:
-  - push
-  - pull_request
+    push:
+        branches:
+            - main
+            - master
+    pull_request:
+
 jobs:
-  test:
-    name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }}
-    runs-on: ${{ matrix.os }}
-    strategy:
-      fail-fast: false
-      matrix:
-        version:
-          - "1.6"
-          - "nightly"
-        os:
-          - ubuntu-latest
-        arch:
-          - x64
-          - x86
-    steps:
-      - uses: actions/checkout@v2
-      - uses: julia-actions/setup-julia@v1
-        with:
-          version: ${{ matrix.version }}
-          arch: ${{ matrix.arch }}
-      - uses: actions/cache@v1
-        env:
-          cache-name: cache-artifacts
-        with:
-          path: ~/.julia/artifacts
-          key: ${{ runner.os }}-test-${{ env.cache-name }}-${{ hashFiles('**/Project.toml') }}
-          restore-keys: |
-            ${{ runner.os }}-test-${{ env.cache-name }}-
-            ${{ runner.os }}-test-
-            ${{ runner.os }}-
-      - uses: julia-actions/julia-buildpkg@v1
-      - name: Install MPI dependencies
-        shell: bash
-        run: |
-          julia -e '
-            using Pkg; Pkg.add("MPI"); using MPI; MPI.install_mpiexecjl()
-          '
+    test:
+        name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }}
+        runs-on: ${{ matrix.os }}
+        strategy:
+            fail-fast: false
+            matrix:
+                version:
+                    - "1.6"
+                    - "nightly"
+                os:
+                    - ubuntu-latest
+                arch:
+                    - x64
+                    - x86
+        steps:
+            - uses: actions/checkout@v4
+            - uses: julia-actions/setup-julia@v2
+              with:
+                  version: ${{ matrix.version }}
+                  arch: ${{ matrix.arch }}
+            - uses: actions/cache@v4
+              env:
+                  cache-name: cache-artifacts
+              with:
+                  path: ~/.julia/artifacts
+                  key: ${{ runner.os }}-test-${{ env.cache-name }}-${{ hashFiles('**/Project.toml') }}
+                  restore-keys: |
+                      ${{ runner.os }}-test-${{ env.cache-name }}-
+                      ${{ runner.os }}-test-
+                      ${{ runner.os }}-
+
+            - uses: julia-actions/julia-buildpkg@v1
+
+            - name: Install MPI dependencies
+              shell: bash
+              run: |
+                  julia -e '
+                    using Pkg; Pkg.add("MPI"); using MPI; MPI.install_mpiexecjl()
+                  '
 
-      - uses: julia-actions/julia-runtest@v1
-        env:
-          JULIA_NUM_THREADS: 4
-      # - name: Execute MPI-parallel tests
-      #   run: |
-      #     julia --project -e '
-      #       using Pkg; Pkg.build(); Pkg.precompile()
-      #       Pkg.add("MPI"); using MPI; MPI.install_mpiexecjl()
-      #       Pkg.test(; test_args=["quick"])
-      #     '
-      #     $HOME/.julia/bin/mpiexecjl -np 8 julia --check-bounds=yes --depwarn=yes --project --color=yes -e 'using Pkg; Pkg.test(coverage=true)'
-      # if: ${{ matrix.payload == 'mpi' }}
-      # continue-on-error: ${{ matrix.version == 'nightly' }}
+            - uses: julia-actions/julia-runtest@v1
+              env:
+                  JULIA_NUM_THREADS: 4
 
-      - uses: julia-actions/julia-processcoverage@v1
-      - uses: codecov/codecov-action@v2
-        with:
-          file: lcov.info
+            - uses: julia-actions/julia-processcoverage@v1
+            - uses: codecov/codecov-action@v5
+              with:
+                  file: lcov.info
+                  token: ${{ secrets.CODECOV_TOKEN }}
 
-      - uses: julia-actions/julia-processcoverage@v1
-      - uses: codecov/codecov-action@v1
-        with:
-          file: lcov.info
-  docs:
-    name: Documentation
-    runs-on: ubuntu-latest
-    steps:
-      - uses: actions/checkout@v2
-      - uses: julia-actions/setup-julia@v1
-        with:
-          version: "1"
-      - run: |
-          julia --project=docs -e '
-            using Pkg
-            Pkg.develop(PackageSpec(path=pwd()))
-            Pkg.instantiate()'
-      - run: |
-          julia --project=docs -e '
-            using Documenter: DocMeta, doctest
-            using MCIntegration
-            DocMeta.setdocmeta!(MCIntegration, :DocTestSetup, :(using MCIntegration); recursive=true)
-            doctest(MCIntegration)'
-      - run: julia --project=docs docs/make.jl
-        env:
-          GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
-          DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }}
+    docs:
+        name: Documentation
+        runs-on: ubuntu-latest
+        steps:
+            - uses: actions/checkout@v4
+            - uses: julia-actions/setup-julia@v2
+              with:
+                  version: "1"
+            - run: |
+                  julia --project=docs -e '
+                    using Pkg
+                    Pkg.develop(PackageSpec(path=pwd()))
+                    Pkg.instantiate()'
+            - run: |
+                  julia --project=docs -e '
+                    using Documenter: DocMeta, doctest
+                    using MCIntegration
+                    DocMeta.setdocmeta!(MCIntegration, :DocTestSetup, :(using MCIntegration); recursive=true)
+                    doctest(MCIntegration)'
+            - run: julia --project=docs docs/make.jl
+              env:
+                  GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
+                  DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }}

src/configuration.jl

@@ -255,8 +255,8 @@ function Configuration(;
 
     # propose and accept shape: number of updates X integrand number X max(integrand number, variable number)
     # the last index will waste some memory, but the dimension is small anyway
-    propose = zeros(Float64, (2, Nd, max(Nd, Nv))) .+ 1.0e-8 # add a small initial value to avoid Inf when inverted
-    accept = zeros(Float64, (2, Nd, max(Nd, Nv)))
+    propose = zeros(Float64, (3, Nd, max(Nd, Nv))) .+ 1.0e-8 # add a small initial value to avoid Inf when inverted
+    accept = zeros(Float64, (3, Nd, max(Nd, Nv)))
 
     return Configuration{Nd - 1,typeof(var),typeof(userdata),typeof(obs),type}(
         seed, MersenneTwister(seed), var, userdata,  # static parameters
@@ -494,6 +494,34 @@ function report(config::Configuration, total_neval=nothing)
         end
     end
     println(bar)
+
+    println(yellow(@sprintf("%-20s %12s %12s %12s", "SwapVariable", "Proposed", "Accepted", "Ratio  ")))
+    for idx = 1:Nd-1 # normalization diagram don't have variable to change
+        for (vi, var) in enumerate(var)
+            if var isa Continuous
+                typestr = "Continuous"
+            elseif var isa Discrete
+                typestr = "Discrete"
+            elseif var isa CompositeVar
+                typestr = "Composite"
+            elseif var isa FermiK
+                typestr = "FermiK"
+            else
+                typestr = "$(typeof(var))"
+                # typestr = split(typestr, ".")[end]
+            end
+            @printf(
+                "  %2d / %-11s:   %11.6f%% %11.6f%% %12.6f\n",
+                idx, typestr,
+                propose[3, idx, vi] / neval * 100.0,
+                accept[3, idx, vi] / neval * 100.0,
+                accept[3, idx, vi] / propose[3, idx, vi]
+            )
+            totalproposed += propose[3, idx, vi]
+        end
+    end
+    println(bar)
+
     println(yellow("Integrand            Visited      ReWeight"))
     @printf("  Norm   :     %12i %12.6f\n", visited[end], reweight[end])
     for idx = 1:Nd-1

src/main.jl

@@ -7,7 +7,7 @@
         block=16, 
         measure::Union{Nothing,Function}=nothing,
         measurefreq::Int=1,
-        nburnin::Int=100,
+        thermal_ratio::Float64=0.1,
         inplace::Bool=false,
         adapt=true,
         gamma=1.0, 
@@ -36,7 +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
+- `thermal_ratio` : Tha thermalization ratio in one Marovov chain. Default is `0.1`.
 - `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`).
@@ -82,7 +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,
+    thermal_ratio::Float64=0.1,
     inplace::Bool=false, # whether to use the inplace version of the integrand
     parallel::Symbol=:nothread, # :thread or :nothread
     print=-1, printio=stdout, timer=[],
@@ -154,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, nburnin, inplace, parallel)
+                        measure, measurefreq, thermal_ratio, 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, nburnin, inplace, parallel)
+                        measure, measurefreq, thermal_ratio, inplace, parallel)
                 end
             end
         end
@@ -236,11 +236,17 @@ end
 function _block!(configs, obsSum, obsSquaredSum, summedConfig,
     solver, progress,
     integrand::Function, nevalperblock, print, timer, debug::Bool,
-    measure::Union{Nothing,Function}, measurefreq, nburnin, inplace, parallel)
+    measure::Union{Nothing,Function}, measurefreq, thermal_ratio, inplace, parallel)
 
     rank = MCUtility.threadid(parallel)
     # println(rank)
 
+    # Ensure rank is within bounds of configs array
+    # When Threads.nthreads() = 1, @threads can still create tasks with threadid() > 1
+    if rank > length(configs)
+        rank = 1  # Fall back to first config when thread ID exceeds array bounds
+    end
+
     config_n = configs[rank] # configuration for the worker with thread id `rank`
     clearStatistics!(config_n) # reset statistics
 
@@ -252,7 +258,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, nburnin = nburnin)
+            measure=measure, measurefreq=measurefreq, thermal_ratio=thermal_ratio)
     else
         error("Solver $solver is not supported!")
     end

src/mcmc/montecarlo.jl

@@ -74,7 +74,7 @@ function montecarlo(config::Configuration{N,V,P,O,T}, integrand::Function, neval
     measurefreq::Int=1,
     measure::Union{Nothing,Function}=nothing,
     idx::Int=1, # the integral to start with
-    nburnin::Int=100
+    thermal_ratio::Float64=0.1
 ) where {N,V,P,O,T}
 
     @assert measurefreq > 0
@@ -131,7 +131,8 @@ function montecarlo(config::Configuration{N,V,P,O,T}, integrand::Function, neval
     # end
     startTime = time()
 
-    for i = 1:(neval+nburnin) 
+    nburnin = Int(floor(neval * thermal_ratio))
+    for i = 1:(neval+nburnin)
         # config.neval += 1
         config.visited[state.curr] += 1
         _update = rand(config.rng, updates) # randomly select an update