Dev

README.md

@@ -27,49 +27,41 @@ function integrand(config)
     end
 end
 
-# Define how to measure the observable
-function measure(config)
-    factor = 1.0 / config.reweight[config.curr]
-    weight = integrand(config)
-    config.observable[config.curr] += weight / abs(weight) * factor #note that config.observable is an array with two elements as discussed below
-end
-
-# MC step of each block
-const blockStep = 1e6
+# MC step of each iteration, and the iteration number. After each iteraction, the program will try to improve the important sampling
+const neval, niter = 1e4, 10
 
-# Define the types variables, the first argument sets the range, the second argument gives the largest change to the variable in one MC update. see the section [variable](#variable) for more details.
-T = MCIntegration.Continuous([0.0, 1.0], 0.5)
+# Define the types variables, the first two arguments set the boundary. see the section [variable](#variable) for more details.
+T = Continuous(0.0, 1.0)
 
 # Define how many (degrees of freedom) variables of each type. 
 # For example, [[n1, n2], [m1, m2], ...] means the first integral involves n1 varibales of type 1, and n2 variables of type2, while the second integral involves m1 variables of type 1 and m2 variables of type 2. 
 dof = [[2,], [3,]]
 
-# Define the container for the observable. It must be a number or an array-like object. In this case, the observable has two elements, corresponds to the results for the two integrals. 
-obs = [0.0, 0.0]
-
 # Define the configuration struct which is container of all kinds of internal data for MC,
-# the second argument is a tuple listing all types of variables, one then specify the degrees of freedom of each variable type in the third argument.  
-config = MCIntegration.Configuration(blockStep, (T,), dof, obs)
+# the first argument is a tuple listing all types of variables, one then specify the degrees of freedom of each variable type in the second argument.  
+config = Configuration((T,), dof)
 
-# perform MC integration. Nblock is the number of independent blocks to estimate the error bar. In MPI mode, the blocks will be sent to different workers. Set "print=n" to control the level of information to print.
-avg, err = MCIntegration.sample(config, integrand, measure; Nblock = 64, print = 1)
+# perform MC integration. Set print>=0 to print more information.
+result = sample(config, integrand; neval=neval, niter=niter, print=-1)
 
