#
# Project : Gardenia
# Source  : sac.jl
# Author  : Li Huang (huangli@caep.cn)
# Status  : Unstable
#
# Last modified: 2024/09/30
#

#=
### *Customized Structs* : *StochAC Solver*
=#

"""
    StochACElement

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

### Members
* Œì‚Çö -> It means the positions of the Œ¥ functions.
* Œì‚Çê -> It means the weights / amplitudes of the Œ¥ functions.
"""
mutable struct StochACElement
    Œì‚Çö :: Array{I64,2}
    Œì‚Çê :: Array{F64,2}
end

"""
    StochACContext

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

### Members
* G·µ•     -> Input data for correlator.
* œÉ¬π     -> Actually 1.0 / œÉ¬π.
* allow  -> Allowable indices.
* grid   -> Imaginary axis grid for input data.
* mesh   -> Real frequency mesh for output spectrum.
* model  -> Default model function.
* kernel -> Default kernel function.
* Aout   -> Calculated spectral function, it is actually ‚ü®n(x)‚ü©.
* Œî      -> Precomputed Œ¥ functions.
* hœÑ     -> Œ±-resolved h(œÑ).
* HŒ±     -> Œ±-resolved Hc.
* UŒ±     -> Œ±-resolved internal energy, it is actually ‚ü®HŒ±‚ü©.
* Œ±‚Çó     -> Vector of the Œ± parameters.
"""
mutable struct StochACContext
    G·µ•     :: Vector{F64}
    œÉ¬π     :: Vector{F64}
    allow  :: Vector{I64}
    grid   :: AbstractGrid
    mesh   :: AbstractMesh
    model  :: Vector{F64}
    kernel :: Array{F64,2}
    Aout   :: Array{F64,2}
    Œî      :: Array{F64,2}
    hœÑ     :: Array{F64,2}
    HŒ±     :: Vector{F64}
    UŒ±     :: Vector{F64}
    Œ±‚Çó     :: Vector{F64}
end

#=
### *Global Drivers*
=#

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

