#
# Project : Gardenia
# Source  : som.jl
# Author  : Li Huang (huangli@caep.cn)
# Status  : Unstable
#
# Last modified: 2024/10/01
#

#=
### *Customized Structs* : *StochOM Solver*
=#

"""
    Box

Rectangle. The field configuration consists of many boxes. They exhibit
various areas (width √ó height). We use the Metropolis important sampling
algorithm to sample them and evaluate their contributions to the spectrum.

### Members
* h -> Height of the box.
* w -> Width of the box.
* c -> Position of the box.
"""
mutable struct Box
    h :: F64
    w :: F64
    c :: F64
end

"""
    StochOMElement

Mutable struct. It is used to record the field configurations, which will
be sampled by Monte Carlo sweeping procedure.

### Members
* C -> Field configuration.
* Œõ -> Contributions of the field configuration to the correlator.
* G -> Reproduced correlator.
* Œî -> Difference between reproduced and raw correlators.
"""
mutable struct StochOMElement
    C :: Vector{Box}
    Œõ :: Array{F64,2}
    G :: Vector{F64}
    Œî :: F64
end

"""
    StochOMContext

Mutable struct. It is used within the StochOM solver only.

### Members
* G·µ•    -> Input data for correlator.
* œÉ¬π    -> Actually 1.0 / œÉ¬π.
* grid  -> Grid for input data.
* mesh  -> Mesh for output spectrum.
* C·µ•    -> It is used to record the field configurations for all attempts.
* Œî·µ•    -> It is used to record the errors for all attempts.
* ùïä·µ•    -> It is used to interpolate the Œõ functions.
"""
mutable struct StochOMContext
    G·µ•   :: Vector{F64}
    œÉ¬π   :: Vector{F64}
    grid :: AbstractGrid
    mesh :: AbstractMesh
    C·µ•   :: Vector{Vector{Box}}
    Œî·µ•   :: Vector{F64}
    ùïä·µ•   :: Vector{CubicSplineInterpolation}
end

#=
### *Global Drivers*
=#

"""
    solve(S::StochOMSolver, rd::RawData)

Solve the analytic continuation problem by the stochastic optimization
method. This solver requires a lot of computational resources to get
reasonable results. It is suitable for both Matsubara and imaginary
time correlators. It is the driver for the StochOM solver.

If the input correlators are bosonic, this solver will return A(œâ) / œâ
via `Aout`, instead of A(œâ). At this time, `Aout` is not compatible with
`Gout`. If the input correlators are fermionic, this solver will return
A(œâ) in `Aout`. Now it is compatible with `Gout`. These behaviors are just
similar to the MaxEnt, StochAC, and StochSK solvers.

Now the StochOM solver supports both continuous and Œ¥-like spectra.

### Arguments
* S -> A StochOMSolver struct.
* rd -> A RawData struct, containing raw data for input correlator.

### Returns
* mesh -> Real frequency mesh, œâ.
* Aout -> Spectral function, A(œâ).
* Gout -> Retarded Green's function, G(œâ).
"""
function solve(S::StochOMSolver, rd::RawData)
    println("[ StochOM ]")
    MC, SC = init(S, rd)

    # Parallel version
    if nworkers() > 1
        #
        println("Using $(nworkers()) workers")
        #
        # Copy configuration dicts
        p1 = deepcopy(PBASE)
        p2 = deepcopy(PStochOM)
        #
        #  Launch the tasks one by one
        ùêπ = Future[]
        for i = 1:nworkers()
            ùëì = @spawnat i + 1 prun(S, p1, p2, MC, SC)
            push!(ùêπ, ùëì)
        end
        #
        # Wait and collect the solutions
        sol = []
        for i = 1:nworkers()
            wait(ùêπ[i])
            push!(sol, fetch(ùêπ[i]))
        end
        #
        # Average the solutions
        Aout = similar(sol[end])
        fill!(Aout, 0.0)
        for i in eachindex(sol)
            @. Aout = Aout + sol[i] / nworkers()
        end
        #
        # Postprocess the solutions
        Gout = last(SC, Aout)
        #
    # Sequential version
    else
        #
        Aout = run(MC, SC)
        Gout = last(SC, Aout)
        #
    end

    return SC.mesh.mesh, Aout, Gout
end

"""
    init(S::StochOMSolver, rd::RawData)

Initialize the StochOM solver and return the StochOMMC and StochOMContext
structs. Please don't call this function directly.

### Arguments
* S -> A StochOMSolver struct.
* rd -> A RawData struct, containing raw data for input correlator.

### Returns
* MC -> A StochOMMC struct.
* SC -> A StochOMContext struct.
"""
function init(S::StochOMSolver, rd::RawData)
    # Prepare input data
    G·µ•, œÉ¬π = init_iodata(S, rd)
    println("Postprocess input data: ", length(œÉ¬π), " points")

    # Prepare grid for input data
    grid = make_grid(rd)
    println("Build grid for input data: ", length(grid), " points")

    # Prepare mesh for output spectrum
    mesh = make_mesh()
    println("Build mesh for spectrum: ", length(mesh), " points")

    # Initialize counters for Monte Carlo engine
    MC = init_mc(S)
    println("Create infrastructure for Monte Carlo sampling")

    # Prepare some key variables
    C·µ•, Œî·µ•, ùïä·µ• = init_context(S, grid)
    SC = StochOMContext(G·µ•, œÉ¬π, grid, mesh, C·µ•, Œî·µ•, ùïä·µ•)
    println("Initialize context for the StochOM solver")

    return MC, SC
end

"""
    run(MC::StochOMMC, SC::StochOMContext)

Perform stochastic optimization simulation, sequential version.

### Arguments
* MC -> A StochOMMC struct.
* SC -> A StochOMContext struct.

### Returns
* Aout -> Spectral function, A(œâ).
"""
function run(MC::StochOMMC, SC::StochOMContext)
    # By default, we should write the analytic continuation results
    # into the external files.
    _fwrite = get_b("fwrite")
    fwrite = isa(_fwrite, Missing) || _fwrite ? true : false

    # Setup essential parameters
    ntry = get_s("ntry")
    nstep = get_s("nstep")

    # Sample and collect data
    println("Start stochastic sampling...")
    for l = 1:ntry
        # Re-initialize the simulation
        SE = init_element(MC, SC)

        # For each attempt, we should perform `nstep √ó N` Monte Carlo
        # updates, where `N` means length of the Markov chain.
        for _ = 1:nstep
            update(MC, SE, SC)
        end

        # Write Monte Carlo statistics
        l % 10 == 0 && fwrite && write_statistics(MC)

        # Accumulate the data
        SC.Œî·µ•[l] = SE.Œî
        SC.C·µ•[l] = deepcopy(SE.C)

        # Show error function for the current attempt
        @printf("try -> %6i (%6i) Œî -> %8.4e \n", l, ntry, SE.Œî)
        flush(stdout)
    end

    # Generate spectral density from Monte Carlo field configuration
    return average(SC)
end

"""
    prun(
        S::StochOMSolver,
        p1::Dict{String,Vector{Any}},
        p2::Dict{String,Vector{Any}},
        MC::StochOMMC,
        SC::StochOMContext
    )

Perform stochastic optimization simulation, parallel version.
The arguments `p1` and `p2` are copies of PBASE and PStochOM, respectively.

### Arguments
* S -> A StochOMSolver struct.
* p1 -> A copy of PBASE.
* p2 -> A copy of PStochOM.
* MC -> A StochOMMC struct.
* SC -> A StochOMContext struct.

### Returns
* Aout -> Spectral function, A(œâ).
"""
function prun(
    S::StochOMSolver,
    p1::Dict{String,Vector{Any}},
    p2::Dict{String,Vector{Any}},
    MC::StochOMMC,
    SC::StochOMContext
    )
    # Revise parameteric dicts
    rev_dict_b(p1)
    rev_dict_s(S, p2)

    # Initialize random number generator again
    MC.rng = MersenneTwister(rand(1:10000) * myid() + 1981)

    # By default, we should write the analytic continuation results
    # into the external files.
    _fwrite = get_b("fwrite")
    fwrite = isa(_fwrite, Missing) || _fwrite ? true : false

    # Setup essential parameters
    ntry = get_s("ntry")
    nstep = get_s("nstep")

    # Sample and collect data
    println("Start stochastic sampling...")
    for l = 1:ntry
        # Re-initialize the simulation
        SE = init_element(MC, SC)

        # For each attempt, we should perform `nstep √ó N` Monte Carlo
        # updates, where `N` means length of the Markov chain.
        for _ = 1:nstep
            update(MC, SE, SC)
        end

        # Write Monte Carlo statistics
        myid() == 2 && l % 10 == 0 && fwrite && write_statistics(MC)

        # Accumulate the data
        SC.Œî·µ•[l] = SE.Œî
        SC.C·µ•[l] = deepcopy(SE.C)

        # Show error function for the current attempt
        @printf("try -> %6i (%6i) Œî -> %8.4e \n", l, ntry, SE.Œî)
        flush(stdout)
    end

    # Generate spectral density from Monte Carlo field configuration
    return average(SC)