-#avg, err are the same shape as obs. In MPI mode, only the root node return meaningful estimates. All other workers simply return nothing
-if isnothing(avg) == false
+# In MPI mode, only the root node return meaningful estimates. All other workers simply return nothing
+if isnothing(result) == false
+    # MCIntegration.summary(result) # uncomment this line to print the summary of the result
+    avg, err = result.mean, result.stdev
     println("Circle area: $(avg[1]) +- $(err[1]) (exact: $(π / 4.0))")
     println("Sphere volume: $(avg[2]) +- $(err[2]) (exact: $(4.0 * π / 3.0 / 8))")
 end
 ```
 
 # Variables
 
-This package defines some common types of variables. Internally, each variable type holds a vector of variables (which is the field named `data`). The actual number of variables in this vector is called the degrees of freedom (dof). Note that different integral may share the same variable types, but have different degrees of freedom. In the above code example, the integral for the circle area and the sphere volume both involve the variable type `Tau`. The former has dof=2, while the latter has dof=3. 
+This package defines some common types of variables. Internally, each variable type holds a vector of variables (which is the field named `data`). The actual number of variables in this vector is called the degrees of freedom (dof). Note that different integral may share the same variable types, but have different degrees of freedom. In the above code example, the integral for the circle area and the sphere volume both involve the variable type `Continuous`. The former has dof=2, while the latter has dof=3. 
 
-Here we briefly list some of the common variables types
+Here we list some of the common variables types
 
-- Continous([start, end], length scale): continuous real-valued variables with specified range and a length scale. The length scale controls the change of the variable in one MC update. A reasonable estimate of the length scale improves the MC efficiency.
+- Continous(start, end): continuous real-valued variables on the domain [start, end). MC will learn the distribution and perform an imporant sampling accordingly.
 
-- Discrete(lower, upper): integer variables in the closed set [lower, upper]. MC will uniformly sample all integers within this range.
+- Discrete(lower, upper): integer variables in the closed set [lower, upper]. MC will learn the distribution and perform an imporant sampling accordingly.
 
 More supported variables types can be found in the [source code](src/variable.jl).
 

src/MCIntegration.jl

@@ -2,6 +2,7 @@ module MCIntegration
 include("utility/utility.jl")
 
 include("montecarlo.jl")
-export montecarlo, Configuration, Diagram, FermiK, BoseK, Tau, TauPair
-
+export Configuration, FermiK, BoseK, Tau, TauPair
+export Continuous, Discrete
+export sample
 end

src/configuration.jl

@@ -15,8 +15,6 @@ mutable struct Configuration
 
  - `para`: user-defined parameter, set to nothing if not needed
 
- - `totalStep`: the total number of updates for this configuration
-
  - `var`: TUPLE of variables, each variable should be derived from the abstract type Variable, see variable.jl for details). Use a tuple rather than a vector improves the performance.
 
  ## integrand properties
@@ -77,35 +75,39 @@ mutable struct Configuration{V,P,O}
     accept::Array{Float64,3} # updates index, integrand index, integrand index 
 
     """
-    Configuration(totalStep, var::V, dof, obs::O; para::P=nothing, state=nothing, reweight=nothing, seed=nothing, neighbor=Vector{Vector{Int}}([])) where {V,P,O}
+    function Configuration(var::V, dof, obs::O=length(dof) == 1 ? 0.0 : zeros(length(dof));
+        para::P=nothing,
+        reweight::Vector{Float64}=ones(length(dof) + 1),
+        seed::Int=rand(Random.RandomDevice(), 1:1000000),
+        neighbor::Union{Vector{Vector{Int}},Vector{Tuple{Int,Int}},Nothing}=nothing
+    ) where {V,P,O}
 
     Create a Configuration struct
 
- # Arguments
+     # Arguments
 
- ## Static parameters
+     ## Static parameters
 
- - `totalStep`: the total number MC steps of each block (one block, one configuration)
-
- - `var`: TUPLE of variables, each variable should be derived from the abstract type Variable, see variable.jl for details). Use a tuple rather than a vector improves the performance.
+     - `var`: TUPLE of variables, each variable should be derived from the abstract type Variable, see variable.jl for details). Use a tuple rather than a vector improves the performance.
 
- - `dof::Vector{Vector{Int}}`: degrees of freedom of each integrand, e.g., [[0, 1], [2, 3]] means the first integrand has zero var#1 and one var#2; while the second integrand has two var#1 and 3 var#2. 
+     - `dof::Vector{Vector{Int}}`: degrees of freedom of each integrand, e.g., [[0, 1], [2, 3]] means the first integrand has zero var#1 and one var#2; while the second integrand has two var#1 and 3 var#2. 
 
- - `obs`: observables that is required to calculate the integrands, will be used in the `measure` function call
-    It is either an array of any type with the common operations like +-*/^ defined. 
+     - `obs`: observables that is required to calculate the integrands, will be used in the `measure` function call.
+        It is either an array of any type with the common operations like +-*/^ defined. 
+        By default, it will be set to 0.0 if there is only one integrand (e.g., length(dof)==1); otherwise, it will be set to zeros(length(dof)).
 
- - `para`: user-defined parameter, set to nothing if not needed
+     - `para`: user-defined parameter, set to nothing if not needed
 
- - `reweight`: reweight factors for each integrands. If not set, then all factors will be initialized with one.
+     - `reweight`: reweight factors for each integrands. If not set, then all factors will be initialized with one.
 
- - `seed`: seed to initialize random numebr generator, also serves as the unique pid of the configuration. If it is nothing, then use RandomDevice() to generate a random seed in [1, 1000_1000]
+     - `seed`: seed to initialize random numebr generator, also serves as the unique pid of the configuration. If it is nothing, then use RandomDevice() to generate a random seed in [1, 1000_1000]
 
-- `neighbor::Vector{Tuple{Int, Int}}` : vector of tuples that defines the neighboring integrands. Two neighboring integrands are directly connected in the Markov chain. 
-    e.g., [(1, 2), (2, 3)] means the integrand 1 and 2 are neighbor, and 2 and 3 are neighbor.  
-   The neighbor vector defines a undirected graph showing how the integrands are connected. Please make sure all integrands are connected.
-   By default, we assume the N integrands are in the increase order, meaning the neighbor will be set to [(N+1, 1), (1, 2), (2, 4), ..., (N-1, N)], where the first N entries are for diagram 1, 2, ..., N and the last entry is for the normalization diagram. Only the first diagram is connected to the normalization diagram.
-   Only highly correlated integrands are not highly correlated should be defined as neighbors. Otherwise, most of the updates between the neighboring integrands will be rejected and wasted.
-"""
+    - `neighbor::Vector{Tuple{Int, Int}}` : vector of tuples that defines the neighboring integrands. Two neighboring integrands are directly connected in the Markov chain. 
+        e.g., [(1, 2), (2, 3)] means the integrand 1 and 2 are neighbor, and 2 and 3 are neighbor.  
+       The neighbor vector defines a undirected graph showing how the integrands are connected. Please make sure all integrands are connected.
+       By default, we assume the N integrands are in the increase order, meaning the neighbor will be set to [(N+1, 1), (1, 2), (2, 4), ..., (N-1, N)], where the first N entries are for diagram 1, 2, ..., N and the last entry is for the normalization diagram. Only the first diagram is connected to the normalization diagram.
+       Only highly correlated integrands are not highly correlated should be defined as neighbors. Otherwise, most of the updates between the neighboring integrands will be rejected and wasted.
+    """
     function Configuration(var::V, dof, obs::O=length(dof) == 1 ? 0.0 : zeros(length(dof));
         para::P=nothing,
         reweight::Vector{Float64}=ones(length(dof) + 1),

src/montecarlo.jl

@@ -33,9 +33,17 @@ sample(config::Configuration, integrand::Function, measure::Function; Nblock=16,
  - `integrand`: function call to evaluate the integrand. It should accept an argument of the type `Configuration`, and return a weight. 
     Internally, MC only samples the absolute value of the weight. Therefore, it is also important to define Main.abs for the weight if its type is user-defined. 
 
-- `measure`: function call to measure. It should accept an argument of the type `Configuration`, then manipulate observables `obs`. 
+- `measure`: function call to measure. It should accept an argument of the type `Configuration`, then manipulate observables `obs`. By default, the function MCIntegration.simple_measure will be used.
 
-- `Nblock`: Number of blocks, each block corresponds to one Configuration struct. The tasks will automatically distributed to multi-process in MPI mode. If the numebr of workers N is larger than Nblock, then Nblock will be set to be N.
+- `neval`: number of evaluations of the integrand per iteration. By default, it is set to 1e4 * length(config.dof).
+
+- `niter`: number of iterations. The reweight factor and the variables will be self-adapted after each iteration. By default, it is set to 10.
+
+- `block`: Number of blocks. Each block will be evaluated by about neval/block times. The results from the blocks will be assumed to be statistically independent, and will be used to estimate the error.
+   The tasks will automatically distributed to multi-process in MPI mode. If the numebr of workers N is larger than block, then block will be set to be N.
+   By default, it is set to 16.
+
+- `alpha`: Learning rate of the reweight factor after each iteraction. Note that alpha <=1, where alpha = 0 means no reweighting.  
 
 - `print`: -1 to not print anything, 0 to print minimal information, >0 to print summary for every `print` seconds
 
@@ -46,20 +54,17 @@ sample(config::Configuration, integrand::Function, measure::Function; Nblock=16,
 - `saveio`: `io` to save
 
 - `timer`: `StopWatch` other than print and save.
-
-- `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=simple_measure;
     neval=1e4 * length(config.dof), # number of evaluations
     niter=10, # number of iterations
     block=16, # number of blocks