Solve the analytic continuation problem by the stochastic analytic
continuation algorithm (K. S. D. Beach's version). This is the driver for
the StochAC solver.

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

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

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

### Returns
* mesh -> Real frequency mesh, œâ.
* Asum -> Final spectral function, A(œâ). Note that it is Œ±-averaged.
* Gout -> Retarded Green's function, G(œâ).
"""
function solve(S::StochACSolver, rd::RawData)
    nmesh = get_b("nmesh")
    nalph = get_a("nalph")

    println("[ StochAC ]")
    MC, SE, SC = init(S, rd)

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

    return SC.mesh.mesh, Asum, Gout
end

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

Initialize the StochAC solver and return the StochACMC, StochACElement,
and StochACContext structs. Please don't call this function directly.

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

### Returns
* MC -> A StochACMC struct.
* SE -> A StochACElement struct.
* SC -> A StochACContext struct.
"""
function init(S::StochACSolver, rd::RawData)
    # Initialize possible constraints.
    # The array allow contains all the possible indices for Œ¥ functions.
    fmesh = calc_fmesh(S)
    allow = constraints(S, fmesh)

    # 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")

    # Initialize Monte Carlo configurations
    SE = init_element(S, MC.rng, allow)
    println("Randomize Monte Carlo configurations")

    # Prepare some key variables
    SC = init_context(SE, G·µ•, œÉ¬π, allow, grid, mesh, fmesh)
    println("Initialize context for the StochAC solver")

    return MC, SE, SC
end

"""
    run(MC::StochACMC, SE::StochACElement, SC::StochACContext)

Perform stochastic analytic continuation simulation, sequential version.

### Arguments
* MC -> A StochACMC struct.
* SE -> A StochACElement struct.
* SC -> A StochACContext struct.

### Returns
* Aout -> Spectral function, A(œâ).
* UŒ± -> Œ±-resolved internal energy.
"""
function run(MC::StochACMC, SE::StochACElement, SC::StochACContext)
    # 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
    nstep = get_a("nstep")
    output_per_steps = get_a("ndump")
    measure_per_steps = 100

    # Warmup the Monte Carlo engine
    println("Start thermalization...")
    warmup(MC, SE, SC)

    # Sample and collect data
    step = 0.0
    println("Start stochastic sampling...")
    for iter = 1:nstep
        sample(MC, SE, SC)

        if iter % measure_per_steps == 0
            step = step + 1.0
            measure(SE, SC)
        end

        if iter % output_per_steps == 0
            prog = round(I64, iter / nstep * 100)
            @printf("step = %9i ", iter)
            @printf("[progress = %3i]\n", prog)
            flush(stdout)
            fwrite && write_statistics(MC)
        end
    end

    return average(step, SC)
end

"""
    prun(
        S::StochACSolver,
        p1::Dict{String,Vector{Any}},
        p2::Dict{String,Vector{Any}},
        MC::StochACMC,
        SE::StochACElement,
        SC::StochACContext
    )

Perform stochastic analytic continuation simulation, parallel version.
The arguments `p1` and `p2` are copies of PBASE and PStochAC, respectively.

### Arguments
* S -> A StochACSolver struct.
* p1 -> A copy of PBASE.
* p2 -> A copy of PStochAC.
* MC -> A StochACMC struct.
* SE -> A StochACElement struct.
* SC -> A StochACContext struct.

### Returns
* Aout -> Spectral function, A(œâ).
* UŒ± -> Œ±-resolved internal energy.
"""
function prun(
    S::StochACSolver,
    p1::Dict{String,Vector{Any}},
    p2::Dict{String,Vector{Any}},
    MC::StochACMC,
    SE::StochACElement,
    SC::StochACContext
    )
    # Revise parameteric dicts
    rev_dict_b(p1)
    rev_dict_a(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
    nstep = get_a("nstep")
    output_per_steps = get_a("ndump")
    measure_per_steps = 100

    # Warmup the Monte Carlo engine
    println("Start thermalization...")
    warmup(MC, SE, SC)

    # Sample and collect data
    step = 0.0
    println("Start stochastic sampling...")
    for iter = 1:nstep
        sample(MC, SE, SC)

        if iter % measure_per_steps == 0
            step = step + 1.0
            measure(SE, SC)
        end

        if iter % output_per_steps == 0
            prog = round(I64, iter / nstep * 100)
            @printf("step = %9i ", iter)
            @printf("[progress = %3i]\n", prog)
            flush(stdout)
            myid() == 2 && fwrite && write_statistics(MC)
        end
    end

    return average(step, SC)
end

"""
    average(step::F64, SC::StochACContext)

Postprocess the results generated during the stochastic analytic
continuation simulations. It will calculate the spectral functions, and
Œ±-resolved internal energies.

### Arguments
* step -> How many steps are there in the Monte Carlo samplings.
* SC   -> A StochACContext struct.

### Returns
* Aout -> Spectral function, A(œâ,Œ±).
* UŒ± -> Œ±-resolved internal energy.
"""
function average(step::F64, SC::StochACContext)
    # Get key parameters
    nmesh = length(SC.mesh)
    nalph = length(SC.Œ±‚Çó)

    # Renormalize the spectral functions
    Aout = zeros(F64, nmesh, nalph)
    for i = 1:nalph
        for j = 1:nmesh
            Aout[j,i] = SC.Aout[j,i] * SC.model[j] / œÄ / step
        end
    end

    # Renormalize the internal energies
    UŒ± = SC.UŒ± / step

    return Aout, UŒ±
end

"""
    last(SC::StochACContext, Aout::Array{F64,2}, UŒ±::Vector{F64})

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

### Arguments
* SC   -> A StochACContext struct.
* Aout -> Œ±-dependent spectral functions.
* UŒ±   -> Œ±-dependent internal energies.

### Returns
* Asum -> Final spectral function (Œ±-averaged), A(œâ).
* G -> Retarded Green's function, G(œâ).
"""
function last(SC::StochACContext, Aout::Array{F64,2}, UŒ±::Vector{F64})
    function fitfun(x, p)
        return @. p[1] * x + p[2]
    end

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

    # Get dimensional parameters
    nmesh, nalph = size(Aout)

    # Try to fit the internal energies to find out optimal Œ±
    guess = [1.0, 1.0]
    fit_l = curve_fit(fitfun, SC.Œ±‚Çó[1:5], log10.(UŒ±[1:5]), guess)
    fit_r = curve_fit(fitfun, SC.Œ±‚Çó[end-4:end], log10.(UŒ±[end-4:end]), guess)
    a, b = fit_l.param
    c, d = fit_r.param
    aopt = (d - b) / (a - c)
    close = argmin( abs.( SC.Œ±‚Çó .- aopt ) )
    println("Fitting parameters [a,b] are: [ $a, $b ]")
    println("Fitting parameters [c,d] are: [ $c, $d ]")
    println("Perhaps the optimal Œ± is: ", aopt)
    fwrite && write_hamiltonian(SC.Œ±‚Çó, UŒ±)

    # Calculate final spectral functions and write them
    Asum = zeros(F64, nmesh)
    for i = close : nalph - 1
        @. Asum = Asum + (UŒ±[i] - UŒ±[i+1]) * Aout[:,i]
    end
    @. Asum = Asum / (UŒ±[close] - UŒ±[end])
    fwrite && write_spectrum(SC.mesh, Asum)
    fwrite && write_spectrum(SC.mesh, SC.Œ±‚Çó, Aout)
    fwrite && write_model(SC.mesh, SC.model)

    # Reproduce input data and write them
    kernel = make_kernel(SC.mesh, SC.grid)
    G = reprod(SC.mesh, kernel, Asum)
    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, Asum)
    else
        _G = kramers(SC.mesh, Asum .* SC.mesh)
    end
    fwrite && write_complete(SC.mesh, _G)

    return Asum, _G
end

#=
### *Core Algorithms*
=#

"""
    warmup(MC::StochACMC, SE::StochACElement, SC::StochACContext)

Warmup the Monte Carlo engine to acheieve thermalized equilibrium. After
that, the Monte Carlo counters will be reset.

### Arguments
* MC -> A StochACMC struct.
* SE -> A StochACElement struct.
* SC -> A StochACContext struct.

### Returns
N/A
"""
function warmup(MC::StochACMC, SE::StochACElement, SC::StochACContext)
    # Set essential parameter
    nwarm = get_a("nwarm")
    output_per_steps = 100

    # Shuffle the Monte Carlo field configuration
    for iter = 1:nwarm
        sample(MC, SE, SC)

        if iter % output_per_steps == 0
            prog = round(I64, iter / nwarm * 100)
            @printf("step = %9i ", iter)
            @printf("swap( 1 ) -> %8.4f ", MC.Macc[1] / MC.Mtry[1])
            @printf("swap(end) -> %8.4f ", MC.Macc[end] / MC.Mtry[end])
            @printf("[progress = %3i]\n", prog)
            flush(stdout)
        end
    end

    # Reset the counters
    fill!(MC.Macc, 0.0)
    fill!(MC.Mtry, 0.0)

    fill!(MC.Sacc, 0.0)
    fill!(MC.Stry, 0.0)
end

"""
    sample(MC::StochACMC, SE::StochACElement, SC::StochACContext)

Perform Monte Carlo sweeps and sample the field configurations.

### Arguments
* MC -> A StochACMC struct.
* SE -> A StochACElement struct.
* SC -> A StochACContext struct.

### Returns
N/A
"""
function sample(MC::StochACMC, SE::StochACElement, SC::StochACContext)
    nalph = get_a("nalph")

    if rand(MC.rng) < 0.9
        if rand(MC.rng) > 0.5
            for i = 1:nalph
                try_move_a(i, MC, SE, SC)
            end
        else
            if rand(MC.rng) > 0.2
                for i = 1:nalph
                    try_move_s(i, MC, SE, SC)
                end
            else
                for i = 1:nalph
                    try_move_p(i, MC, SE, SC)
                end
            end
        end
    else
        if nalph > 1
            try_move_x(MC, SE, SC)
        end
    end
end

"""
    measure(SE::StochACElement, SC::StochACContext)

Accumulate the Œ±-resolved spectral functions and internal energies.

### Arguments
* SE -> A StochACElement struct.
* SC -> A StochACContext struct.

### Returns
N/A
"""
function measure(SE::StochACElement, SC::StochACContext)
    nalph = get_a("nalph")

    # Loop over each Œ± parameter
    for ia = 1:nalph
        da = view(SE.Œì‚Çê, :, ia)
        dp = view(SE.Œì‚Çö, :, ia)
        SC.Aout[:,ia] = SC.Aout[:,ia] .+ SC.Œî[:,dp] * da
        SC.UŒ±[ia] = SC.UŒ±[ia] + SC.HŒ±[ia]
    end
end

#=
### *Service Functions*
=#

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

Preprocess the input data (`rd`).

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

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

See also: [`RawData`](@ref).
"""
function init_iodata(S::StochACSolver, rd::RawData)
    G = make_data(rd)
    G·µ• = G.value # G·µ• = abs.(G.value)
    œÉ¬π = 1.0 ./ sqrt.(G.covar)

    return G·µ•, œÉ¬π
end

"""
    init_mc(S::StochACSolver)

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

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

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

See also: [`StochACMC`](@ref).
"""
function init_mc(S::StochACSolver)
    nalph = get_a("nalph")
    #
    seed = rand(1:100000000)
    rng = MersenneTwister(seed)
    #
    Macc = zeros(F64, nalph)
    Mtry = zeros(F64, nalph)
    Sacc = zeros(F64, nalph)
    Stry = zeros(F64, nalph)
    #
    MC = StochACMC(rng, Macc, Mtry, Sacc, Stry)

    return MC
end

"""
    init_element(
        S::StochACSolver,
        rng::AbstractRNG,
        allow::Vector{I64}
    )

Randomize the configurations for future Monte Carlo sampling. It will
return a StochACElement struct.

### Arguments
* S     -> A StochACSolver struct.
* rng   -> Random number generator.
* allow -> Allowed positions for the Œ¥ peaks.

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

See also: [`StochACElement`](@ref).
"""
function init_element(
    S::StochACSolver,
    rng::AbstractRNG,
    allow::Vector{I64}
    )
    nalph = get_a("nalph")
    ngamm = get_a("ngamm")

    Œì‚Çö = rand(rng, allow, (ngamm, nalph))
    Œì‚Çê = rand(rng, F64, (ngamm, nalph))

    for j = 1:nalph
        Œì‚±º = view(Œì‚Çê, :, j)
        s = sum(Œì‚±º)
        @. Œì‚±º = Œì‚±º / s
    end

    SE = StochACElement(Œì‚Çö, Œì‚Çê)

    return SE
end

"""
    init_context(
        SE::StochACElement,
        G·µ•::Vector{F64},
        œÉ¬π::Vector{F64},
        allow::Vector{I64},
        grid::AbstractGrid,
        mesh::AbstractMesh,
        fmesh::AbstractMesh
    )

Try to create a StochACContext struct, which contains some key variables,
including grid, mesh, input correlator and the corresponding standard
deviation, kernel matrix, spectral function, and Œ±-resolved Hamiltonian.

### Arguments
* SE -> A StochACElement struct.
* G·µ• -> Input correlator. It will be changed in this function.
* œÉ¬π -> Standard deviation for input correlator.
* allow -> Allowable indices for Œ¥-like peaks.
* grid -> Imaginary axis grid for input data.
* mesh -> Real frequency mesh for output spectrum.
* fmesh -> Very fine mesh in [wmin, wmax].

### Returns
* SC -> A StochACContext struct.
"""
function init_context(
    SE::StochACElement,
    G·µ•::Vector{F64},
    œÉ¬π::Vector{F64},
    allow::Vector{I64},
    grid::AbstractGrid,
    mesh::AbstractMesh,
    fmesh::AbstractMesh
    )
    # Get parameters
    nmesh = get_b("nmesh")
    nalph = get_a("nalph")

    # Allocate memory for spectral function, A(œâ,Œ±)
    Aout = zeros(F64, nmesh, nalph)

    # Prepare some key variables
    # Only flat model is valid for the StochAC solver.
    model = make_model(mesh)

    # Precompute Œ¥ functions
    œï = calc_phi(mesh, model)
    Œî = calc_delta(fmesh, œï)

    # Build kernel matrix
    kernel = make_kernel(fmesh, grid)

    # In order to accelerate the calculations, the singular space of the
    # kernel function is used. At first, we preform singular value
    # decomposition for K/œÉ:
    #     K/œÉ = U S V·µÄ
    # Then
    #     (G - KA)/œÉ = G/œÉ - K/œÉA
    #                = UU'(G/œÉ - USV·µÄA)
    #                = U(U'G/œÉ - U'USV·µÄA)
    #                = U(U'G/œÉ - SV·µÄA)
    #                = U(G' - K'A)
    # In the StochAC solver, let G·µ• ‚Üí G', kernel ‚Üí K'. Then new œá¬≤ is
    # calculated by
    #     |G' - K'A|¬≤
    # instead of
    #     |G - KA|¬≤/œÉ¬≤

    # Singular value decomposition of K/œÉ
    U, V, S = make_singular_space(Diagonal(œÉ¬π) * kernel)

    # Get new kernel matrix
    kernel = Diagonal(S) * V'

    # Get new (input) correlator
    G·µ• = U' *  (G·µ• .* œÉ¬π)

    # Precompute hamiltonian
    hœÑ, HŒ±, UŒ± = calc_hamil(SE.Œì‚Çö, SE.Œì‚Çê, kernel, G·µ•)

    # Precompute Œ± parameters
    Œ±‚Çó = calc_alpha()

    return StochACContext(G·µ•, œÉ¬π, allow, grid, mesh, model,
                        kernel, Aout, Œî, hœÑ, HŒ±, UŒ±, Œ±‚Çó)
end

"""
    calc_fmesh(S::StochACSolver)

Try to calculate very fine (dense) mesh in [wmin, wmax], which is used
internally to build the kernel function. Note that this mesh could be
non-uniform. If the file `fmesh.inp` exists, the code will try to load
it to initialize the mesh. Or else the code will try to create a linear
mesh automatically.

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

### Returns
* fmesh -> A very fine, perhaps non-uniform mesh in [wmin, wmax].

See also: [`LinearMesh`](@ref), [`DynamicMesh`](@ref).
"""
function calc_fmesh(S::StochACSolver)
    wmin = get_b("wmin")
    wmax = get_b("wmax")
    nfine = get_a("nfine")

    # Filename for the predefined mesh
    # This file should contain at least `nfine` lines
    fn = "fmesh.inp"

    # If the file `fmesh.inp` exists, we will use it to build the mesh.
    if isfile(fn)
        mesh = zeros(F64, nfine)
        #
        open(fn, "r") do fin
            for i = 1:nfine
                arr = line_to_array(fin)
                mesh[i] = parse(F64, arr[2])
            end
        end
        #
        fmesh = DynamicMesh(mesh)
    # Or else we will return a linear mesh directly.
    else
        fmesh = LinearMesh(nfine, wmin, wmax)
    end

    return fmesh
end

"""
    calc_phi(am::AbstractMesh, model::Vector{F64})

Try to calculate œï(œâ) function. `am` is the mesh for calculated spectrum,
and `model` means the default model function.

For the definition of the œï(œâ) function, see Eq.(16) in arXiv:0403055. It
creates a smooth mapping from ‚Ñù to [0,1].

### Arguments
See above explanations.

### Returns
* œï -> The œï(œâ) function.

See also: [`calc_delta`](@ref).
"""
function calc_phi(am::AbstractMesh, model::Vector{F64})
    œï = cumsum(model .* am.weight)
    return œï
end

"""
    calc_delta(fmesh::AbstractMesh, œï::Vector{F64})

Precompute the Œî functions. `fmesh` is a very dense mesh in [wmin, wmax]
and `œï` is the œï function.

Here we just use f(x) = Œ∑ / (x¬≤ + Œ∑¬≤) to approximate the Œ¥ function, where
Œ∑ is a small parameter.

### Arguments
See above explanations.

### Returns
* Œî -> The Œî(œâ) function.

See also: [`calc_phi`](@ref).
"""
function calc_delta(fmesh::AbstractMesh, œï::Vector{F64})
    nmesh = length(œï)
    #
    nfine = length(fmesh)
    wmax = fmesh.wmax
    wmin = fmesh.wmin
    #
    Œ∑‚ÇÅ = 0.001
    Œ∑‚ÇÇ = 0.001 ^ 2.0

    Œî = zeros(F64, nmesh, nfine)
    s = similar(œï)
    for i = 1:nfine
        # We should convert the mesh `fmesh` from [wmin,wmax] to [0,1].
        ùë• = (fmesh[i] - wmin) / (wmax - wmin)
        @. s = (œï - ùë•) ^ 2.0 + Œ∑‚ÇÇ
        @. Œî[:,i] = Œ∑‚ÇÅ / s
    end

    return Œî
end

"""
    calc_hamil(
        Œì‚Çö::Array{I64,2},
        Œì‚Çê::Array{I64,2},
        kernel::Matrix{F64},
        G·µ•::Vector{F64}
    )

Initialize h(œÑ) and H(Œ±) using Eq.(35) and Eq.(36), respectively. `Œì‚Çö`
and `Œì‚Çê` represent n(x), `kernel` means the kernel function, `G·µ•` is the
correlator. Note that `kernel` and `G·µ•` have been rotated into singular
space. Please see comments in `init()` for more details.

### Arguments
See above explanations.

### Returns
* hœÑ -> Œ±-resolved h(œÑ).
* HŒ± -> Œ±-resolved Hc.
* UŒ± -> Œ±-resolved internal energy, it is actually ‚ü®HŒ±‚ü©.

See also: [`calc_htau`](@ref).
"""
function calc_hamil(
    Œì‚Çö::Array{I64,2},
    Œì‚Çê::Array{F64,2},
    kernel::Matrix{F64},
    G·µ•::Vector{F64}
    )
    nalph = get_a("nalph")
    ngrid = length(G·µ•) # It is not equal to get_b("ngrid") any more!

    hœÑ = zeros(F64, ngrid, nalph)
    HŒ± = zeros(F64, nalph)
    UŒ± = zeros(F64, nalph)

    for i = 1:nalph
        hœÑ[:,i] = calc_htau(Œì‚Çö[:,i], Œì‚Çê[:,i], kernel, G·µ•)
        HŒ±[i] = dot(hœÑ[:,i], hœÑ[:,i])
    end

    return hœÑ, HŒ±, UŒ±
end

"""
    calc_htau(
        Œì‚Çö::Vector{I64},
        Œì‚Çê::Vector{F64},
        kernel::Matrix{F64},
        G·µ•:Vector{F64}
    )

Try to calculate Œ±-dependent h(œÑ) via Eq.(36). `Œì‚Çö` and `Œì‚Çê` represent
n(x), `kernel` means the kernel function, `G·µ•` is the correlator. Note
that `kernel` and `G·µ•` have been rotated into singular space. Please
see comments in `init_context()` for more details.

### Arguments
See above explanations.

### Returns
* hœÑ -> Œ±-resolved h(œÑ).

See also: [`calc_hamil`](@ref).
"""
function calc_htau(Œì‚Çö::Vector{I64}, Œì‚Çê::Vector{F64},
                   kernel::Matrix{F64},
                   G·µ•::Vector{F64})
    hœÑ = similar(G·µ•)
    #
    for i in eachindex(G·µ•)
        hœÑ[i] = dot(Œì‚Çê, view(kernel, i, Œì‚Çö)) - G·µ•[i]
    end
    #
    return hœÑ
end

"""
    calc_alpha()

Generate a list for the Œ± parameters.

### Arguments
N/A

### Returns
* Œ±‚Çó -> List of the Œ± parameters.
"""
function calc_alpha()
    nalph = get_a("nalph")
    alpha = get_a("alpha")
    ratio = get_a("ratio")

    Œ±‚Çó = collect(alpha * (ratio ^ (x - 1)) for x in 1:nalph)

    return Œ±‚Çó
end

"""
    constraints(S::StochACSolver, fmesh::AbstractMesh)

Try to implement the constrained stochastic analytic continuation
method. This function will return a collection. It contains all the
allowable indices. Be careful, the constrained stochastic analytic
continuation method is compatible with the self-adaptive mesh.

### Arguments
* S     -> A StochACSolver struct.
* fmesh -> Very dense mesh for the Œ¥ peaks.

### Returns
* allow -> Allowable indices.

See also: [`StochACSolver`](@ref).
"""
function constraints(S::StochACSolver, fmesh::AbstractMesh)
    exclude = get_b("exclude")
    nfine = get_a("nfine")
    @assert nfine == length(fmesh)

    allow = I64[]

    # Go through the fine mesh and check every mesh point.
    # Is is excluded?
    for i in eachindex(fmesh)
        is_excluded = false
        #
        if !isa(exclude, Missing)
            for j in eachindex(exclude)
                if exclude[j][1] ‚â§ fmesh[i] ‚â§ exclude[j][2]
                    is_excluded = true
                    break
                end
            end
        end
        #
        if !is_excluded
            push!(allow, i)
        end
    end

    return allow
end

"""
    try_move_s(
        i::I64,
        MC::StochACMC,
        SE::StochACElement,
        SC::StochACContext
    )

Select one Œ¥ function randomly and then change its position.

### Arguments
* i -> Index for Œ± parameters.
* MC -> A StochACMC struct.
* SE -> A StochACElement struct.
* SC -> A StochACContext struct.

### Returns
N/A

See also: [`try_move_p`](@ref).
"""
function try_move_s(
    i::I64,
    MC::StochACMC,
    SE::StochACElement,
    SC::StochACContext
    )
    # Get current number of Œ¥ functions
    ngamm = get_a("ngamm")

    # Choose one Œ¥ function
    Œ≥ = rand(MC.rng, 1:ngamm)

    # Extract weight for the Œ¥ function
    a = SE.Œì‚Çê[Œ≥,i]

    # Choose new position for the Œ¥ function
    p = rand(MC.rng, SC.allow)

    # Try to calculate the change of Hc using Eq.~(42).
    hc = view(SC.hœÑ, :, i)
    K‚Çö = view(SC.kernel, :, p)
    K·µß = view(SC.kernel, :, SE.Œì‚Çö[Œ≥,i])
    #
    Œ¥hc = a * (K‚Çö - K·µß)
    Œ¥H = dot(Œ¥hc, 2.0 * hc + Œ¥hc)

    # Apply Metropolis algorithm
    MC.Mtry[i] = MC.Mtry[i] + 1
    if Œ¥H ‚â§ 0.0 || exp(-SC.Œ±‚Çó[i] * Œ¥H) > rand(MC.rng)
        # Update Monte Carlo configurations
        SE.Œì‚Çö[Œ≥,i] = p

        # Update h(œÑ)
        @. hc = hc + Œ¥hc

        # Update Hc
        SC.HŒ±[i] = SC.HŒ±[i] + Œ¥H

        # Update Monte Carlo counter
        MC.Macc[i] = MC.Macc[i] + 1
    end
end

"""
    try_move_p(
        i::I64,
        MC::StochACMC,
        SE::StochACElement,
        SC::StochACContext
    )

Select two Œ¥ functions randomly and then change their positions.

### Arguments
* i -> Index for Œ± parameters.
* MC -> A StochACMC struct.
* SE -> A StochACElement struct.
* SC -> A StochACContext struct.

### Returns
N/A

See also: [`try_move_s`](@ref).
"""
function try_move_p(
    i::I64,
    MC::StochACMC,
    SE::StochACElement,
    SC::StochACContext
    )
    # Get current number of Œ¥ functions
    ngamm = get_a("ngamm")
    #
    if ngamm < 2
        return
    end

    # Choose two Œ¥ functions, they are labelled as Œ≥‚ÇÅ and Œ≥‚ÇÇ, respectively.
    Œ≥‚ÇÅ = 1
    Œ≥‚ÇÇ = 1
    while Œ≥‚ÇÅ == Œ≥‚ÇÇ
        Œ≥‚ÇÅ = rand(MC.rng, 1:ngamm)
        Œ≥‚ÇÇ = rand(MC.rng, 1:ngamm)
    end

    # Extract weights for the two Œ¥ functions (a‚ÇÅ and a‚ÇÇ)
    a‚ÇÅ = SE.Œì‚Çê[Œ≥‚ÇÅ,i]
    a‚ÇÇ = SE.Œì‚Çê[Œ≥‚ÇÇ,i]

    # Choose new positions for the two Œ¥ functions (p‚ÇÅ and p‚ÇÇ).
    # Note that their old positions are SE.Œì‚Çö[Œ≥‚ÇÅ,i] and SE.Œì‚Çö[Œ≥‚ÇÇ,i].
    p‚ÇÅ = rand(MC.rng, SC.allow)
    p‚ÇÇ = rand(MC.rng, SC.allow)

    # Try to calculate the change of Hc using Eq.~(42).
    hc = view(SC.hœÑ, :, i)
    K‚ÇÅ = view(SC.kernel, :, p‚ÇÅ)
    K‚ÇÇ = view(SC.kernel, :, p‚ÇÇ)
    K‚ÇÉ = view(SC.kernel, :, SE.Œì‚Çö[Œ≥‚ÇÅ,i])
    K‚ÇÑ = view(SC.kernel, :, SE.Œì‚Çö[Œ≥‚ÇÇ,i])
    #
    Œ¥hc = a‚ÇÅ * (K‚ÇÅ - K‚ÇÉ) + a‚ÇÇ * (K‚ÇÇ - K‚ÇÑ)
    Œ¥H = dot(Œ¥hc, 2.0 * hc + Œ¥hc)

    # Apply Metropolis algorithm
    MC.Mtry[i] = MC.Mtry[i] + 1
    if Œ¥H ‚â§ 0.0 || exp(-SC.Œ±‚Çó[i] * Œ¥H) > rand(MC.rng)
        # Update Monte Carlo configurations
        SE.Œì‚Çö[Œ≥‚ÇÅ,i] = p‚ÇÅ
        SE.Œì‚Çö[Œ≥‚ÇÇ,i] = p‚ÇÇ

        # Update h(œÑ)
        @. hc = hc + Œ¥hc

        # Update Hc
        SC.HŒ±[i] = SC.HŒ±[i] + Œ¥H

        # Update Monte Carlo counter
        MC.Macc[i] = MC.Macc[i] + 1
    end
end

"""
    try_move_a(
        i::I64,
        MC::StochACMC,
        SE::StochACElement,
        SC::StochACContext
    )

Select two Œ¥ functions randomly and then change their weights.

### Arguments
* i -> Index for Œ± parameters.
* MC -> A StochACMC struct.
* SE -> A StochACElement struct.
* SC -> A StochACContext struct.

### Returns
N/A

See also: [`try_move_x`](@ref).
"""
function try_move_a(
    i::I64,
    MC::StochACMC,
    SE::StochACElement,
    SC::StochACContext
    )
    # Get current number of Œ¥ functions
    ngamm = get_a("ngamm")
    #
    if ngamm < 2
        return
    end

    # Choose two Œ¥ functions, they are labelled as Œ≥‚ÇÅ and Œ≥‚ÇÇ, respectively.
    Œ≥‚ÇÅ = 1
    Œ≥‚ÇÇ = 1
    while Œ≥‚ÇÅ == Œ≥‚ÇÇ
        Œ≥‚ÇÅ = rand(MC.rng, 1:ngamm)
        Œ≥‚ÇÇ = rand(MC.rng, 1:ngamm)
    end

    # Extract weights for the two Œ¥ functions (a‚ÇÉ and a‚ÇÑ), then try to
    # calculate new weights for them (a‚ÇÅ and a‚ÇÇ).
    a‚ÇÅ = 0.0
    a‚ÇÇ = 0.0
    a‚ÇÉ = SE.Œì‚Çê[Œ≥‚ÇÅ,i]
    a‚ÇÑ = SE.Œì‚Çê[Œ≥‚ÇÇ,i]
    Œ¥a = 0.0
    while true
        Œ¥a = rand(MC.rng) * (a‚ÇÉ + a‚ÇÑ) - a‚ÇÉ
        a‚ÇÅ = a‚ÇÉ + Œ¥a
        a‚ÇÇ = a‚ÇÑ - Œ¥a
        if a‚ÇÅ > 0 && a‚ÇÇ > 0
            break
        end
    end

    # Try to calculate the change of Hc using Eq.~(42).
    hc = view(SC.hœÑ, :, i)
    K‚ÇÅ = view(SC.kernel, :, SE.Œì‚Çö[Œ≥‚ÇÅ,i])
    K‚ÇÇ = view(SC.kernel, :, SE.Œì‚Çö[Œ≥‚ÇÇ,i])
    #
    Œ¥hc = Œ¥a * (K‚ÇÅ - K‚ÇÇ)
    Œ¥H = dot(Œ¥hc, 2.0 * hc + Œ¥hc)

    # Apply Metropolis algorithm
    MC.Mtry[i] = MC.Mtry[i] + 1
    if Œ¥H ‚â§ 0.0 || exp(-SC.Œ±‚Çó[i] * Œ¥H) > rand(MC.rng)
        # Update Monte Carlo configurations
        SE.Œì‚Çê[Œ≥‚ÇÅ,i] = a‚ÇÅ
        SE.Œì‚Çê[Œ≥‚ÇÇ,i] = a‚ÇÇ

        # Update h(œÑ)
        @. hc = hc + Œ¥hc

        # Update Hc
        SC.HŒ±[i] = SC.HŒ±[i] + Œ¥H

        # Update Monte Carlo counter
        MC.Macc[i] = MC.Macc[i] + 1
    end
end

"""
    try_move_x(
        MC::StochACMC,
        SE::StochACElement,
        SC::StochACContext
    )

Try to exchange field configurations between two adjacent layers. Because
this function involves two layers, so it doesn't need the argument `i`.

### Arguments
* MC -> A StochACMC struct.
* SE -> A StochACElement struct.
* SC -> A StochACContext struct.

### Returns
N/A

See also: [`try_move_a`](@ref).
"""
function try_move_x(
    MC::StochACMC,
    SE::StochACElement,
    SC::StochACContext
    )
    # Get number of Œ± parameters
    nalph = get_a("nalph")

    # Select two adjacent layers (two adjacent Œ± parameters)
    i = rand(MC.rng, 1:nalph)
    j = rand(MC.rng) > 0.5 ? i + 1 : i - 1
    i == 1 && (j = i + 1)
    i == nalph && (j = i - 1)

    # Calculate change of Hc
    Œ¥Œ± = SC.Œ±‚Çó[i] - SC.Œ±‚Çó[j]
    Œ¥H = SC.HŒ±[i] - SC.HŒ±[j]

    # Apply Metropolis algorithm
    MC.Stry[i] = MC.Stry[i] + 1
    MC.Stry[j] = MC.Stry[j] + 1
    if exp(Œ¥Œ± * Œ¥H) > rand(MC.rng)
        # Update Monte Carlo configurations
        SE.Œì‚Çö[:,i], SE.Œì‚Çö[:,j] = SE.Œì‚Çö[:,j], SE.Œì‚Çö[:,i]
        SE.Œì‚Çê[:,i], SE.Œì‚Çê[:,j] = SE.Œì‚Çê[:,j], SE.Œì‚Çê[:,i]

        # Update h(œÑ) and Hc
        SC.hœÑ[:,i], SC.hœÑ[:,j] = SC.hœÑ[:,j], SC.hœÑ[:,i]
        SC.HŒ±[i], SC.HŒ±[j] = SC.HŒ±[j], SC.HŒ±[i]

        # Update Monte Carlo counters
        MC.Sacc[i] = MC.Sacc[i] + 1
        MC.Sacc[j] = MC.Sacc[j] + 1
    end
end