end

"""
    average(SC::StochOMContext)

Postprocess the collected results after the stochastic optimization
simulations. It will generate the spectral functions.

### Arguments
* SC -> A StochOMContext struct.

### Returns
* Aom -> Spectral function, A(œâ).
"""
function average(SC::StochOMContext)
    # By default, we should write the analytic continuation results
    # into the external files.
    _fwrite = get_b("fwrite")
    fwrite = isa(_fwrite, Missing) || _fwrite ? true : false

    nmesh = get_b("nmesh")
    ntry  = get_s("ntry")

    # Calculate the median of SC.Œî·µ•
    dev_ave = median(SC.Œî·µ•)

    # Determine the Œ±good parameter, which is used to filter the
    # calculated spectra.
    Œ±good = 1.2
    if count(x -> x < dev_ave / Œ±good, SC.Œî·µ•) == 0
        Œ±good = 1.0
    end

    # Accumulate the final spectrum
    Aom = zeros(F64, nmesh)
    passed = I64[]
    for l = 1:ntry
        # Filter the reasonable spectra
        if SC.Œî·µ•[l] < dev_ave / Œ±good
            # Generate the spectrum, and add it to Aom.
            for w = 1:nmesh
                _omega = SC.mesh[w]
                # Scan all boxes
                for r = 1:length(SC.C·µ•[l])
                    R = SC.C·µ•[l][r]
                    # Yes, this box contributes. The point, _omega, is
                    # covered by this box.
                    if R.c - 0.5 * R.w ‚â§ _omega ‚â§ R.c + 0.5 * R.w
                        Aom[w] = Aom[w] + R.h
                    end
                end
            end
            #
            # Record which spectrum is used
            append!(passed, l)
        end
    end

    # Normalize the spectrum
    Lgood = count(x -> x < dev_ave / Œ±good, SC.Œî·µ•)
    @assert Lgood == length(passed)
    @. Aom = Aom / Lgood
    @printf("Median œá¬≤ : %16.12e\n", dev_ave)
    @printf("Accepted configurations : %5i\n", Lgood)

    # Write indices of selected solutions
    if nworkers() > 1
        myid() == 2 && fwrite && write_passed(passed, dev_ave, Œ±good)
    else
        fwrite && write_passed(passed, dev_ave, Œ±good)
    end

    return Aom
end

"""
    last(SC::StochOMContext, Aout::Vector{F64})

It will process and write the calculated results by the StochOM solver,
including final spectral function and reproduced correlator.

### Arguments
* SC   -> A StochOMContext struct.
* Aout -> Spectral function, A(œâ).

### Returns
* G -> Retarded Green's function, G(œâ).
"""
function last(SC::StochOMContext, Aout::Vector{F64})
    # By default, we should write the analytic continuation results
    # into the external files.
    _fwrite = get_b("fwrite")
    fwrite = isa(_fwrite, Missing) || _fwrite ? true : false

    # Write the spectral function
    fwrite && write_spectrum(SC.mesh, Aout)

    # Reproduce input data and write them
    kernel = make_kernel(SC.mesh, SC.grid)
    G = reprod(SC.mesh, kernel, Aout)
    fwrite && write_backward(SC.grid, G)

    # Calculate full response function on real axis and write them
    if get_b("ktype") == "fermi"
        _G = kramers(SC.mesh, Aout)
    else
        _G = kramers(SC.mesh, Aout .* SC.mesh)
    end
    fwrite && write_complete(SC.mesh, _G)

    return _G
end

#=
### *Core Algorithms*
=#

"""
    update(MC::StochOMMC, SE::StochOMElement, SC::StochOMContext)

Using the Metropolis algorithm to update the field configuration, i.e, a
collection of hundreds of boxes. Be careful, this function only updates
the Monte Carlo configurations (in other words, `SE`). It doesn't record
them. Measurements are done in `run()` and `prun()`. This is the reason
why this function is named as `update()`, instead of `sample()`.

### Arguments
* MC -> A StochOMMC struct. It containts some counters.
* SE -> A StochOMElement struct. It contains Monte Carlo configurations.
* SC -> A StochOMContext struct. It contains grid, mesh, and G·µ•.

### Returns
N/A
"""
function update(MC::StochOMMC, SE::StochOMElement, SC::StochOMContext)
    Tmax = 100 # Length of the Markov chain
    nbox = get_s("nbox")

    ST = deepcopy(SE)

    # The Markov chain is divided into two stages
    T1 = rand(MC.rng, 1:Tmax)
    d1 = rand(MC.rng, F64)
    d2 = 1.0 + rand(MC.rng, F64)
    #
    # The first stage
    for _ = 1:T1
        update_type = rand(MC.rng, 1:7)

        @cswitch update_type begin
            @case 1
                if 1 ‚â§ length(ST.C) ‚â§ nbox - 1
                    try_insert(MC, ST, SC, d1)
                end
                break

            @case 2
                if length(ST.C) ‚â• 2
                    try_remove(MC, ST, SC, d1)
                end
                break

            @case 3
                if length(ST.C) ‚â• 1
                    try_shift(MC, ST, SC, d1)
                end
                break

            @case 4
                if length(ST.C) ‚â• 1
                    try_width(MC, ST, SC, d1)
                end
                break

            @case 5
                if length(ST.C) ‚â• 2
                    try_height(MC, ST, SC, d1)
                end
                break

            @case 6
                if 1 ‚â§ length(ST.C) ‚â§ nbox - 1
                    try_split(MC, ST, SC, d1)
                end
                break

            @case 7
                if length(ST.C) ‚â• 2
                    try_merge(MC, ST, SC, d1)
                end
                break
        end

    end
    #
    # The second stage
    for _ = T1+1:Tmax
        update_type = rand(MC.rng, 1:7)

        @cswitch update_type begin
            @case 1
                if 1 ‚â§ length(ST.C) ‚â§ nbox - 1
                    try_insert(MC, ST, SC, d2)
                end
                break

            @case 2
                if length(ST.C) ‚â• 2
                    try_remove(MC, ST, SC, d2)
                end
                break

            @case 3
                if length(ST.C) ‚â• 1
                    try_shift(MC, ST, SC, d2)
                end
                break

            @case 4
                if length(ST.C) ‚â• 1
                    try_width(MC, ST, SC, d2)
                end
                break

            @case 5
                if length(ST.C) ‚â• 2
                    try_height(MC, ST, SC, d2)
                end
                break

            @case 6
                if 1 ‚â§ length(ST.C) ‚â§ nbox - 1
                    try_split(MC, ST, SC, d2)
                end
                break

            @case 7
                if length(ST.C) ‚â• 2
                    try_merge(MC, ST, SC, d2)
                end
                break
        end
    end

    if ST.Œî < SE.Œî
        SE.C = deepcopy(ST.C)
        SE.Œõ .= ST.Œõ
        SE.G .= ST.G
        SE.Œî  = ST.Œî
    end
end

#=
### *Service Functions*
=#

"""
    init_iodata(S::StochOMSolver, rd::RawData)

Preprocess the input data (`rd`).

### Arguments
* S -> A StochOMSolver struct.
* rd -> A RawData struct, which contains essential input data.

### Returns
* G·µ• -> Input correlator.
* œÉ¬π -> 1.0 / œÉ¬π.

See also: [`RawData`](@ref), [`GreenData`](@ref).
"""
function init_iodata(S::StochOMSolver, rd::RawData)
    G = make_data(rd)
    G·µ• = G.value
    œÉ¬π = 1.0 ./ G.error

    return G·µ•, œÉ¬π
end

"""
    init_mc(S::StochOMSolver)

Try to create a StochOMMC struct. Some counters for Monte Carlo updates
are initialized here.

### Arguments
* S -> A StochOMSolver struct.

### Returns
* MC -> A StochOMMC struct.

See also: [`StochOMMC`](@ref).
"""
function init_mc(S::StochOMSolver)
    seed = rand(1:100000000)
    rng = MersenneTwister(seed)
    #
    Macc = zeros(I64, 7)
    Mtry = zeros(I64, 7)
    #
    MC = StochOMMC(rng, Macc, Mtry)

    return MC
