Merge pull request #10 from numericalEFT/dev

add offset to variables

docs/make.jl

@@ -1,32 +1,33 @@
 using MCIntegration
 using Documenter
 
-DocMeta.setdocmeta!(MCIntegration, :DocTestSetup, :(using MCIntegration); recursive=true)
+DocMeta.setdocmeta!(MCIntegration, :DocTestSetup, :(using MCIntegration); recursive = true)
 
 makedocs(;
-    modules=[MCIntegration],
-    authors="Kun Chen, Xiansheng Cai",
-    repo="https://github.com/numericaleft/MCIntegration.jl/blob/{commit}{path}#{line}",
-    sitename="MCIntegration.jl",
-    format=Documenter.HTML(;
-        prettyurls=get(ENV, "CI", "false") == "true",
-        canonical="https://numericaleft.github.io/MCIntegration.jl",
-        assets=String[],
+    modules = [MCIntegration],
+    authors = "Kun Chen, Xiansheng Cai",
+    repo = "https://github.com/numericaleft/MCIntegration.jl/blob/{commit}{path}#{line}",
+    sitename = "MCIntegration.jl",
+    format = Documenter.HTML(;
+        prettyurls = get(ENV, "CI", "false") == "true",
+        canonical = "https://numericaleft.github.io/MCIntegration.jl",
+        assets = String[]
     ),
-    pages=[
+    pages = [
         "Home" => "index.md",
         "Manual" => Any[
             "man/important_sampling.md"
         ],
         "Library" => Any[
-                map(s -> "lib/$(s)", sort(readdir(joinpath(@__DIR__, "src/lib"))))
-                # "Internals" => map(s -> "lib/$(s)", sort(readdir(joinpath(@__DIR__, "src/lib"))))
+            "lib/montecarlo.md",
+        # map(s -> "lib/$(s)", sort(readdir(joinpath(@__DIR__, "src/lib"))))
+        # "Internals" => map(s -> "lib/$(s)", sort(readdir(joinpath(@__DIR__, "src/lib"))))
         ]
-    ],
+    ]
 )
 
 deploydocs(;
-    repo="github.com/numericalEFT/MCIntegration.jl",
+    repo = "github.com/numericalEFT/MCIntegration.jl"
 )
 
 # using Lehmann

example/bubble.jl

@@ -1,11 +1,22 @@
 # This example demonstrated how to calculate the bubble diagram of free electrons using the Monte Carlo module
 
-using QuantumStatistics, LinearAlgebra, Random, Printf, BenchmarkTools, InteractiveUtils, Parameters
+using LinearAlgebra, Random, Printf, BenchmarkTools, InteractiveUtils, Parameters
+using ElectronGas
+using StaticArrays
+using Lehmann
+using MCIntegration
 # using ProfileView
 
-const Steps = 1e8
+const Steps = 1e6
 
-include("parameter.jl")
+# include("parameter.jl")
+beta = 25.0
+rs = 1.0
+const basic = Parameter.rydbergUnit(1 / beta, rs, 3)
+const β = basic.β
+const kF = basic.kF
+const me = basic.me
+const spin = basic.spin
 
 @with_kw struct Para
     n::Int = 0 # external Matsubara frequency
@@ -21,17 +32,17 @@ function integrand(config)
 
     T, K, Ext = config.var[1], config.var[2], config.var[3]
     k = K[1]
-    Tin, Tout = 0.0, T[1] 
+    Tin, Tout = T[1], T[2]
     extidx = Ext[1]
     q = para.extQ[extidx] # external momentum
     kq = k + q
-    τ = (Tout - Tin) / β
-    ω1 = (dot(k, k) - kF^2) / (2me) * β
-    g1 = Spectral.kernelFermiT(τ, ω1)
-    ω2 = (dot(kq, kq) - kF^2) / (2me) * β
-    g2 = Spectral.kernelFermiT(-τ, ω2)
+    τ = (Tout - Tin)
+    ω1 = (dot(k, k) - kF^2) / (2me)
+    g1 = Spectral.kernelFermiT(τ, ω1, β)
+    ω2 = (dot(kq, kq) - kF^2) / (2me)
+    g2 = Spectral.kernelFermiT(-τ, ω2, β)
     phase = 1.0 / (2π)^3
-    return g1 * g2 * spin * phase * cos(2π * para.n * τ)
+    return g1 * g2 * spin * phase * cos(2π * para.n * τ / β) / β
 end
 
 function measure(config)
@@ -45,17 +56,17 @@ end
 function run(steps)
 
     para = Para()
-    @unpack extQ, Qsize = para 
+    @unpack extQ, Qsize = para
 
-    T = MonteCarlo.Tau(β, β / 2.0)
-    K = MonteCarlo.FermiK(3, kF, 0.2 * kF, 10.0 * kF)
-    Ext = MonteCarlo.Discrete(1, length(extQ)) # external variable is specified
+    T = MCIntegration.Tau(β, β / 2.0)
+    K = MCIntegration.FermiK(3, kF, 0.2 * kF, 10.0 * kF)
+    Ext = MCIntegration.Discrete(1, length(extQ)) # external variable is specified
 
-    dof = [[1, 1, 1],] # degrees of freedom of the normalization diagram and the bubble
+    dof = [[2, 1, 1],] # degrees of freedom of the normalization diagram and the bubble
     obs = zeros(Float64, Qsize) # observable for the normalization diagram and the bubble
 
-    config = MonteCarlo.Configuration(steps, (T, K, Ext), dof, obs; para=para)
-    avg, std = MonteCarlo.sample(config, integrand, measure; print=0, Nblock=16)
+    config = MCIntegration.Configuration(steps, (T, K, Ext), dof, obs; para = para)
+    avg, std = MCIntegration.sample(config, integrand, measure; print = 0, Nblock = 16)
     # @profview MonteCarlo.sample(config, integrand, measure; print=0, Nblock=1)
     # sleep(100)
 
@@ -64,10 +75,9 @@ function run(steps)
 
         for (idx, q) in enumerate(extQ)
             q = q[1]
-            p, err = TwoPoint.LindhardΩnFiniteTemperature(dim, q, n, kF, β, me, spin)
-            @printf("%10.6f  %10.6f ± %10.6f  %10.6f ± %10.6f\n", q / kF, avg[idx], std[idx], p, err)
+            p = Polarization.Polarization0_ZeroTemp(q, para.n, basic) * spin
+            @printf("%10.6f  %10.6f ± %10.6f  %10.6f\n", q / basic.kF, avg[idx], std[idx], p)
         end
-
     end
 end
 

src/montecarlo.jl

@@ -47,9 +47,10 @@ sample(config::Configuration, integrand::Function, measure::Function; Nblock=16,
 
 - `reweight = config.totalStep/10`: the MC steps before reweighting the integrands. Set to -1 if reweighting is not wanted.
 """
-function sample(config::Configuration, integrand::Function, measure::Function; Nblock = 16, print = 0, printio = stdout, save = 0, saveio = nothing, timer = [], reweight = config.totalStep / 10)
+function sample(config::Configuration, integrand::Function, measure::Function;
+    Nblock = 16, print = 0, printio = stdout, save = 0, saveio = nothing, timer = [], reweight = config.totalStep / 10)
 
-    println(reweight)
+    # println(reweight)
 
     ############ initialized timer ####################################
     if print > 0
@@ -100,14 +101,8 @@ function sample(config::Configuration, integrand::Function, measure::Function; N
     end
 
     #################### collect statistics  ####################################
-    if typeof(obsSum) <: AbstractArray
-        MPI.Reduce!(obsSum, MPI.SUM, root, comm) # root node gets the sum of observables from all blocks
-        MPI.Reduce!(obsSquaredSum, MPI.SUM, root, comm) # root node gets the squared sum of observables from all blocks
-    else
-        result = [obsSum, obsSquaredSum]  # MPI.Reduce works for array only
-        MPI.Reduce!(result, MPI.SUM, root, comm) # root node gets the sum of observables from all blocks
-        obsSum, obsSquaredSum = result
-    end
+    MPIreduce(obsSum)
+    MPIreduce(obsSquaredSum)
     summary = reduceStat(summary, root, comm) # root node gets the summed MC information
 
     if MPI.Comm_rank(comm) == root
@@ -178,8 +173,19 @@ function doReweight(config)
     for (vi, v) in enumerate(config.visited)
         if v > 1000
             config.reweight[vi] *= avgstep / v
+            if config.reweight[vi] < 1e-10
+                config.reweight[vi] = 1e-10
+            end
         end
     end
+    # renoormalize all reweight to be (0.0, 1.0)
+    config.reweight .= config.reweight ./ sum(config.reweight)
+    # dample reweight factor to avoid rapid, destabilizing changes
+    # reweight factor close to 1.0 will not be changed much
+    # reweight factor close to zero will be amplified significantly
+    # Check Eq. (19) of https://arxiv.org/pdf/2009.05112.pdf for more detail
+    α = 2.0
+    config.reweight = @. ((1 - config.reweight) / log(1 / config.reweight))^α
 end
 
 function printSummary(summary, neighbor, var)

src/sampler.jl

@@ -65,6 +65,7 @@ Propose to generate new Fermi K in [Kf-δK, Kf+δK)
 - `newK`:  vector of dimension of d=2 or 3
 """
 function create!(K::FermiK{D}, idx::Int, config) where {D}
+    @assert idx > K.offset
     (idx >= length(K.data) - 1) && error("$idx overflow!")
     rng = config.rng
     ############ Simple Way ########################
@@ -81,12 +82,17 @@ function create!(K::FermiK{D}, idx::Int, config) where {D}
     if D == 3 # dimension 3
         θ = π * rand(rng)
         # newK .= Kamp .* Mom(cos(ϕ) * sin(θ), sin(ϕ) * sin(θ), cos(θ))
-        K[idx] = @SVector [Kamp * cos(ϕ) * sin(θ), Kamp * sin(ϕ) * sin(θ), Kamp * cos(θ)]
+        # K[idx] = @SVector [Kamp * cos(ϕ) * sin(θ), Kamp * sin(ϕ) * sin(θ), Kamp * cos(θ)]
+        K.data[1, idx] = Kamp * cos(ϕ) * sin(θ)
+        K.data[2, idx] = Kamp * sin(ϕ) * sin(θ)
+        K.data[3, idx] = Kamp * cos(θ)
         return 2 * K.δk * 2π * π * (sin(θ) * Kamp^2)
         # prop density of KAmp in [Kf-dK, Kf+dK), prop density of Phi
         # prop density of Theta, Jacobian
     else  # DIM==2
-        K[idx] = @SVector [Kamp * cos(ϕ), Kamp * sin(ϕ)]
+        # K[idx] = @SVector [Kamp * cos(ϕ), Kamp * sin(ϕ)]
+        K.data[1, idx] = Kamp * cos(ϕ)
+        K.data[2, idx] = Kamp * sin(ϕ)
         return 2 * K.δk * 2π * Kamp
         # prop density of KAmp in [Kf-dK, Kf+dK), prop density of Phi, Jacobian
     end
@@ -102,6 +108,7 @@ Propose to remove an existing Fermi K in [Kf-δK, Kf+δK)
 - `oldK`:  vector of dimension of d=2 or 3
 """
 function remove!(K::FermiK{D}, idx::Int, config) where {D}
+    @assert idx > K.offset
     (idx >= length(K.data) - 1) && error("$idx overflow!")
     ############## Simple Way #########################
     # for i in 1:DIM
@@ -135,6 +142,7 @@ removeRollback!(K::FermiK{D}, idx::Int, config) where {D} = nothing
 Propose to shift oldK to newK. Work for generic momentum vector
 """
 function shift!(K::FermiK{D}, idx::Int, config) where {D}
+    @assert idx > K.offset
     (idx >= length(K.data) - 1) && error("$idx overflow!")
     K[end] = K[idx]  # save current K
 
@@ -151,19 +159,29 @@ function shift!(K::FermiK{D}, idx::Int, config) where {D}
             # sample uniformly on sphere, check http://corysimon.github.io/articles/uniformdistn-on-sphere/ 
             θ = acos(1 - 2 * rand(rng))
             Kamp = sqrt(K[idx][1]^2 + K[idx][2]^2 + K[idx][3]^2)
-            K[idx] = @SVector [Kamp * cos(ϕ) * sin(θ), Kamp * sin(ϕ) * sin(θ), Kamp * cos(θ)]
+            # K[idx] = @SVector [Kamp * cos(ϕ) * sin(θ), Kamp * sin(ϕ) * sin(θ), Kamp * cos(θ)]
+            K.data[1, idx] = Kamp * cos(ϕ) * sin(θ)
+            K.data[2, idx] = Kamp * sin(ϕ) * sin(θ)
+            K.data[3, idx] = Kamp * cos(θ)
             return 1.0
         else # D=2
             Kamp = sqrt(K[idx][1]^2 + K[idx][2]^2)
-            K = @SVector [Kamp * cos(ϕ), Kamp * sin(ϕ)]
+            # K = @SVector [Kamp * cos(ϕ), Kamp * sin(ϕ)]
+            K.data[1, idx] = Kamp * cos(ϕ)
+            K.data[2, idx] = Kamp * sin(ϕ)
             return 1.0
         end
     else
         Kc, dk = K[idx], K.δk
         if (D == 3)
-            K[idx] = @SVector [Kc[1] + (rand(rng) - 0.5) * dk, Kc[2] + (rand(rng) - 0.5) * dk, Kc[3] + (rand(rng) - 0.5) * dk]
+            # K[idx] = @SVector [Kc[1] + (rand(rng) - 0.5) * dk, Kc[2] + (rand(rng) - 0.5) * dk, Kc[3] + (rand(rng) - 0.5) * dk]
+            K.data[1, idx] = Kc[1] + (rand(rng) - 0.5) * dk
+            K.data[2, idx] = Kc[2] + (rand(rng) - 0.5) * dk
+            K.data[3, idx] = Kc[3] + (rand(rng) - 0.5) * dk
         else # D=2
-            K[idx] = @SVector [Kc[1] + (rand(rng) - 0.5) * dk, Kc[2] + (rand(rng) - 0.5) * dk]
+            # K[idx] = @SVector [Kc[1] + (rand(rng) - 0.5) * dk, Kc[2] + (rand(rng) - 0.5) * dk]
+            K.data[1, idx] = Kc[1] + (rand(rng) - 0.5) * dk
+            K.data[2, idx] = Kc[2] + (rand(rng) - 0.5) * dk
         end
         # K[idx] += (rand(rng, D) .- 0.5) .* K.δk
         return 1.0

src/statistics.jl

@@ -8,7 +8,7 @@ mutable struct SummaryStat
     accept::Array{Float64,3}
 end
 
-function addStat(config, summary=nothing)
+function addStat(config, summary = nothing)
     if isnothing(summary)
         return SummaryStat(config.step, config.totalStep, config.visited, config.reweight, config.propose, config.accept)
     else
@@ -37,4 +37,23 @@ function reduceStat(summary, root, comm)
     else
         return summary
     end
+end
+
+function MPIreduce(data)
+    comm = MPI.COMM_WORLD
+    Nworker = MPI.Comm_size(comm)  # number of MPI workers
+    rank = MPI.Comm_rank(comm)  # rank of current MPI worker
+    root = 0 # rank of the root worker
+
+    if Nworker == 1 #no parallelization
+        return data
+    end
+    if typeof(data) <: AbstractArray
+        MPI.Reduce!(data, MPI.SUM, root, comm) # root node gets the sum of observables from all blocks
+        return data
+    else
+        result = [data,]  # MPI.Reduce works for array only
+        MPI.Reduce!(result, MPI.SUM, root, comm) # root node gets the sum of observables from all blocks
+        return result[1]
+    end
 end

src/updates.jl

@@ -14,13 +14,14 @@ function changeIntegrand(config, integrand)
 
     # create/remove variables if there are more/less degrees of freedom
     for vi = 1:length(config.var)
+        offset = config.var[vi].offset
         if (currdof[vi] < newdof[vi]) # more degrees of freedom
             for pos = currdof[vi]+1:newdof[vi]
-                prop *= create!(config.var[vi], pos, config)
+                prop *= create!(config.var[vi], pos + offset, config)
             end
         elseif (currdof[vi] > newdof[vi]) # less degrees of freedom
             for pos = newdof[vi]+1:currdof[vi]
-                prop *= remove!(config.var[vi], pos, config)
+                prop *= remove!(config.var[vi], pos + offset, config)
             end
         end
     end
@@ -44,13 +45,14 @@ function changeIntegrand(config, integrand)
 
         ############ Redo changes to config.var #############
         for vi = 1:length(config.var)
+            offset = config.var[vi].offset
             if (currdof[vi] < newdof[vi]) # more degrees of freedom
                 for pos = currdof[vi]+1:newdof[vi]
-                    createRollback!(config.var[vi], pos, config)
+                    createRollback!(config.var[vi], pos + offset, config)
                 end
             elseif (currdof[vi] > newdof[vi]) # less degrees of freedom
                 for pos = newdof[vi]+1:currdof[vi]
-                    removeRollback!(config.var[vi], pos, config)
+                    removeRollback!(config.var[vi], pos + offset, config)
                 end
             end
         end
@@ -68,7 +70,7 @@ function changeVariable(config, integrand)
     vi = rand(config.rng, 1:length(currdof)) # update the variable type of the index vi
     var = config.var[vi]
     (currdof[vi] <= 0) && return # return if the var has zero degree of freedom
-    idx = rand(config.rng, 1:currdof[vi]) # randomly choose one var to update
+    idx = var.offset + rand(config.rng, 1:currdof[vi]) # randomly choose one var to update
 
     # oldvar = copy(var[idx])
     currAbsWeight = config.absWeight