-    alpha=0.5, # learning rate
-    beta=1.0, # learning rate
+    alpha=1.0, # learning rate of the reweight factor
     print=0, printio=stdout, save=0, saveio=nothing, timer=[])
 
     # println(reweight)
 
-    if beta > 1.0
+    if alpha > 1.0
         @warn(red("beta should be less than 1.0"))
     end
 
@@ -153,7 +158,7 @@ function sample(config::Configuration, integrand::Function, measure::Function=si
             # average(results)
 
             ################### self-learning ##########################################
-            doReweight!(summedConfig, beta)
+            doReweight!(summedConfig, alpha)
             for var in summedConfig.var
                 train!(var)
             end
@@ -209,7 +214,7 @@ function montecarlo(config::Configuration, integrand::Function, measure::Functio
     end
 
     ########### MC simulation ##################################
-    if (print >= 0)
+    if (print > 0)
         println(green("Seed $(config.seed) Start Simulation ..."))
     end
     startTime = time()
@@ -249,7 +254,7 @@ function montecarlo(config::Configuration, integrand::Function, measure::Functio
         end
     end
 
-    if (print >= 0)
+    if (print > 0)
         println(green("Seed $(config.seed) End Simulation. Cost $(time() - startTime) seconds."))
     end
 
@@ -268,14 +273,14 @@ function simple_measure(config, integrand)
     end
 end
 
-function doReweight!(config, beta)
+function doReweight!(config, alpha)
     avgstep = sum(config.visited)
     for (vi, v) in enumerate(config.visited)
         # if v > 1000
         if v <= 1
-            config.reweight[vi] *= (avgstep)^beta
+            config.reweight[vi] *= (avgstep)^alpha
         else
-            config.reweight[vi] *= (avgstep / v)^beta
+            config.reweight[vi] *= (avgstep / v)^alpha
         end
     end
     # renoormalize all reweight to be (0.0, 1.0)

src/statistics.jl

@@ -1,3 +1,18 @@
+"""
+struct Result{O,C}
+
+    the returned result of the MC integration.
+
+# Members
+
+- `mean`: mean of the MC integration
+- `stdev`: standard deviation of the MC integration
+- `chi2`: chi-square of the MC integration
+- `neval`: number of evaluations of the integrand
+- `dof`: degrees of freedom of the MC integration (number of iterations - 1)
+- `config`: configuration of the MC integration from the last iteration
+- `iteractions`: list of tuples [(data, error, Configuration), ...] from each iteration
+"""
 struct Result{O,C}
     mean::O
     stdev::O
@@ -35,6 +50,14 @@ function tostring(mval, merr; pm="±")
     # return "$val $pm $(round(merr, sigdigits=error_digits))"
 end
 
+"""
+function summary(result::Result, pick::Union{Function,AbstractVector}=obs -> real(first(obs)))
+
+    print the summary of the result. 
+    It will first print the configuration from the last iteration, then print the weighted average and standard deviation of the picked observable from each iteration.
+    The pick function is used to select one of the observable to be printed. The return value of pick function must be a Number.
+
+"""
 function summary(result::Result, pick::Union{Function,AbstractVector}=obs -> real(first(obs)))
     summary(result.config)
 
@@ -81,6 +104,16 @@ function MPIreduce(data)
     end
 end
 
+"""
+
+function average(history, max=length(history))
+
+average the history[1:max]. Return the mean, standard deviation and chi2 of the history.
+
+# Arguments
+- `history`: a list of tuples, such as [(data, error, Configuration), ...]
+- `max`: the number of data to average over
+"""
 function average(history, max=length(history))
     @assert max > 0
     if max == 1

src/utility/utility.jl

@@ -64,6 +64,11 @@ function locate(accumulation::AbstractVector, p::Number)
     # error("p=$p is out of the upper bound $(accumulation[end])")
 end
 
+"""
+function smooth(dist::AbstractVector, factor=6)
+
+    Smooth the distribution by averaging two nearest neighbor. The average ratio is given by 1 : factor : 1 for the elements which are not on the boundary.
+"""
 function smooth(dist::AbstractVector, factor=6)
     if length(dist) <= 1
         return dist
@@ -77,7 +82,17 @@ function smooth(dist::AbstractVector, factor=6)
     return new
 end
 
-function rescale(dist::AbstractVector, alpha=1.0)
+"""
+function rescale(dist::AbstractVector, alpha=1.5)
+
+    rescale the dist array to avoid overreacting to atypically large number.
+    There are three steps:
+    1. dist will be first normalize to [0, 1].
+    2. Then the values that are close to 1.0 will not be changed much, while that close to zero will be amplified to a value controlled by alpha.
+    3. In the end, the rescaled dist array will be normalized to [0, 1].
+    Check Eq. (19) of https://arxiv.org/pdf/2009.05112.pdf for more detail
+"""
+function rescale(dist::AbstractVector, alpha=1.5)
     if length(dist) == 1
         return dist
     end

test/montecarlo.jl

@@ -180,14 +180,11 @@ function TestDiscrete(totalstep)
     end
 
     X = MCIntegration.Discrete(1, 3, adapt=true)
-    dof = [[1,],] # number of X variable of tthe integrand
+    dof = [[1,],] # number of X variable of the integrand
     config = MCIntegration.Configuration((X,), dof)
 
     result = MCIntegration.sample(config, integrand; neval=totalstep, niter=10, block=64, print=-1)
     # println(MCIntegration.summary(result))
-    # println(result.config.var[1].histogram)
-    # println(result.config.var[1].distribution)
-    # println(result.config.var[1].accumulation)
     return result.mean, result.stdev
 end