end

"""
    init_element(MC::StochOMMC, SC::StochOMContext)

Try to initialize a StochOMElement struct. In other words, we should
randomize the configurations for future Monte Carlo sampling here.

### Arguments
* MC -> A StochOMMC struct.
* SC -> A StochOMContext struct.

### Returns
* SE -> A StochOMElement struct.

See also: [`StochOMElement`](@ref).
"""
function init_element(MC::StochOMMC, SC::StochOMContext)
    wmin = get_b("wmin")
    wmax = get_b("wmax")
    nbox = get_s("nbox")
    sbox = get_s("sbox")
    wbox = get_s("wbox")

    # Generate weights randomly
    _Know = rand(MC.rng, 2:nbox)
    _weight = zeros(F64, _Know)
    for i = 1:_Know
        _weight[i] = rand(MC.rng, F64)
    end
    _weight[end] = 1.0

    # Sort weights, make sure the sum of weights is always 1.0.
    sort!(_weight)
    weight = diff(_weight)
    insert!(weight, 1, _weight[1])
    sort!(weight)

    # Make sure that each weight is larger than sbox.
    plus_count = 1
    minus_count = _Know
    while weight[plus_count] < sbox
        while weight[minus_count] < 2 * sbox
            minus_count = minus_count - 1
        end
        weight[plus_count] = weight[plus_count] + sbox
        weight[minus_count] = weight[minus_count] - sbox
        plus_count = plus_count + 1
    end

    # Create some boxes with random c, w, and h.
    C = Box[]
    Œõ = zeros(F64, length(SC.G·µ•), nbox)
    Œî = 0.0
    #
    for k = 1:_Know
        c = wmin + wbox / 2.0 + (wmax - wmin - wbox) * rand(MC.rng, F64)
        w = wbox + (min(2.0 * (c - wmin), 2.0 * (wmax - c)) - wbox) * rand(MC.rng, F64)
        while !constraints(c - w/2, c + w/2)
            c = wmin + wbox / 2.0 + (wmax - wmin - wbox) * rand(MC.rng, F64)
            w = wbox + (min(2.0 * (c - wmin), 2.0 * (wmax - c)) - wbox) * rand(MC.rng, F64)
        end
        h = weight[k] / w
        R = Box(h, w, c)
        push!(C, R)
        Œõ[:,k] .= eval_lambda(R, SC.grid, SC.ùïä·µ•)
    end
    #
    # Calculate Green's function and relative error using boxes
    G = calc_green(Œõ, _Know)
    Œî = calc_error(G, SC.G·µ•, SC.œÉ¬π)

    return StochOMElement(C, Œõ, G, Œî)
end

"""
    init_context(S::StochOMSolver, grid::AbstractGrid)

Try to initialize the key members of a StochOMContext struct.

### Arguments
* S -> A StochOMSolver struct.
* grid -> Grid for input data.

### Returns
* C·µ• -> Field configurations for all attempts.
* Œî·µ• -> Errors for all attempts.
* ùïä·µ• -> Interpolators for the Œõ functions.

See also: [`StochOMContext`](@ref).
"""
function init_context(S::StochOMSolver, grid::AbstractGrid)
    wmin = get_b("wmin")
    wmax = get_b("wmax")
    ntry = get_s("ntry")
    nbox = get_s("nbox")

    # If we increase nmesh gradually, perhaps we could get more precise
    # interpolants ùïä·µ•.
    nmesh = get_b("nmesh") # nmesh = 101, 201, 301, ...
    ngrid = get_b("ngrid")
    @assert ngrid == length(grid)

    # Initialize errors
    Œî·µ• = zeros(F64, ntry)

    # Initialize field configurations (boxes)
    C·µ• = []
    for _ = 1:ntry
        C = Box[]
        for _ = 1:nbox
            push!(C, Box(0.0, 0.0, 0.0))
        end
        push!(C·µ•, C)
    end

    # Initialize interpolants ùïä·µ•
    # It is useful only when the input data is in imaginary time axis.
    ùïä·µ• = Vector{CubicSplineInterpolation}(undef, ngrid)
    #
    if get_b("grid") in ("ftime", "fpart", "btime", "bpart")
        # Create linear mesh for the interpolants
        am = LinearMesh(nmesh, wmin, wmax)

        # Calculate the interpolants at the nodes
        #
        # Initialize memory
        # See below remarks for the Œõ function
        Œõ = zeros(F64, ngrid, nmesh)
        #
        # Just evaluate the interpolants
        for m in eachindex(am)
            if m > 1
                # Create linear mesh for the integrand
                cm = LinearMesh(nmesh, wmin, am[m])
                #
                # Calculate the integrand, i.e., the kernel.
                K = make_kernel(cm, grid)
                #
                # Calculate the integral using trapz rule. Perhaps more
                # precise algorithms should be used.
                for i = 1:ngrid
                    Œõ[i,m] = trapz(cm, K[i,:])
                end
            end
        end

        # Create CubicSplineInterpolation structs in time grid œÑ
        for i = 1:ngrid
            ùïä·µ•[i] = CubicSplineInterpolation(Œõ[i,:], am.mesh)
        end
    end

    return C·µ•, Œî·µ•, ùïä·µ•
end

#=
*Remarks* :

Here, we would like to calculate the ``\Lambda`` function:

```math
\begin{equation}
\Lambda(\omega_n) = \int^{\infty}_{-\infty}
    d\omega~K(\omega_n,\omega) A(\omega).
\end{equation}
```

For a given box ``R``, its contribution to the ``\Lambda`` function is

```math
\begin{equation}
Œõ_{R}(\omega_n) = h \int^{c+w/2}_{c-w/2}
    d\omega~K(\omega_n,\omega),
\end{equation}
```

where ``h``, ``w``, and ``c`` denote the height, width, and center of the
box ``R``. Next, we will show you how to calculate ``\Lambda_R(\omega_n)``.
We can use `sympy` to derive the integral formula.

---

**A** For fermionic Matsubara frequency system.

```math
\begin{equation}
Œõ_{R}(\omega_n) = h \int^{c+w/2}_{c-w/2}
    d\omega~\frac{1}{i\omega_n - \omega}.
\end{equation}
```

The python codes are as follows:

```python
>>> from sympy import *
>>> w, x, a, b, h = symbols("w x a b h")
>>> integrate(h/(w - x), (x, a, b))
h*log(a - w) - h*log(b - w)
```

So, we have

```math
\begin{equation}
Œõ_{R}(\omega_n) = h \log
    \left(\frac{i\omega_n - c + w/2}{i\omega_n - c - w/2}\right).
\end{equation}
```

---

**B** For bosonic Matsubara frequency system.

```math
\begin{equation}
Œõ_{R}(\omega_n) = h \int^{c+w/2}_{c-w/2}
    d\omega~\frac{\omega}{i\omega_n - \omega}.
\end{equation}
```

The python codes are as follows:

```python
>>> from sympy import *
>>> w, x, a, b, h = symbols("w x a b h")
>>> f = integrate(h*x/(w - x), (x, a, b))
>>> simplify(f)
h*(a - b + w*log(a - w) - w*log(b - w))
```

So, we have

```math
\begin{equation}
Œõ_{R}(\omega_n) = h * \left[-w + i\omega_n * \log
    \left(\frac{i\omega_n - c + w/2}{i\omega_n - c - w/2}\right)
\right].
\end{equation}
```

---

**C** For bosonic Matsubara frequency system (symmetric version).

```math
\begin{equation}
Œõ_{R}(\omega_n) = h \int^{c+w/2}_{c-w/2}
    d\omega~\frac{-2\omega^2}{\omega_n^2 + \omega^2}.
\end{equation}
```

We can utilize the following formula:

```math
\begin{equation}
\int \frac{x^2}{a^2 + x^2} dx =
    x - a~\tan^{-1}\left(\frac{x}{a}\right)
\end{equation}
```

Thus,

```math
\begin{equation}
Œõ_{R}(\omega_n) = -2h
    \left[
        \omega - \omega_n \tan^{-1} \left(\frac{\omega}{\omega_n}\right)
    \right]
    \bigg|^{c+w/2}_{c-w/2}
\end{equation}
```

Finally, we have

```math
\begin{equation}
Œõ_{R}(\omega_n) = -2h
    \left[
        w +
        \omega_n \tan^{-1} \left(\frac{c - w/2}{\omega_n}\right)
        -
        \omega_n \tan^{-1} \left(\frac{c + w/2}{\omega_n}\right)
    \right]
\end{equation}
```

---

**D** For imaginary time axis.

If the kernel function ``K`` is defined in imaginary time axis, then the
``\Lambda`` and ``\Lambda_R`` functions become:

```math
\begin{equation}
\Lambda(\tau) = \int^{\infty}_{-\infty}
    d\omega~K(\tau,\omega) A(\omega),
\end{equation}
```

```math
\begin{equation}
Œõ_{R}(\tau) = h \int^{c+w/2}_{c-w/2}
    d\omega~K(\tau,\omega).
\end{equation}
```

We just evaluate the following integral numerically

```math
\begin{equation}
ùïä_R(\tau,Œ©) = \int^{Œ©}_{\omega_{\text{min}}}
    d\omega~K(\tau,\omega).
\end{equation}
```

We notice that `` \Omega \in (\omega_{\text{min}},\omega_{\text{max}}]``,
and ``\Omega`` is defined in a dense linear mesh. Then we reach

```math
\begin{equation}
Œõ_{R}(\tau) = h \left[ ùïä_R(\tau,c + w/2) - ùïä_R(\tau,c - w/2) \right].
\end{equation}
```

In the present implementation, Eq.(13) is evaluated by trapz algorithm,
and Eq.(14) is evaluate using cubic spline interpolation.

We have implemented the above formulas in `eval_lambda()`.
=#

"""
    eval_lambda(
        r::Box,
        grid::FermionicMatsubaraGrid,
        ùïä::Vector{<:AbstractInterpolation}
    )

Try to calculate the contribution of a given box `r` to the Œõ function.
This function works for FermionicMatsubaraGrid only. Because there is an
analytic expression for this case, ùïä is useless.

Actually, ùïä is undefined here. See init_context().

### Arguments
* r    -> A box or rectangle.
* grid -> Imaginary axis grid for input data.
* ùïä    -> An interpolant.

### Returns
* Œõ -> Œõ(iœâ‚Çô) function, 1D function.

See also: [`FermionicMatsubaraGrid`](@ref).
"""
function eval_lambda(
    r::Box,
    grid::FermionicMatsubaraGrid,
    ùïä::Vector{<:AbstractInterpolation}
    )
    # Get left and right boundaries of the given box
    e‚ÇÅ = r.c - 0.5 * r.w
    e‚ÇÇ = r.c + 0.5 * r.w

    # Evaluate Œõ
    iw = im * grid.œâ
    Œõ = @. r.h * log((iw - e‚ÇÅ) / (iw - e‚ÇÇ))

    return vcat(real(Œõ), imag(Œõ))
end

"""
    eval_lambda(
        r::Box,
        grid::FermionicFragmentMatsubaraGrid,
        ùïä::Vector{<:AbstractInterpolation}
    )

Try to calculate the contribution of a given box `r` to the Œõ function.
This function works for FermionicFragmentMatsubaraGrid only. Because there
is an analytic expression for this case, ùïä is useless.

Actually, ùïä is undefined here. See init_context().

### Arguments
* r    -> A box or rectangle.
* grid -> Imaginary axis grid for input data.
* ùïä    -> An interpolant.

### Returns
* Œõ -> Œõ(iœâ‚Çô) function, 1D function.

See also: [`FermionicFragmentMatsubaraGrid`](@ref).
"""
function eval_lambda(
    r::Box,
    grid::FermionicFragmentMatsubaraGrid,
    ùïä::Vector{<:AbstractInterpolation}
    )
    # Get left and right boundaries of the given box
    e‚ÇÅ = r.c - 0.5 * r.w
    e‚ÇÇ = r.c + 0.5 * r.w

    # Evaluate Œõ
    iw = im * grid.œâ
    Œõ = @. r.h * log((iw - e‚ÇÅ) / (iw - e‚ÇÇ))

    return vcat(real(Œõ), imag(Œõ))
end

"""
    eval_lambda(
        r::Box,
        grid::FermionicImaginaryTimeGrid,
        ùïä::Vector{<:AbstractInterpolation}
    )

Try to calculate the contribution of a given box `r` to the Œõ function.
This function works for FermionicImaginaryTimeGrid only. Since there is
no analytic expressions for this case, the cubic spline interpolation
algorithm is adopted. Here, ùïä is initialized in init_context().

### Arguments
* r    -> A box or rectangle.
* grid -> Imaginary axis grid for input data.
* ùïä    -> An interpolant.

### Returns
* Œõ -> Œõ(œÑ) function, 1D function.

See also: [`FermionicImaginaryTimeGrid`](@ref).
"""
function eval_lambda(
    r::Box,
    grid::FermionicImaginaryTimeGrid,
    ùïä::Vector{<:AbstractInterpolation}
    )
    # Get left and right boundaries of the given box
    e‚ÇÅ = r.c - 0.5 * r.w
    e‚ÇÇ = r.c + 0.5 * r.w

    # Initialize Œõ function
    ntime = grid.ntime
    Œõ = zeros(F64, ntime)

    # ùïä·µ¢(e‚ÇÇ): integral boundary is from wmin to e‚ÇÇ
    # ùïä·µ¢(e‚ÇÅ): integral boundary is from wmin to e‚ÇÅ
    for i = 1:ntime
        Œõ[i] = ( ùïä[i](e‚ÇÇ) - ùïä[i](e‚ÇÅ) ) * r.h
    end

    return Œõ
end

"""
    eval_lambda(
        r::Box,
        grid::FermionicFragmentTimeGrid,
        ùïä::Vector{<:AbstractInterpolation}
    )

Try to calculate the contribution of a given box `r` to the Œõ function.
This function works for FermionicFragmentTimeGrid only. Since there is
no analytic expressions for this case, the cubic spline interpolation
algorithm is adopted. Here, ùïä is initialized in init_context().

### Arguments
* r    -> A box or rectangle.
* grid -> Imaginary axis grid for input data.
* ùïä    -> An interpolant.

### Returns
* Œõ -> Œõ(œÑ) function, 1D function.

See also: [`FermionicFragmentTimeGrid`](@ref).
"""
function eval_lambda(
    r::Box,
    grid::FermionicFragmentTimeGrid,
    ùïä::Vector{<:AbstractInterpolation}
    )
    # Get left and right boundaries of the given box
    e‚ÇÅ = r.c - 0.5 * r.w
    e‚ÇÇ = r.c + 0.5 * r.w

    # Initialize Œõ function
    ntime = grid.ntime
    Œõ = zeros(F64, ntime)

    # ùïä·µ¢(e‚ÇÇ): integral boundary is from wmin to e‚ÇÇ
    # ùïä·µ¢(e‚ÇÅ): integral boundary is from wmin to e‚ÇÅ
    for i = 1:ntime
        Œõ[i] = ( ùïä[i](e‚ÇÇ) - ùïä[i](e‚ÇÅ) ) * r.h
    end

    return Œõ
end

"""
    eval_lambda(
        r::Box,
        grid::BosonicMatsubaraGrid,
        ùïä::Vector{<:AbstractInterpolation}
    )

Try to calculate the contribution of a given box `r` to the Œõ function.
This function works for BosonicMatsubaraGrid only. Because there is an
analytic expression for this case, ùïä is useless.

Actually, ùïä is undefined here. See init_context().

### Arguments
* r    -> A box or rectangle.
* grid -> Imaginary axis grid for input data.
* ùïä    -> An interpolant.

### Returns
* Œõ -> Œõ(iœâ‚Çô) function, 1D function.

See also: [`BosonicMatsubaraGrid`](@ref).
"""
function eval_lambda(
    r::Box,
    grid::BosonicMatsubaraGrid,
    ùïä::Vector{<:AbstractInterpolation}
    )
    # Get type of bosonic kernel
    ktype = get_b("ktype")

    # Get left and right boundaries of the given box
    e‚ÇÅ = r.c - 0.5 * r.w
    e‚ÇÇ = r.c + 0.5 * r.w

    # Evaluate Œõ
    if ktype == "bsymm"
        Œõ = @. atan( e‚ÇÅ / grid.œâ ) - atan( e‚ÇÇ / grid.œâ )
        Œõ = -2.0 * r.h * (r.w .+ grid.œâ .* Œõ)
        return Œõ
    else
        iw = im * grid.œâ
        Œõ = @. r.h * (-r.w + iw * log((iw - e‚ÇÅ) / (iw - e‚ÇÇ)))
        return vcat(real(Œõ), imag(Œõ))
    end
end

"""
    eval_lambda(
        r::Box,
        grid::BosonicFragmentMatsubaraGrid,
        ùïä::Vector{<:AbstractInterpolation}
    )

Try to calculate the contribution of a given box `r` to the Œõ function.
This function works for BosonicFragmentMatsubaraGrid only. Because there
is an analytic expression for this case, ùïä is useless.

Actually, ùïä is undefined here. See init_context().

### Arguments
* r    -> A box or rectangle.
* grid -> Imaginary axis grid for input data.
* ùïä    -> An interpolant.

### Returns
* Œõ -> Œõ(iœâ‚Çô) function, 1D function.

See also: [`BosonicFragmentMatsubaraGrid`](@ref).
"""
function eval_lambda(
    r::Box,
    grid::BosonicFragmentMatsubaraGrid,
    ùïä::Vector{<:AbstractInterpolation}
    )
    # Get type of bosonic kernel
    ktype = get_b("ktype")

    # Get left and right boundaries of the given box
    e‚ÇÅ = r.c - 0.5 * r.w
    e‚ÇÇ = r.c + 0.5 * r.w

    # Evaluate Œõ
    if ktype == "bsymm"
        Œõ = @. atan( e‚ÇÅ / grid.œâ ) - atan( e‚ÇÇ / grid.œâ )
        Œõ = -2.0 * r.h * (r.w .+ grid.œâ .* Œõ)
        return Œõ
    else
        iw = im * grid.œâ
        Œõ = @. r.h * (-r.w + iw * log((iw - e‚ÇÅ) / (iw - e‚ÇÇ)))
        return vcat(real(Œõ), imag(Œõ))
    end
end

"""
    eval_lambda(
        r::Box,
        grid::BosonicImaginaryTimeGrid,
        ùïä::Vector{<:AbstractInterpolation}
    )

Try to calculate the contribution of a given box `r` to the Œõ function.
This function works for BosonicImaginaryTimeGrid only. Since there is
no analytic expressions for this case, the cubic spline interpolation
algorithm is adopted. Here, ùïä is initialized in init_context().

### Arguments
* r    -> A box or rectangle.
* grid -> Imaginary axis grid for input data.
* ùïä    -> An interpolant.

### Returns
* Œõ -> Œõ(œÑ) function, 1D function.

See also: [`BosonicImaginaryTimeGrid`](@ref).
"""
function eval_lambda(
    r::Box,
    grid::BosonicImaginaryTimeGrid,
    ùïä::Vector{<:AbstractInterpolation}
    )
    # Get left and right boundaries of the given box
    e‚ÇÅ = r.c - 0.5 * r.w
    e‚ÇÇ = r.c + 0.5 * r.w

    # Initialize Œõ function
    ntime = grid.ntime
    Œõ = zeros(F64, ntime)

    # ùïä·µ¢(e‚ÇÇ): integral boundary is from wmin to e‚ÇÇ
    # ùïä·µ¢(e‚ÇÅ): integral boundary is from wmin to e‚ÇÅ
    for i = 1:ntime
        Œõ[i] = ( ùïä[i](e‚ÇÇ) - ùïä[i](e‚ÇÅ) ) * r.h
    end

    return Œõ
end

"""
    eval_lambda(
        r::Box,
        grid::BosonicFragmentTimeGrid,
        ùïä::Vector{<:AbstractInterpolation}
    )

Try to calculate the contribution of a given box `r` to the Œõ function.
This function works for BosonicFragmentTimeGrid only. Since there is
no analytic expressions for this case, the cubic spline interpolation
algorithm is adopted. Here, ùïä is initialized in init_context().

### Arguments
* r    -> A box or rectangle.
* grid -> Imaginary axis grid for input data.
* ùïä    -> An interpolant.

### Returns
* Œõ -> Œõ(œÑ) function, 1D function.

See also: [`BosonicFragmentTimeGrid`](@ref).
"""
function eval_lambda(
    r::Box,
    grid::BosonicFragmentTimeGrid,
    ùïä::Vector{<:AbstractInterpolation}
    )
    # Get left and right boundaries of the given box
    e‚ÇÅ = r.c - 0.5 * r.w
    e‚ÇÇ = r.c + 0.5 * r.w

    # Initialize Œõ function
    ntime = grid.ntime
    Œõ = zeros(F64, ntime)

    # ùïä·µ¢(e‚ÇÇ): integral boundary is from wmin to e‚ÇÇ
    # ùïä·µ¢(e‚ÇÅ): integral boundary is from wmin to e‚ÇÅ
    for i = 1:ntime
        Œõ[i] = ( ùïä[i](e‚ÇÇ) - ùïä[i](e‚ÇÅ) ) * r.h
    end

    return Œõ
end

"""
    calc_error(G::Vector{F64}, G·µ•::Vector{F64}, œÉ¬π::Vector{F64})

Try to calculate œá¬≤. Here `G·µ•` and `œÉ¬π` denote the raw correlator and
related standard deviation. `G` means the reproduced correlator.

### Arguments
See above explanations.

### Returns
* Œî -> œá¬≤, distance between reconstructed and raw correlators.

See also: [`calc_green`](@ref).
"""
function calc_error(G::Vector{F64}, G·µ•::Vector{F64}, œÉ¬π::Vector{F64})
    return sum( ( (G .- G·µ•) .* œÉ¬π ) .^ 2.0 )
end

"""
    calc_green(Œõ::Array{F64,2}, nk::I64)

Try to reconstruct the correlator via the field configuration. Now this
function is called by init_element(). But perhaps we can use it in last().

### Arguments
* Œõ -> The Œõ function. See above remarks.
* nk -> Current number of boxes.

### Returns
* G -> Reconstructed Green's function.

See also: [`calc_error`](@ref).
"""
function calc_green(Œõ::Array{F64,2}, nk::I64)
    ngrid, nbox = size(Œõ)
    @assert nk ‚â§ nbox

    G = zeros(F64, ngrid)
    for k = 1:nk
        for g = 1:ngrid
            G[g] = G[g] + Œõ[g,k]
        end
    end

    return G
end

"""
    calc_norm(C::Vector{Box})

Calculate the total area of all boxes. Now this function is not used.

### Arguments
* C -> The current Monte Carlo field configuration.

### Returns
* norm -> Area of all boxes.
"""
function calc_norm(C::Vector{Box})
    norm = sum(map(x -> x.h * x.w, C))
    return norm
end

"""
    constraints(e‚ÇÅ::F64, e‚ÇÇ::F64)

This function is used to judge whether a given box overlapes with the
forbidden zone. Here `e‚ÇÅ` and `e‚ÇÇ` denote the left and right boundaries
of the box.

### Arguments
See above explanations.

### Returns
* ex -> Boolean, whether a given box is valid.

"""
function constraints(e‚ÇÅ::F64, e‚ÇÇ::F64)
    exclude = get_b("exclude")
    @assert e‚ÇÅ ‚â§ e‚ÇÇ

    if !isa(exclude, Missing)
        for i in eachindex(exclude)
            if e‚ÇÅ ‚â§ exclude[i][1] ‚â§ e‚ÇÇ ‚â§ exclude[i][2]
                return false
            end
            if exclude[i][1] ‚â§ e‚ÇÅ ‚â§ exclude[i][2] ‚â§ e‚ÇÇ
                return false
            end
            if exclude[i][1] ‚â§ e‚ÇÅ ‚â§ e‚ÇÇ ‚â§ exclude[i][2]
                return false
            end
            if e‚ÇÅ ‚â§ exclude[i][1] ‚â§ exclude[i][2] ‚â§ e‚ÇÇ
                return false
            end
        end
    end

    return true
end

"""
    try_insert(
        MC::StochOMMC,
        SE::StochOMElement,
        SC::StochOMContext,
        dacc::F64
    )

Insert a new box into the field configuration.

### Arguments
* MC -> A StochOMMC struct. It containts some counters.
* SE -> A StochOMElement struct. It contains Monte Carlo configurations.
* SC -> A StochOMContext struct. It contains grid, mesh, and G·µ•.
* dacc -> A predefined parameter used to calculate transition probability.

### Returns
N/A
"""
function try_insert(
    MC::StochOMMC,
    SE::StochOMElement,
    SC::StochOMContext,
    dacc::F64
    )
    sbox = get_s("sbox")
    wbox = get_s("wbox")
    wmin = get_b("wmin")
    wmax = get_b("wmax")
    csize = length(SE.C)

    # Choose a box randomly
    t = rand(MC.rng, 1:csize)

    # Check area of this box
    R = SE.C[t]
    if R.h * R.w ‚â§ 2.0 * sbox
        return
    end

    # Determine minimum and maximum areas of the new box
    dx_min = sbox
    dx_max = R.h * R.w - sbox
    if dx_max ‚â§ dx_min
        return
    end

    # Determine parameters for the new box
    r‚ÇÅ = rand(MC.rng, F64)
    r‚ÇÇ = rand(MC.rng, F64)
    #
    c = (wmin + wbox / 2.0) + (wmax - wmin - wbox) * r‚ÇÅ
    #
    w_new_max = 2.0 * min(wmax - c, c - wmin)
    dx = Pdx(dx_min, dx_max, MC.rng)
    #
    h = dx / w_new_max + (dx / wbox - dx / w_new_max) * r‚ÇÇ
    w = dx / h

    # Rnew will be used to update Box t, while Radd is the new box.
    if !constraints(c - w/2, c + w/2)
        return
    end
    Rnew = Box(R.h - dx / R.w, R.w, R.c)
    Radd = Box(h, w, c)

    # Calculate update for Œõ
    G‚ÇÅ = SE.Œõ[:,t]
    G‚ÇÇ = eval_lambda(Rnew, SC.grid, SC.ùïä·µ•)
    G‚ÇÉ = eval_lambda(Radd, SC.grid, SC.ùïä·µ•)

    # Calculate new Œî function, it is actually the error function.
    Œî = calc_error(SE.G - G‚ÇÅ + G‚ÇÇ + G‚ÇÉ, SC.G·µ•, SC.œÉ¬π)

    # Apply the Metropolis algorithm
    if rand(MC.rng, F64) < ((SE.Œî/Œî) ^ (1.0 + dacc))
        # Update box t
        SE.C[t] = Rnew

        # Add new Box
        push!(SE.C, Radd)

        # Update Œî, G, and Œõ.
        SE.Œî = Œî
        @. SE.G = SE.G - G‚ÇÅ + G‚ÇÇ + G‚ÇÉ
        @. SE.Œõ[:,t] = G‚ÇÇ
        @. SE.Œõ[:,csize+1] = G‚ÇÉ

        # Update the counter
        MC.Macc[1] = MC.Macc[1] + 1
    end

    # Update the counter
    MC.Mtry[1] = MC.Mtry[1] + 1
end

"""
    try_remove(
        MC::StochOMMC,
        SE::StochOMElement,
        SC::StochOMContext,
        dacc::F64
    )

Remove an old box from the field configuration.

### Arguments
* MC -> A StochOMMC struct. It containts some counters.
* SE -> A StochOMElement struct. It contains Monte Carlo configurations.
* SC -> A StochOMContext struct. It contains grid, mesh, and G·µ•.
* dacc -> A predefined parameter used to calculate transition probability.

### Returns
N/A
"""
function try_remove(
    MC::StochOMMC,
    SE::StochOMElement,
    SC::StochOMContext,
    dacc::F64
    )
    csize = length(SE.C)

    # Choose two boxes randomly
    # Box t‚ÇÅ will be removed, while box t‚ÇÇ will be modified.
    t‚ÇÅ = rand(MC.rng, 1:csize)
    t‚ÇÇ = rand(MC.rng, 1:csize)
    #
    while t‚ÇÅ == t‚ÇÇ
        t‚ÇÇ = rand(MC.rng, 1:csize)
    end
    #
    if t‚ÇÅ < t‚ÇÇ
        t‚ÇÅ, t‚ÇÇ = t‚ÇÇ, t‚ÇÅ
    end

    # Get box t‚ÇÅ and box t‚ÇÇ
    R‚ÇÅ = SE.C[t‚ÇÅ]
    R‚ÇÇ = SE.C[t‚ÇÇ]
    R‚Çë = SE.C[end]

    # Generate new box t‚ÇÇ
    dx = R‚ÇÅ.h * R‚ÇÅ.w
    R‚ÇÇ‚Çô = Box(R‚ÇÇ.h + dx / R‚ÇÇ.w, R‚ÇÇ.w, R‚ÇÇ.c)

    # Calculate update for Œõ
    G‚ÇÅ = SE.Œõ[:,t‚ÇÅ]
    G‚ÇÇ = SE.Œõ[:,t‚ÇÇ]
    G‚Çë = SE.Œõ[:,csize]
    G‚ÇÇ‚Çô = eval_lambda(R‚ÇÇ‚Çô, SC.grid, SC.ùïä·µ•)

    # Calculate new Œî function, it is actually the error function.
    Œî = calc_error(SE.G - G‚ÇÅ - G‚ÇÇ + G‚ÇÇ‚Çô, SC.G·µ•, SC.œÉ¬π)

    # Apply the Metropolis algorithm
    if rand(MC.rng, F64) < ((SE.Œî/Œî) ^ (1.0 + dacc))
        # Update box t‚ÇÇ
        SE.C[t‚ÇÇ] = R‚ÇÇ‚Çô

        # Backup the last box in box t‚ÇÅ
        if t‚ÇÅ < csize
            SE.C[t‚ÇÅ] = R‚Çë
        end

        # Delete the last box, since its value has been stored in t‚ÇÅ.
        pop!(SE.C)

        # Update Œî, G, and Œõ.
        SE.Œî = Œî
        @. SE.G = SE.G - G‚ÇÅ - G‚ÇÇ + G‚ÇÇ‚Çô
        @. SE.Œõ[:,t‚ÇÇ] = G‚ÇÇ‚Çô
        if t‚ÇÅ < csize
            @. SE.Œõ[:,t‚ÇÅ] = G‚Çë
        end

        # Update the counter
        MC.Macc[2] = MC.Macc[2] + 1
    end

    # Update the counter
    MC.Mtry[2] = MC.Mtry[2] + 1
end

"""
    try_shift(
        MC::StochOMMC,
        SE::StochOMElement,
        SC::StochOMContext,
        dacc::F64
    )

Change the position of given box in the field configuration.

### Arguments
* MC -> A StochOMMC struct. It containts some counters.
* SE -> A StochOMElement struct. It contains Monte Carlo configurations.
* SC -> A StochOMContext struct. It contains grid, mesh, and G·µ•.
* dacc -> A predefined parameter used to calculate transition probability.

### Returns
N/A
"""
function try_shift(
    MC::StochOMMC,
    SE::StochOMElement,
    SC::StochOMContext,
    dacc::F64
    )
    wmin = get_b("wmin")
    wmax = get_b("wmax")
    csize = length(SE.C)

    # Choose a box randomly
    t = rand(MC.rng, 1:csize)

    # Retreive the box t
    R = SE.C[t]

    # Determine left and right boundaries for the center of the box
    dx_min = wmin + R.w / 2.0 - R.c
    dx_max = wmax - R.w / 2.0 - R.c
    if dx_max ‚â§ dx_min
        return
    end

    # Calculate Œ¥c and generate shifted box
    Œ¥c = Pdx(dx_min, dx_max, MC.rng)
    if !constraints(R.c + Œ¥c - R.w/2, R.c + Œ¥c + R.w/2)
        return
    end
    R‚Çô = Box(R.h, R.w, R.c + Œ¥c)

    # Calculate update for Œõ
    G‚ÇÅ = SE.Œõ[:,t]
    G‚ÇÇ = eval_lambda(R‚Çô, SC.grid, SC.ùïä·µ•)

    # Calculate new Œî function, it is actually the error function.
    Œî = calc_error(SE.G - G‚ÇÅ + G‚ÇÇ, SC.G·µ•, SC.œÉ¬π)

    # Apply the Metropolis algorithm
    if rand(MC.rng, F64) < ((SE.Œî/Œî) ^ (1.0 + dacc))
        # Update box t
        SE.C[t] = R‚Çô

        # Update Œî, G, and Œõ.
        SE.Œî = Œî
        @. SE.G = SE.G - G‚ÇÅ + G‚ÇÇ
        @. SE.Œõ[:,t] = G‚ÇÇ

        # Update the counter
        MC.Macc[3] = MC.Macc[3] + 1
    end

    # Update the counter
    MC.Mtry[3] = MC.Mtry[3] + 1
end

"""
    try_width(
        MC::StochOMMC,
        SE::StochOMElement,
        SC::StochOMContext,
        dacc::F64
    )

Change the width and height of given box in the field configuration. Note
that the box's area is kept.

### Arguments
* MC -> A StochOMMC struct. It containts some counters.
* SE -> A StochOMElement struct. It contains Monte Carlo configurations.
* SC -> A StochOMContext struct. It contains grid, mesh, and G·µ•.
* dacc -> A predefined parameter used to calculate transition probability.

### Returns
N/A
"""
function try_width(
    MC::StochOMMC,
    SE::StochOMElement,
    SC::StochOMContext,
    dacc::F64
    )
    wbox = get_s("wbox")
    wmin = get_b("wmin")
    wmax = get_b("wmax")
    csize = length(SE.C)

    # Choose a box randomly
    t = rand(MC.rng, 1:csize)

    # Retreive the box t
    R = SE.C[t]

    # Determine left and right boundaries for the width of the box
    weight = R.h * R.w
    dx_min = wbox - R.w
    dx_max = min(2.0 * (R.c - wmin), 2.0 * (wmax - R.c)) - R.w
    if dx_max ‚â§ dx_min
        return
    end

    # Calculate Œ¥w and generate new box
    dw = Pdx(dx_min, dx_max, MC.rng)
    w = R.w + dw
    h = weight / w
    c = R.c
    if !constraints(c - w/2, c + w/2)
        return
    end
    R‚Çô = Box(h, w, c)

    # Calculate update for Œõ
    G‚ÇÅ = SE.Œõ[:,t]
    G‚ÇÇ = eval_lambda(R‚Çô, SC.grid, SC.ùïä·µ•)

    # Calculate new Œî function, it is actually the error function.
    Œî = calc_error(SE.G - G‚ÇÅ + G‚ÇÇ, SC.G·µ•, SC.œÉ¬π)

    # Apply the Metropolis algorithm
    if rand(MC.rng, F64) < ((SE.Œî/Œî) ^ (1.0 + dacc))
        # Update box t
        SE.C[t] = R‚Çô

        # Update Œî, G, and Œõ.
        SE.Œî = Œî
        @. SE.G = SE.G - G‚ÇÅ + G‚ÇÇ
        @. SE.Œõ[:,t] = G‚ÇÇ

        # Update the counter
        MC.Macc[4] = MC.Macc[4] + 1
    end

    # Update the counter
    MC.Mtry[4] = MC.Mtry[4] + 1
end

"""
    try_height(
        MC::StochOMMC,
        SE::StochOMElement,
        SC::StochOMContext,
        dacc::F64
    )

Change the heights of two given boxes in the field configuration.

### Arguments
* MC -> A StochOMMC struct. It containts some counters.
* SE -> A StochOMElement struct. It contains Monte Carlo configurations.
* SC -> A StochOMContext struct. It contains grid, mesh, and G·µ•.
* dacc -> A predefined parameter used to calculate transition probability.

### Returns
N/A
"""
function try_height(
    MC::StochOMMC,
    SE::StochOMElement,
    SC::StochOMContext,
    dacc::F64
    )
    sbox  = get_s("sbox")
    csize = length(SE.C)

    # Choose two boxes randomly
    t‚ÇÅ = rand(MC.rng, 1:csize)
    t‚ÇÇ = rand(MC.rng, 1:csize)
    #
    while t‚ÇÅ == t‚ÇÇ
        t‚ÇÇ = rand(MC.rng, 1:csize)
    end

    # Get box t‚ÇÅ and box t‚ÇÇ
    R‚ÇÅ = SE.C[t‚ÇÅ]
    R‚ÇÇ = SE.C[t‚ÇÇ]

    # Determine left and right boundaries for the height of the box t‚ÇÅ
    w‚ÇÅ = R‚ÇÅ.w
    w‚ÇÇ = R‚ÇÇ.w
    h‚ÇÅ = R‚ÇÅ.h
    h‚ÇÇ = R‚ÇÇ.h
    dx_min = sbox / w‚ÇÅ - h‚ÇÅ
    dx_max = (h‚ÇÇ - sbox / w‚ÇÇ) * w‚ÇÇ / w‚ÇÅ
    if dx_max ‚â§ dx_min
        return
    end

    # Calculate Œ¥h and generate new box t‚ÇÅ and box t‚ÇÇ
    dh = Pdx(dx_min, dx_max, MC.rng)
    R‚ÇÅ‚Çô = Box(R‚ÇÅ.h + dh, R‚ÇÅ.w, R‚ÇÅ.c)
    R‚ÇÇ‚Çô = Box(R‚ÇÇ.h - dh * w‚ÇÅ / w‚ÇÇ, R‚ÇÇ.w, R‚ÇÇ.c)

    # Calculate update for Œõ
    G‚ÇÅA = SE.Œõ[:,t‚ÇÅ]
    G‚ÇÅB = eval_lambda(R‚ÇÅ‚Çô, SC.grid, SC.ùïä·µ•)
    G‚ÇÇA = SE.Œõ[:,t‚ÇÇ]
    G‚ÇÇB = eval_lambda(R‚ÇÇ‚Çô, SC.grid, SC.ùïä·µ•)

    # Calculate new Œî function, it is actually the error function.
    Œî = calc_error(SE.G - G‚ÇÅA + G‚ÇÅB - G‚ÇÇA + G‚ÇÇB, SC.G·µ•, SC.œÉ¬π)

    # Apply the Metropolis algorithm
    if rand(MC.rng, F64) < ((SE.Œî/Œî) ^ (1.0 + dacc))
        # Update box t‚ÇÅ and box t‚ÇÇ
        SE.C[t‚ÇÅ] = R‚ÇÅ‚Çô
        SE.C[t‚ÇÇ] = R‚ÇÇ‚Çô

        # Update Œî, G, and Œõ.
        SE.Œî = Œî
        @. SE.G = SE.G - G‚ÇÅA + G‚ÇÅB - G‚ÇÇA + G‚ÇÇB
        @. SE.Œõ[:,t‚ÇÅ] = G‚ÇÅB
        @. SE.Œõ[:,t‚ÇÇ] = G‚ÇÇB

        # Update the counter
        MC.Macc[5] = MC.Macc[5] + 1
    end

    # Update the counter
    MC.Mtry[5] = MC.Mtry[5] + 1
end

"""
    try_split(
        MC::StochOMMC,
        SE::StochOMElement,
        SC::StochOMContext,
        dacc::F64
    )

Split a given box into two boxes in the field configuration.

### Arguments
* MC -> A StochOMMC struct. It containts some counters.
* SE -> A StochOMElement struct. It contains Monte Carlo configurations.
* SC -> A StochOMContext struct. It contains grid, mesh, and G·µ•.
* dacc -> A predefined parameter used to calculate transition probability.

### Returns
N/A
"""
function try_split(
    MC::StochOMMC,
    SE::StochOMElement,
    SC::StochOMContext,
    dacc::F64
    )
    wbox = get_s("wbox")
    sbox = get_s("sbox")
    wmin = get_b("wmin")
    wmax = get_b("wmax")
    csize = length(SE.C)

    # Choose a box randomly
    t = rand(MC.rng, 1:csize)

    # Retreive the box t
    R‚ÇÅ = SE.C[t]
    if R‚ÇÅ.w ‚â§ 2 * wbox || R‚ÇÅ.w * R‚ÇÅ.h ‚â§ 2.0 * sbox
        return
    end

    # Determine height for new boxes (h and h)
    h = R‚ÇÅ.h

    # Determine width for new boxes (w‚ÇÅ and w‚ÇÇ)
    w‚ÇÅ = wbox + (R‚ÇÅ.w - 2.0 * wbox) * rand(MC.rng, F64)
    w‚ÇÇ = R‚ÇÅ.w - w‚ÇÅ
    if w‚ÇÅ > w‚ÇÇ
        w‚ÇÅ, w‚ÇÇ = w‚ÇÇ, w‚ÇÅ
    end

    # Determine center for new boxes (c‚ÇÅ + Œ¥c‚ÇÅ and c‚ÇÇ + Œ¥c‚ÇÇ)
    c‚ÇÅ = R‚ÇÅ.c - R‚ÇÅ.w / 2.0 + w‚ÇÅ / 2.0
    c‚ÇÇ = R‚ÇÅ.c + R‚ÇÅ.w / 2.0 - w‚ÇÇ / 2.0
    dx_min = wmin + w‚ÇÅ / 2.0 - c‚ÇÅ
    dx_max = wmax - w‚ÇÅ / 2.0 - c‚ÇÅ
    if dx_max ‚â§ dx_min
        return
    end
    Œ¥c‚ÇÅ = Pdx(dx_min, dx_max, MC.rng)
    Œ¥c‚ÇÇ = -1.0 * w‚ÇÅ * Œ¥c‚ÇÅ / w‚ÇÇ
    if !constraints(c‚ÇÅ + Œ¥c‚ÇÅ - w‚ÇÅ/2, c‚ÇÅ + Œ¥c‚ÇÅ + w‚ÇÅ/2) ||
       !constraints(c‚ÇÇ + Œ¥c‚ÇÇ - w‚ÇÇ/2, c‚ÇÇ + Œ¥c‚ÇÇ + w‚ÇÇ/2)
        return
    end

    if (c‚ÇÅ + Œ¥c‚ÇÅ ‚â• wmin + w‚ÇÅ / 2.0) &&
       (c‚ÇÅ + Œ¥c‚ÇÅ ‚â§ wmax - w‚ÇÅ / 2.0) &&
       (c‚ÇÇ + Œ¥c‚ÇÇ ‚â• wmin + w‚ÇÇ / 2.0) &&
       (c‚ÇÇ + Œ¥c‚ÇÇ ‚â§ wmax - w‚ÇÇ / 2.0)

        # Generate two new boxes
        R‚ÇÇ = Box(h, w‚ÇÅ, c‚ÇÅ + Œ¥c‚ÇÅ)
        R‚ÇÉ = Box(h, w‚ÇÇ, c‚ÇÇ + Œ¥c‚ÇÇ)

        # Calculate update for Œõ
        G‚ÇÅ = SE.Œõ[:,t]
        G‚ÇÇ = eval_lambda(R‚ÇÇ, SC.grid, SC.ùïä·µ•)
        G‚ÇÉ = eval_lambda(R‚ÇÉ, SC.grid, SC.ùïä·µ•)

        # Calculate new Œî function, it is actually the error function.
        Œî = calc_error(SE.G - G‚ÇÅ + G‚ÇÇ + G‚ÇÉ, SC.G·µ•, SC.œÉ¬π)

        # Apply the Metropolis algorithm
        if rand(MC.rng, F64) < ((SE.Œî/Œî) ^ (1.0 + dacc))
            # Remove old box t and insert two new boxes
            SE.C[t] = R‚ÇÇ
            push!(SE.C, R‚ÇÉ)

            # Update Œî, G, and Œõ.
            SE.Œî = Œî
            @. SE.G = SE.G - G‚ÇÅ + G‚ÇÇ + G‚ÇÉ
            @. SE.Œõ[:,t] = G‚ÇÇ
            @. SE.Œõ[:,csize+1] = G‚ÇÉ

            # Update the counter
            MC.Macc[6] = MC.Macc[6] + 1
        end
    end

    # Update the counter
    MC.Mtry[6] = MC.Mtry[6] + 1
end

"""
    try_merge(
        MC::StochOMMC,
        SE::StochOMElement,
        SC::StochOMContext,
        dacc::F64
    )

Merge two given boxes into one box in the field configuration.

### Arguments
* MC -> A StochOMMC struct. It containts some counters.
* SE -> A StochOMElement struct. It contains Monte Carlo configurations.
* SC -> A StochOMContext struct. It contains grid, mesh, and G·µ•.
* dacc -> A predefined parameter used to calculate transition probability.

### Returns
N/A
"""
function try_merge(
    MC::StochOMMC,
    SE::StochOMElement,
    SC::StochOMContext,
    dacc::F64
    )
    wmin = get_b("wmin")
    wmax = get_b("wmax")
    csize = length(SE.C)

    # Choose two boxes randomly
    # Box t‚ÇÇ will be removed, while box t‚ÇÅ will be modified.
    t‚ÇÅ = rand(MC.rng, 1:csize)
    t‚ÇÇ = rand(MC.rng, 1:csize)
    #
    while t‚ÇÅ == t‚ÇÇ
        t‚ÇÇ = rand(MC.rng, 1:csize)
    end
    #
    if t‚ÇÅ > t‚ÇÇ
        t‚ÇÅ, t‚ÇÇ = t‚ÇÇ, t‚ÇÅ
    end

    # Get box t‚ÇÅ and box t‚ÇÇ
    R‚ÇÅ = SE.C[t‚ÇÅ]
    R‚ÇÇ = SE.C[t‚ÇÇ]

    # Determine h, w, and c for new box
    weight = R‚ÇÅ.h * R‚ÇÅ.w + R‚ÇÇ.h * R‚ÇÇ.w
    w_new = 0.5 * (R‚ÇÅ.w + R‚ÇÇ.w)
    h_new = weight / w_new
    c_new = R‚ÇÅ.c + (R‚ÇÇ.c - R‚ÇÅ.c) * R‚ÇÇ.h * R‚ÇÇ.w / weight

    # Determine left and right boundaries for the center of the new box
    dx_min = wmin + w_new / 2.0 - c_new
    dx_max = wmax - w_new / 2.0 - c_new
    if dx_max ‚â§ dx_min
        return
    end

    # Calculate Œ¥c and generate new box
    Œ¥c = Pdx(dx_min, dx_max, MC.rng)
    if !constraints(c_new + Œ¥c - w_new/2, c_new + Œ¥c + w_new/2)
        return
    end
    R‚Çô = Box(h_new, w_new, c_new + Œ¥c)

    # Calculate update for Œõ
    G‚ÇÅ = SE.Œõ[:,t‚ÇÅ]
    G‚ÇÇ = SE.Œõ[:,t‚ÇÇ]
    G‚Çë = SE.Œõ[:,csize]
    G‚Çô = eval_lambda(R‚Çô, SC.grid, SC.ùïä·µ•)

    # Calculate new Œî function, it is actually the error function.
    Œî = calc_error(SE.G - G‚ÇÅ - G‚ÇÇ + G‚Çô, SC.G·µ•, SC.œÉ¬π)

    # Apply the Metropolis algorithm
    if rand(MC.rng, F64) < ((SE.Œî/Œî) ^ (1.0 + dacc))
        # Update box t‚ÇÅ with new box
        SE.C[t‚ÇÅ] = R‚Çô

        # Delete box t‚ÇÇ
        if t‚ÇÇ < csize
            SE.C[t‚ÇÇ] = SE.C[end]
        end
        pop!(SE.C)

        # Update Œî, G, and Œõ.
        SE.Œî = Œî
        @. SE.G = SE.G - G‚ÇÅ - G‚ÇÇ + G‚Çô
        @. SE.Œõ[:,t‚ÇÅ] = G‚Çô
        if t‚ÇÇ < csize
            @. SE.Œõ[:,t‚ÇÇ] = G‚Çë
        end

        # Update the counter
        MC.Macc[7] = MC.Macc[7] + 1
    end

    # Update the counter
    MC.Mtry[7] = MC.Mtry[7] + 1
end

#=
*Remarks* : *Probability Density Function*

Every proposed elementary update is parametrized by a real number
``\delta\xi \in [\delta\xi_{min}:\delta\xi_{max}]``. A concrete meaning
of ``\delta\xi`` depends on the elementary update in question. For
instance, ``\delta\xi`` can be a shift of the centre of an existing
rectangle, or the weight of a rectangle to be added. In general, ``\delta
\xi`` are defined so that larger ``|\delta\xi|`` correspond to more
prominent changes in the configuration. `StochOM` randomly generates
values of ``\delta\xi`` according to the following probability density
function:

```math
\begin{equation}
\mathcal{P}(\delta\xi) = N
\exp\left(-\gamma \frac{|\delta\xi|}{X}\right),
\end{equation}
```

```math
\begin{equation}
X \equiv \max(|\delta\xi_{min}|, |\delta\xi_{max}|),
\end{equation}
```

```math
\begin{equation}
N = \frac{\gamma}{X}
\left[
\text{sign}(\delta\xi_{min})\left(e^{-\gamma|\delta\xi_{min}|/X}-1\right) +
\text{sign}(\delta\xi_{max})\left(1-e^{-\gamma|\delta\xi_{max}|/X}\right)
\right]^{-1}
\end{equation}
```

User can change parameter ``\gamma`` to control non-uniformity of the PDF.
=#

"""
    Pdx(xmin::F64, xmax::F64, rng::AbstractRNG)

Try to calculate the value of Œ¥Œæ for every elementary update according to
the probability density function. The actual meaning of Œ¥Œæ depends on the
elementary update.

### Arguments
* xmin -> Minimum value of Œ¥Œæ.
* xmax -> Maximum value of Œ¥Œæ
* rng -> Random number generator.

### Returns
* N -> Value of Œ¥Œæ.
"""
function Pdx(xmin::F64, xmax::F64, rng::AbstractRNG)
    xmin_abs = abs(xmin)
    xmax_abs = abs(xmax)
    X = max(xmin_abs, xmax_abs)

    Œ≥ = 2.0
    Œ≥_X = Œ≥ / X

    Œ∑ = rand(rng, F64)
    ùëÅ  = (1 - Œ∑) * copysign(expm1(-Œ≥_X * xmin_abs), xmin)
    ùëÅ +=      Œ∑  * copysign(expm1(-Œ≥_X * xmax_abs), xmax)

    return copysign( log1p(-abs(ùëÅ)) / Œ≥_X, ùëÅ)
end
