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

#=
### *Customized Structs* : *StochPX Solver*
=#

"""
    StochPXElement

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

For the off-diagonal elements of the matrix-valued Green's function, the
signs of the poles (ğ�•Š) could be negative (-1.0). However, for the other
cases, ğ�•Š is always positive (+1.0).

### Members
* P  -> It means the positions of the poles.
* A  -> It means the weights / amplitudes of the poles.
* ğ�•Š  -> It means the signs of the poles.
"""
mutable struct StochPXElement
    P  :: Vector{I64}
    A  :: Vector{F64}
    ğ�•Š  :: Vector{F64}
end

"""
    StochPXContext

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

Note that Ï‡Â² denotes the goodness-of-fit functional, and Gáµ§ denotes the
reproduced correlator. They should be always compatible with P, A, and ğ�•Š
in the StochPXElement struct.

### Members
* Gáµ¥    -> Input data for correlator.
* Gáµ§    -> Generated correlator.
* ÏƒÂ¹    -> Actually 1.0 / ÏƒÂ¹.
* allow -> Allowable indices.
* grid  -> Grid for input data.
* mesh  -> Mesh for output spectrum.
* fmesh -> Very dense mesh for the poles.
* Î›     -> Precomputed kernel matrix.
* Î˜     -> Artificial inverse temperature.
* Ï‡Â²    -> Goodness-of-fit functional for the current configuration.
* Ï‡Â²áµ¥   -> Vector of goodness-of-fit functional.
* Páµ¥    -> Vector of poles' positions.
* Aáµ¥    -> Vector of poles' amplitudes.
* ğ�•Šáµ¥    -> Vector of poles' signs.
"""
mutable struct StochPXContext
    Gáµ¥    :: Vector{F64}
    Gáµ§    :: Vector{F64}
    ÏƒÂ¹    :: Vector{F64}
    allow :: Vector{I64}
    grid  :: AbstractGrid
    mesh  :: AbstractMesh
    fmesh :: AbstractMesh
    Î›     :: Array{F64,2}
    Î˜     :: F64
    Ï‡Â²    :: F64
    Ï‡Â²áµ¥   :: Vector{F64}
    Páµ¥    :: Vector{Vector{I64}}
    Aáµ¥    :: Vector{Vector{F64}}
    ğ�•Šáµ¥    :: Vector{Vector{F64}}
end

#=
### *Global Drivers*
=#

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

Solve the analytic continuation problem by the stochastic pole expansion.
It is the driver for the StochPX solver. Note that this solver is still
`experimental`. It is useful for analytic continuation of Matsubara data.
In other words, this solver should not be used to handle the imagnary-time
Green's function.

Similar to the BarRat and NevanAC solvers, this solver always returns A(Ï‰).

### Arguments
* S -> A StochPXSolver 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::StochPXSolver, rd::RawData)
    ktype = get_b("ktype")
    ngrid = get_b("ngrid")
    nmesh = get_b("nmesh")

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

    # Parallel version
    if nworkers() > 1
        #
        println("Using $(nworkers()) workers")
        #
        # Copy configuration dicts
        p1 = deepcopy(PBASE)
        p2 = deepcopy(PStochPX)
        #
        # 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)
        Gout = zeros(C64, nmesh)
        if ktype == "bsymm"
            Gáµ£ = zeros(F64, ngrid)
        else
            Gáµ£ = zeros(F64, 2 * ngrid)
        end
        for i in eachindex(sol)
            a, b, c = sol[i]
            @. Aout = Aout + a / nworkers()
            @. Gout = Gout + b / nworkers()
            @. Gáµ£ = Gáµ£ + c / nworkers()
        end
        #
        # Postprocess the solutions
        last(SC, Aout, Gout, Gáµ£)
        #
    # Sequential version
    else
        #
        Aout, Gout, Gáµ£ = run(MC, SE, SC)
        last(SC, Aout, Gout, Gáµ£)
        #
    end

    return SC.mesh.mesh, Aout, Gout
end

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

Initialize the StochPX solver and return the StochPXMC, StochPXElement,
and StochPXContext structs. Please don't call this function directly.

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

### Returns
* MC -> A StochPXMC struct.
* SE -> A StochPXElement struct.
* SC -> A StochPXContext struct.
"""
function init(S::StochPXSolver, rd::RawData)
    # Initialize possible constraints.
    # The array arrow contains all the possible indices for poles.
    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
    Gáµ§, Î›, Î˜, Ï‡Â², Ï‡Â²áµ¥, Páµ¥, Aáµ¥, ğ�•Šáµ¥ = init_context(S, SE, grid, fmesh, Gáµ¥)
    SC = StochPXContext(Gáµ¥, Gáµ§, ÏƒÂ¹, allow, grid, mesh, fmesh,
                        Î›, Î˜, Ï‡Â², Ï‡Â²áµ¥, Páµ¥, Aáµ¥, ğ�•Šáµ¥)
    println("Initialize context for the StochPX solver")

    return MC, SE, SC
end

"""
    run(MC::StochPXMC, SE::StochPXElement, SC::StochPXContext)

Perform stochastic pole expansion simulation, sequential version.

### Arguments
* MC -> A StochPXMC struct.
* SE -> A StochPXElement struct.
* SC -> A StochPXContext struct.

### Returns
* Aout -> Spectral function, A(Ï‰).
* Gout -> Retarded Green's function, G(Ï‰).
* Gáµ£ -> Reproduced Green's function, G(iÏ‰â‚™).
"""
function run(MC::StochPXMC, SE::StochPXElement, SC::StochPXContext)
    # 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_x("ntry")
    nstep = get_x("nstep")

    # Warmup the Monte Carlo engine
    println("Start thermalization...")
    for _ = 1:nstep
        sample(1, MC, SE, SC)
    end

    # Sample and collect data
    println("Start stochastic sampling...")
    for t = 1:ntry
        # Reset Monte Carlo counters
        reset_mc(MC)

        # Reset Monte Carlo field configuration
        reset_element(MC.rng, SC.allow, SE)

        # Reset Gáµ§ and Ï‡Â² in SC (StochPXContext)
        reset_context(t, SE, SC)

        # Apply simulated annealing algorithm
        for _ = 1:nstep
            sample(t, MC, SE, SC)
        end

        # Write Monte Carlo statistics
        fwrite && write_statistics(MC)

        # Show the best Ï‡Â² (the smallest) for the current attempt
        @printf("try = %6i -> [Ï‡Â² = %9.4e]\n", t, SC.Ï‡Â²áµ¥[t])
        flush(stdout)
    end

    # Write pole expansion coefficients
    fwrite && write_pole(SC.Páµ¥, SC.Aáµ¥, SC.ğ�•Šáµ¥, SC.Ï‡Â²áµ¥, SC.fmesh)

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

"""
    prun(
        S::StochPXSolver,
        p1::Dict{String,Vector{Any}},
        p2::Dict{String,Vector{Any}},
        MC::StochPXMC,
        SE::StochPXElement,
        SC::StochPXContext
    )

Perform stochastic pole expansion simulation, parallel version.
The arguments `p1` and `p2` are copies of PBASE and PStochPX, respectively.

### Arguments
* S -> A StochPXSolver struct.
* p1 -> A copy of PBASE.
* p2 -> A copy of PStochPX.
* MC -> A StochPXMC struct.
* SE -> A StochPXElement struct.
* SC -> A StochPXContext struct.

### Returns
* Aout -> Spectral function, A(Ï‰).
* Gout -> Retarded Green's function, G(Ï‰).
* Gáµ£ -> Reproduced Green's function, G(iÏ‰â‚™).
"""
function prun(
    S::StochPXSolver,
    p1::Dict{String,Vector{Any}},
    p2::Dict{String,Vector{Any}},
    MC::StochPXMC,
    SE::StochPXElement,
    SC::StochPXContext
    )
    # Revise parameteric dicts
    # We have to make sure that all processes share the same parameters.
    rev_dict_b(p1)
    rev_dict_x(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_x("ntry")
    nstep = get_x("nstep")

    # Warmup the Monte Carlo engine
    println("Start thermalization...")
    for _ = 1:nstep
        sample(1, MC, SE, SC)
    end

    # Sample and collect data
    println("Start stochastic sampling...")
    for t = 1:ntry
        # Reset Monte Carlo counters
        reset_mc(MC)

        # Reset Monte Carlo field configuration
        reset_element(MC.rng, SC.allow, SE)

        # Reset Gáµ§ and Ï‡Â² in SC (StochPXContext)
        reset_context(t, SE, SC)

        # Apply simulated annealing algorithm
        for _ = 1:nstep
            sample(t, MC, SE, SC)
        end

        # Write Monte Carlo statistics
        myid() == 2 && fwrite && write_statistics(MC)

        # Show the best Ï‡Â² (the smallest) for the current attempt
        @printf("try = %6i -> [Ï‡Â² = %9.4e]\n", t, SC.Ï‡Â²áµ¥[t])
        flush(stdout)
    end

    # Write pole expansion coefficients
    myid() == 2 && fwrite && write_pole(SC.Páµ¥, SC.Aáµ¥, SC.ğ�•Šáµ¥, SC.Ï‡Â²áµ¥, SC.fmesh)

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

"""
    average(SC::StochPXContext)

Postprocess the results generated during the stochastic pole expansion
simulations. It will calculate the spectral functions, real frequency
Green's function, and imaginary frequency Green's function.

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

### Returns
* Aout -> Spectral function, A(Ï‰).
* Gout -> Retarded Green's function, G(Ï‰).
* Gáµ£ -> Reproduced Green's function, G(iÏ‰â‚™).
"""
function average(SC::StochPXContext)
    # 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
    nmesh = get_b("nmesh")
    method = get_x("method")
    ntry = get_x("ntry")

    # Allocate memory
    # Gout: real frequency Green's function, G(Ï‰).
    # Gáµ£: imaginary frequency Green's function, G(iÏ‰â‚™)
    ngrid, _ = size(SC.Î›)
    Gout = zeros(C64, nmesh)
    Gáµ£ = zeros(F64, ngrid)

    # Choose the best solution
    if method == "best"
        # The Ï‡Â² of the best solution should be the smallest.
        p = argmin(SC.Ï‡Â²áµ¥)
        @printf("Best solution: try = %6i -> [Ï‡Â² = %9.4e]\n", p, SC.Ï‡Â²áµ¥[p])
        #
        # Calculate G(Ï‰)
        Gout = calc_green(p, SC, true)
        #
        # Calculate G(iÏ‰â‚™)
        Gáµ£ = calc_green(p, SC, false)
    #
    # Collect the `good` solutions and calculate their average.
    else
        # Calculate the median of SC.Ï‡Â²áµ¥
        chi2_med = median(SC.Ï‡Â²áµ¥)
        chi2_ave = mean(SC.Ï‡Â²áµ¥)

        # Determine the Î±good parameter, which is used to filter the
        # calculated spectra.
        Î±good = 1.2
        if count(x -> x < chi2_med / Î±good, SC.Ï‡Â²áµ¥) â‰¤ ntry / 10
            Î±good = 1.0
        end

        # Go through all the solutions
        c = 0.0 # A counter
        passed = I64[]
        for i = 1:ntry
            if SC.Ï‡Â²áµ¥[i] < chi2_med / Î±good
                # Calculate and accumulate G(Ï‰)
                G = calc_green(i, SC, true)
                @. Gout = Gout + G
                #
                # Calculate and accumulate G(iÏ‰â‚™)
                G = calc_green(i, SC, false)
                @. Gáµ£ = Gáµ£ + G
                #
                # Increase the counter
                c = c + 1.0
                append!(passed, i)
            end
        end
        #
        # Normalize the final results
        @. Gout = Gout / c
        @. Gáµ£ = Gáµ£ / c
        println("Mean value of Ï‡Â²: $(chi2_ave)")
        println("Median value of Ï‡Â²: $(chi2_med)")
        println("Accumulate $(round(I64,c)) solutions to get the spectral density")
        #
        # Write indices of selected solutions
        if nworkers() > 1
            myid() == 2 && fwrite && write_passed(passed, chi2_med, Î±good)
        else
            fwrite && write_passed(passed, chi2_med, Î±good)
        end
    #
    end

    return -imag.(Gout) / Ï€, Gout, Gáµ£
end

"""
    last(
        SC::StochPXContext,
        Aout::Vector{F64},
        Gout::Vector{C64},
        Gáµ£::Vector{F64}
    )

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

### Arguments
* SC   -> A StochPXContext struct.
* Aout -> Spectral function, A(Ï‰).
* Gout -> Retarded Green's function, G(Ï‰).
* Gáµ£   -> Reconstructed Green's function, G(iÏ‰â‚™).

### Returns
N/A
"""
function last(
    SC::StochPXContext,
    Aout::Vector{F64},
    Gout::Vector{C64},
    Gáµ£::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, A(Ï‰).
    fwrite && write_spectrum(SC.mesh, Aout)

    # Reproduce input data and write them, G(iÏ‰â‚™).
    fwrite && write_backward(SC.grid, Gáµ£)

    # Write full response function on real axis, G(Ï‰).
    fwrite && write_complete(SC.mesh, Gout)
end

#=
### *Core Algorithms*
=#

"""
    sample(
        t::I64,
        MC::StochPXMC,
        SE::StochPXElement,
        SC::StochPXContext
    )

Try to search the configuration space to locate the minimum by using the
simulated annealing algorithm. Here, `t` means the t-th attempt.

### Arguments
* t -> Counter for attempts.
* MC -> A StochPXMC struct.
* SE -> A StochPXElement struct.
* SC -> A StochPXContext struct.

### Returns
N/A
"""
function sample(
    t::I64,
    MC::StochPXMC,
    SE::StochPXElement,
    SC::StochPXContext
    )
    # Try to change positions of poles
    if rand(MC.rng) < 0.5
        if rand(MC.rng) < 0.9
            try_move_s(t, MC, SE, SC)
        else
            try_move_p(t, MC, SE, SC)
        end
    # Try to change amplitudes of poles
    else
        if rand(MC.rng) < 0.5
            try_move_a(t, MC, SE, SC)
        else
            try_move_x(t, MC, SE, SC)
        end
    end
end

"""
    measure(t::I64, SE::StochPXElement, SC::StochPXContext)

Store Monte Carlo field configurations (positions, amplitudes, and signs
of many poles) for the `t`-th attempt. In other words, the current field
configuration (recorded in `SE`) will be saved in `SC`.

Note that not all configurations for the `t`-th attempt will be saved.
Only the solution that exhibits the smallest Ï‡Â² will be saved. For the
`t`-th attempt, the StochPX solver will do `nstep` Monte Carlo updates.
It will calculate all Ï‡Â², and try to figure out the smallest one. Then
the corresponding configuration (solution) will be saved.

### Arguments
* t -> Counter for the attemps.
* SE -> A StochPXElement struct.
* SC -> A StochPXContext struct.

### Returns
N/A
"""
function measure(t::I64, SE::StochPXElement, SC::StochPXContext)
    SC.Ï‡Â²áµ¥[t] = SC.Ï‡Â²
    #
    @. SC.Páµ¥[t] = SE.P
    @. SC.Aáµ¥[t] = SE.A
    @. SC.ğ�•Šáµ¥[t] = SE.ğ�•Š
end

#=
### *Service Functions*
=#

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

Preprocess the input data (`rd`).

### Arguments
* S -> A StochPXSolver 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::StochPXSolver, rd::RawData)
    G = make_data(rd)
    Gáµ¥ = G.value # Gáµ¥ = abs.(G.value)
    ÏƒÂ¹ = 1.0 ./ sqrt.(G.covar)

    return Gáµ¥, ÏƒÂ¹
end

"""
    init_mc(S::StochPXSolver)

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

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

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

See also: [`StochPXMC`](@ref).
"""
function init_mc(S::StochPXSolver)
    seed = rand(1:100000000)
    rng = MersenneTwister(seed)
    #
    Sacc = 0
    Stry = 0
    #
    Pacc = 0
    Ptry = 0
    #
    Aacc = 0
    Atry = 0
    #
    Xacc = 0
    Xtry = 0
    #
    MC = StochPXMC(rng, Sacc, Stry, Pacc, Ptry, Aacc, Atry, Xacc, Xtry)

    return MC
end

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

Randomize the configurations for future Monte Carlo sampling. It will
return a StochPXElement struct. Note that `allow` is generated in the
`constraints()` function.

### Arguments
* S     -> A StochPXSolver struct.
* rng   -> Random number generator.
* allow -> Allowed positions for the poles.

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

See also: [`StochPXElement`](@ref).
"""
function init_element(
    S::StochPXSolver,
    rng::AbstractRNG,
    allow::Vector{I64}
    )
    offdiag = get_b("offdiag")
    npole = get_x("npole")

    if offdiag
        # We just assume that the numbers of poles for the positive and
        # negative parts are equal.
        @assert iseven(npole)
        hpole = npole Ã· 2

        # We have to make sure that both negative weights and positive
        # weights are allowed for the generated poles.
        @assert any(x -> x > 0.0, allow)
        @assert any(x -> x < 0.0, allow)

        # Initialize P, A, and ğ�•Š
        allowâ‚Š = filter(x -> x > 0.0, allow)
        allowâ‚‹ = filter(x -> x < 0.0, allow)
        Pâ‚Š = rand(rng, allowâ‚Š, hpole)
        Pâ‚‹ = rand(rng, allowâ‚‹, hpole)
        Aâ‚Š = rand(rng, F64, hpole)
        Aâ‚‹ = rand(rng, F64, hpole)
        ğ�•Šâ‚Š = ones(F64, hpole)
        ğ�•Šâ‚‹ = ones(F64, hpole) * (-1.0)

        # We have to make sure âˆ‘áµ¢ Aáµ¢ = 1
        s = sum(Aâ‚Š)
        @. Aâ‚Š = Aâ‚Š / s
        s = sum(Aâ‚‹)
        @. Aâ‚‹ = Aâ‚‹ / s

        # Merge positive and negative parts
        P = vcat(Pâ‚Š, Pâ‚‹)
        A = vcat(Aâ‚Š, Aâ‚‹)
        ğ�•Š = vcat(ğ�•Šâ‚Š, ğ�•Šâ‚‹)
    else
        # Initialize P, A, and ğ�•Š
        P = rand(rng, allow, npole)
        A = rand(rng, F64, npole)
        ğ�•Š = ones(F64, npole)

        # We have to make sure âˆ‘áµ¢ Aáµ¢ = 1
        s = sum(A)
        @. A = A / s
    end

    SE = StochPXElement(abs.(P), A, ğ�•Š)

    return SE
end

"""
    init_context(
        S::StochPXSolver,
        SE::StochPXElement,
        grid::AbstractGrid,
        fmesh::AbstractMesh,
        Gáµ¥::Vector{F64}
    )

Try to initialize the key members of a StochPXContext struct. It will try
to return some key variables, which should be used to construct the
StochPXContext struct.

### Arguments
* S     -> A StochPXSolver struct.
* SE    -> A StochPXElement struct.
* grid  -> Grid for input correlator.
* fmesh -> Fine mesh in [wmin, wmax], used to build the kernel matrix Î›.
* Gáµ¥    -> Preprocessed input correlator.

### Returns
* Gáµ§  -> Reconstructed correlator.
* Î›   -> Precomputed kernel matrix.
* Î˜   -> Artificial inverse temperature.
* Ï‡Â²  -> Current goodness-of-fit functional.
* Ï‡Â²áµ¥ -> Vector of goodness-of-fit functional.
* Páµ¥  -> Vector of poles' positions.
* Aáµ¥  -> Vector of poles' amplitudes.
* ğ�•Šáµ¥  -> Vector of poles' signs.

See also: [`StochPXContext`](@ref).
"""
function init_context(
    S::StochPXSolver,
    SE::StochPXElement,
    grid::AbstractGrid,
    fmesh::AbstractMesh,
    Gáµ¥::Vector{F64}
    )
    # Extract some parameters
    ntry = get_x("ntry")
    npole = get_x("npole")
    Î˜ = get_x("theta")

    # Prepare the kernel matrix Î›. It is used to speed up the simulation.
    # Note that Î› depends on the type of kernel.
    Î› = calc_lambda(grid, fmesh, Gáµ¥)

    # We have to make sure that the starting Gáµ§ and Ï‡Â² are consistent with
    # the current Monte Carlo configuration fields.
    Gáµ§ = calc_green(SE.P, SE.A, SE.ğ�•Š, Î›)
    Ï‡Â² = calc_chi2(Gáµ§, Gáµ¥)

    # Ï‡Â²áµ¥ is initialized by a large number. Later it will be updated by
    # the smallest Ï‡Â² during the simulation.
    Ï‡Â²áµ¥ = zeros(F64, ntry)
    @. Ï‡Â²áµ¥ = 1e10

    # P, A, and ğ�•Š should be always compatible with Ï‡Â². They are updated
    # in the `measure()` function.
    Páµ¥ = Vector{I64}[]
    Aáµ¥ = Vector{F64}[]
    ğ�•Šáµ¥ = Vector{F64}[]
    #
    for _ = 1:ntry
        push!(Páµ¥,  ones(I64, npole))
        push!(Aáµ¥, zeros(F64, npole))
        push!(ğ�•Šáµ¥, zeros(F64, npole))
    end

    return Gáµ§, Î›, Î˜, Ï‡Â², Ï‡Â²áµ¥, Páµ¥, Aáµ¥, ğ�•Šáµ¥
end

"""
    reset_mc(MC::StochPXMC)

Reset the counters in StochPXMC struct.

### Arguments
* MC -> A StochPXMC struct.

### Returns
N/A
"""
function reset_mc(MC::StochPXMC)
    MC.Sacc = 0
    MC.Stry = 0
    #
    MC.Pacc = 0
    MC.Ptry = 0
    #
    MC.Aacc = 0
    MC.Atry = 0
    #
    MC.Xacc = 0
    MC.Xtry = 0
end

"""
    reset_element(
        rng::AbstractRNG,
        allow::Vector{I64},
        SE::StochPXElement
    )

Reset the Monte Carlo field configurations (i.e. positions and amplitudes
of the poles). Note that the signs of the poles should not be changed.
Be careful, the corresponding Ï‡Â² (goodness-of-fit functional) and Gáµ§
(generated correlator) will be updated in the `reset_context()` function.

### Arguments
* rng   -> Random number generator.
* allow -> Allowed positions for the poles.
* SE -> A StochPXElement struct.

### Returns
N/A
"""
function reset_element(
    rng::AbstractRNG,
    allow::Vector{I64},
    SE::StochPXElement
    )
    offdiag = get_b("offdiag")
    npole = get_x("npole")

    # For off-diagonal elements
    if offdiag
        # The number of poles must be even.
        # Here, `hpole` means half number of poles.
        hpole = npole Ã· 2

        # How many poles should be changed
        if hpole â‰¤ 5
            if 4 â‰¤ hpole â‰¤ 5
                hselect = 2
            else
                hselect = 1
            end
        else
            hselect = hpole Ã· 5
        end
        @assert hselect â‰¤ hpole

        # Which poles should be changed
        # For poles that with positive weights
        selectedâ‚Š = rand(rng, 1:hpole, hselect)
        unique!(selectedâ‚Š)
        hselectâ‚Š = length(selectedâ‚Š)
        #
        # For poles that with negative weights
        selectedâ‚‹ = rand(rng, hpole+1:npole, hselect)
        unique!(selectedâ‚‹)
        hselectâ‚‹ = length(selectedâ‚‹)

        # Change poles' positions
        if rand(rng) < 0.9
            allowâ‚Š = filter(x -> x > 0.0, allow)
            Pâ‚Š = rand(rng, allowâ‚Š, hselectâ‚Š)
            @. SE.P[selectedâ‚Š] = abs.(Pâ‚Š)
            #
            allowâ‚‹ = filter(x -> x < 0.0, allow)
            Pâ‚‹ = rand(rng, allowâ‚‹, hselectâ‚‹)
            @. SE.P[selectedâ‚‹] = abs.(Pâ‚‹)
        # Change poles' amplitudes
        else
            # For positive-weight poles
            Aâ‚�â‚Š = SE.A[selectedâ‚Š]
            sâ‚�â‚Š = sum(Aâ‚�â‚Š)
            #
            Aâ‚‚â‚Š = rand(rng, F64, hselectâ‚Š)
            sâ‚‚â‚Š = sum(Aâ‚‚â‚Š)
            @. Aâ‚‚â‚Š = Aâ‚‚â‚Š / sâ‚‚â‚Š * sâ‚�â‚Š
            #
            @. SE.A[selectedâ‚Š] = Aâ‚‚â‚Š

            # For negative-weight poles
            Aâ‚�â‚‹ = SE.A[selectedâ‚‹]
            sâ‚�â‚‹ = sum(Aâ‚�â‚‹)
            #
            Aâ‚‚â‚‹ = rand(rng, F64, hselectâ‚‹)
            sâ‚‚â‚‹ = sum(Aâ‚‚â‚‹)
            @. Aâ‚‚â‚‹ = Aâ‚‚â‚‹ / sâ‚‚â‚‹ * sâ‚�â‚‹
            #
            @. SE.A[selectedâ‚‹] = Aâ‚‚â‚‹
        end
    # For diagonal elements
    else
        # How many poles should be changed
        if npole â‰¤ 5
            if 4 â‰¤ npole â‰¤ 5
                nselect = 2
            else
                nselect = 1
            end
        else
            nselect = npole Ã· 5
        end
        @assert nselect â‰¤ npole

        # Which poles should be changed
        selected = rand(rng, 1:npole, nselect)
        unique!(selected)
        nselect = length(selected)

        # Change poles' positions
        if rand(rng) < 0.5
            P = rand(rng, allow, nselect)
            @. SE.P[selected] = P
        # Change poles' amplitudes
        else
            Aâ‚� = SE.A[selected]
            sâ‚� = sum(Aâ‚�)
            #
            Aâ‚‚ = rand(rng, F64, nselect)
            sâ‚‚ = sum(Aâ‚‚)
            @. Aâ‚‚ = Aâ‚‚ / sâ‚‚ * sâ‚�
            #
            @. SE.A[selected] = Aâ‚‚
        end
    end
end

"""
    reset_context(t::I64, SE::StochPXElement, SC::StochPXContext)

Recalculate imaginary frequency Green's function Gáµ§ and goodness-of-fit
functional Ï‡Â² by new Monte Carlo field configurations for the `t`-th
attempts. They must be consistent with each other.

Some key variables in `SC` are also updated as well. Perhaps we should
develop a smart algorhtm to update Î˜ here.

### Arguments
* t -> Counter for attempts.
* SE -> A StochPXElement struct.
* SC -> A StochPXContext struct.

### Returns
N/A
"""
function reset_context(t::I64, SE::StochPXElement, SC::StochPXContext)
    SC.Î˜ = get_x("theta")
    SC.Gáµ§ = calc_green(SE.P, SE.A, SE.ğ�•Š, SC.Î›)
    SC.Ï‡Â² = calc_chi2(SC.Gáµ§, SC.Gáµ¥)
    SC.Ï‡Â²áµ¥[t] = 1e10
end

"""
    calc_fmesh(S::StochPXSolver)

Try to calculate very fine (dense) mesh in [wmin, wmax], which is used
internally to represent the possible positions of poles. 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 generate
a linear mesh automatically.

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

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

See also: [`LinearMesh`](@ref), [`DynamicMesh`](@ref).
"""
function calc_fmesh(S::StochPXSolver)
    wmin = get_b("wmin")
    wmax = get_b("wmax")
    nfine = get_x("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

#=
*Remarks* :

Here, we would like to discuss the kernel matrix and sum-rule in the
stochastic pole expansion algorithm.

**For Fermionic Green's Functions**

The pole expansion of fermionic Green's function reads:

```math
\begin{equation}
G(i\omega_n) = \sum^{N_p}_{i = 1} \frac{A_i}{i\omega_n - P_i}.
\end{equation}
```

The kernel matrix reads:

```math
\begin{equation}
Î›(\omega_n, \omega) = \frac{1}{i\omega_n - \omega}.
\end{equation}
```

The sum-rule reads:

```math
\begin{equation}
\sum^{N_p}_{i=1} A_i = 1.
\end{equation}
```

---

**For Bosonic Green's Functions**

Supposed that ``\chi_0 \equiv -G(i\omega_n = 0)``, which must be a
positive real number, then the pole expansion of bosonic Green's
function reads:

```math
\begin{equation}
G(i\omega_n) = \sum^{N_p}_{i = 1} \frac{A_i}{i\omega_n - P_i}
= \sum^{N_p}_{i=1}
\frac{\chi_0 P_i}{i\omega_n - P_i}
\frac{A_i}{\chi_0 P_i}.
\end{equation}
```

The kernel matrix reads:

```math
\begin{equation}
Î›(\omega_n, \omega) = \frac{\chi_0 \omega}{i\omega_n - \omega}.
\end{equation}
```

Be careful, when ``\omega_n = 0``, ``\Lambda(0,\omega) = -\chi_0``.

The sum-rule reads:

```math
\begin{equation}
\sum^{N_p}_{i=1} \tilde{A}_i = 1,
\end{equation}
```

where

```math
\begin{equation}
\tilde{A}_i = \frac{A_i}{\chi_0 P_i}.
\end{equation}
```

---

**For Correlator of Hermitian Operator**

The pole expansion for correlator of Hermitian operator reads:

```math
\begin{equation}
G(i\omega_n) = \sum^{N_p}_{i = 1} A_i
\left(\frac{1}{i\omega_n - P_i} - \frac{1}{i\omega_n + P_i}\right)
= \sum^{N_p}_{i =1}
\frac{-\chi_0 P^2_i}{\omega^2_n + P^2_i} \frac{2A_i}{\chi_0 P_i}.
\end{equation}
```

Here, ``\chi_0 \equiv -G(i\omega_n = 0)``, which is a positive real
number. And ``\forall i, A_i \ge 0, P_i \ge 0``.

The kernel matrix reads:

```math
\begin{equation}
Î›(\omega_n, \omega) = \frac{-\chi_0 \omega^2}{\omega^2_n + \omega^2}.
\end{equation}
```

Be careful, when ``\omega_n = 0``, ``\Lambda(0,\omega) = -\chi_0``.

The sum-rule reads:

```math
\begin{equation}
\sum^{N_p}_{i=1} \tilde{A}_i = 1,
\end{equation}
```

where

```math
\begin{equation}
\tilde{A}_i = \frac{2A_i}{\chi_0 P_i}.
\end{equation}
```

---

=#

"""
    calc_lambda(
        grid::AbstractGrid,
        fmesh::AbstractMesh,
        Gáµ¥::Vector{F64}
    )

Precompute the kernel matrix Î› (Î› â‰¡ 1 / (iÏ‰â‚™ - Ïµ)). It is the essential
driver function. Note that Î› depends on the kernel's type (`ktype`).

### Arguments
* grid  -> Imaginary axis grid for input data.
* fmesh -> Very dense mesh in [wmin, wmax].
* Gáµ¥ -> Preprocessed input correlator.

### Returns
* Î› -> The kernel matrix, a 2D array.
"""
function calc_lambda(
    grid::AbstractGrid,
    fmesh::AbstractMesh,
    Gáµ¥::Vector{F64}
    )
    ktype = get_b("ktype")
    Ï‡â‚€ = -Gáµ¥[1]

    @cswitch ktype begin
        #
        @case "fermi"
            Î› = calc_lambda(grid, fmesh)
            break
        #
        @case "boson"
            Î› = calc_lambda(grid, fmesh, Ï‡â‚€, false)
            break
        #
        @case "bsymm"
            Î› = calc_lambda(grid, fmesh, Ï‡â‚€, true)
            break
        #
    end

    return Î›
end

"""
    calc_lambda(grid::AbstractGrid, fmesh::AbstractMesh)

Precompute the kernel matrix Î› (Î› â‰¡ 1 / (iÏ‰â‚™ - Ïµ)). It is a service
function and is for the fermionic systems.

### Arguments
* grid  -> Imaginary axis grid for input data.
* fmesh -> Very dense mesh in [wmin, wmax].

### Returns
* Î› -> The kernel matrix, a 2D array.
"""
function calc_lambda(grid::AbstractGrid, fmesh::AbstractMesh)
    ngrid = get_b("ngrid")
    nfine = get_x("nfine")

    _Î› = zeros(C64, ngrid, nfine)
    #
    for i in eachindex(grid)
        iÏ‰â‚™ = im * grid[i]
        for j in eachindex(fmesh)
            _Î›[i,j] = 1.0 / (iÏ‰â‚™ - fmesh[j])
        end
    end
    #
    Î› = vcat(real(_Î›), imag(_Î›))

    return Î›
end

"""
    calc_lambda(
        grid::AbstractGrid,
        fmesh::AbstractMesh,
        Ï‡â‚€::F64,
        bsymm::Bool
    )

Precompute the kernel matrix Î›. Here, `Ï‡â‚€` is actually -G(iÏ‰â‚™ = 0). And
the argument `bsymm` is used to distinguish two different bosonic kernels.
If `bsymm` is false, it means that the kernel is `boson`. If `bsymm` is
true, the kernel is `bsymm`. This function is for the bosonic systems.

### Arguments
* grid  -> Imaginary axis grid for input data.
* fmesh -> Very dense mesh in [wmin, wmax].
* Ï‡â‚€    -> -G(0).
* bsymm -> Type of bosonic kernel.

### Returns
* Î› -> The kernel matrix, a 2D array.
"""
function calc_lambda(
    grid::AbstractGrid,
    fmesh::AbstractMesh,
    Ï‡â‚€::F64,
    bsymm::Bool
    )
    ngrid = get_b("ngrid")
    nfine = get_x("nfine")

    # For standard bosonic kernel matrix
    if bsymm == false

        _Î› = zeros(C64, ngrid, nfine)
        #
        for i in eachindex(grid)
            iÏ‰â‚™ = im * grid[i]
            for j in eachindex(fmesh)
                _Î›[i,j] = Ï‡â‚€ * fmesh[j] / (iÏ‰â‚™ - fmesh[j])
            end
        end
        #
        # Special treatment for iÏ‰â‚™ = 0
        for j in eachindex(fmesh)
            _Î›[1,j] = -Ï‡â‚€
        end

        Î› = vcat(real(_Î›), imag(_Î›))

    # For symmetric bosonic kernel matrix
    else

        _Î› = zeros(F64, ngrid, nfine)
        #
        for i in eachindex(grid)
            Ï‰â‚™ = grid[i]
            for j in eachindex(fmesh)
                _Î›[i,j] = -Ï‡â‚€ * (fmesh[j] ^ 2.0) / (Ï‰â‚™ ^ 2.0 + fmesh[j] ^ 2.0)
            end
        end
        #
        # Special treatment for Ï‰â‚™ = 0
        for j in eachindex(fmesh)
            _Î›[1,j] = -Ï‡â‚€
        end
        #
        Î› = copy(_Î›)

    end

    return Î›
end

"""
    calc_green(t::I64, SC::StochPXContext, real_axis::Bool)

Reconstruct Green's function at imaginary axis or real axis by using the
pole expansion. It is a driver function. If `real_axis = true`, it will
returns G(Ï‰), or else G(iÏ‰â‚™).

### Arguments
* t -> Index of the current attempt.
* SC -> A StochPXContext struct.
* real_axis -> Working at real axis (true) or imaginary axis (false)?

### Returns
* G -> Reconstructed Green's function, G(Ï‰) or G(iÏ‰â‚™).
"""
function calc_green(t::I64, SC::StochPXContext, real_axis::Bool)
    ktype = get_b("ktype")
    ntry = get_x("ntry")
    @assert t â‰¤ ntry

    # Calculate G(iÏ‰â‚™)
    if real_axis == false
        return calc_green(SC.Páµ¥[t], SC.Aáµ¥[t], SC.ğ�•Šáµ¥[t], SC.Î›)
    end

    # Calculate G(Ï‰). Now we don't need SC.Î›.
    Ï‡â‚€ = -SC.Gáµ¥[1]
    @cswitch ktype begin
        @case "fermi"
            G = calc_green(SC.Páµ¥[t], SC.Aáµ¥[t], SC.ğ�•Šáµ¥[t], SC.mesh, SC.fmesh)
            break

        @case "boson"
            G = calc_green(SC.Páµ¥[t], SC.Aáµ¥[t], SC.ğ�•Šáµ¥[t], SC.mesh, SC.fmesh, Ï‡â‚€, false)
            break

        @case "bsymm"
            G = calc_green(SC.Páµ¥[t], SC.Aáµ¥[t], SC.ğ�•Šáµ¥[t], SC.mesh, SC.fmesh, Ï‡â‚€, true)
            break
    end

    return G
end

"""
    calc_green(
        P::Vector{I64},
        A::Vector{F64},
        ğ�•Š::Vector{F64},
        Î›::Array{F64,2}
    )

Reconstruct Green's function at imaginary axis by the pole expansion.

### Arguments
* P -> Positions of poles.
* A -> Amplitudes of poles.
* ğ�•Š -> Signs of poles.
* Î› -> Kernel matrix Î›.

### Returns
* G -> Reconstructed Green's function, G(iÏ‰â‚™).
"""
function calc_green(
    P::Vector{I64},
    A::Vector{F64},
    ğ�•Š::Vector{F64},
    Î›::Array{F64,2}
    )
    # Note that here `ngrid` is equal to 2 Ã— ngrid sometimes.
    ngrid, _ = size(Î›)

    G = zeros(F64, ngrid)
    for i = 1:ngrid
        G[i] = dot(A .* ğ�•Š, Î›[i,P])
    end

    return G
end

"""
    calc_green(
        P::Vector{I64},
        A::Vector{F64},
        ğ�•Š::Vector{F64},
        mesh::AbstractMesh,
        fmesh::AbstractMesh
    )

Reconstruct Green's function at real axis by the pole expansion. It is
for the fermionic systems only.

### Arguments
* P     -> Positions of poles.
* A     -> Amplitudes of poles.
* ğ�•Š     -> Signs of poles.
* mesh  -> Real frequency mesh for spectral functions.
* fmesh -> Very dense real frequency mesh for poles.

### Returns
* G -> Retarded Green's function, G(Ï‰).
"""
function calc_green(
    P::Vector{I64},
    A::Vector{F64},
    ğ�•Š::Vector{F64},
    mesh::AbstractMesh,
    fmesh::AbstractMesh
    )
    Î· = get_x("eta")
    nmesh = length(mesh)

    iÏ‰â‚™ = mesh.mesh .+ im * Î·
    G = zeros(C64, nmesh)
    for i in eachindex(mesh)
        G[i] = sum( @. (A * ğ�•Š) / (iÏ‰â‚™[i] - fmesh.mesh[P]) )
    end

    return G
end

"""
    calc_green(
        P::Vector{I64},
        A::Vector{F64},
        ğ�•Š::Vector{F64},
        mesh::AbstractMesh,
        fmesh::AbstractMesh,
        Ï‡â‚€::F64,
        bsymm::Bool
    )

Reconstruct Green's function at real axis by the pole expansion. Here,
`Ï‡â‚€` is actually -G(iÏ‰â‚™ = 0). And the argument `bsymm` is used to
distinguish two different bosonic kernels. If `bsymm` is false, it means
that the kernel is `boson`. If `bsymm` is true, the kernel is `bsymm`.
It is for the bosonic systems only.

### Arguments
* P     -> Positions of poles.
* A     -> Amplitudes of poles.
* ğ�•Š     -> Signs of poles.
* mesh  -> Real frequency mesh for spectral functions.
* fmesh -> Very dense real frequency mesh for poles.
* Ï‡â‚€    -> -G(0).
* bsymm -> Type of bosonic kernel.

### Returns
* G -> Retarded Green's function, G(Ï‰).
"""
function calc_green(
    P::Vector{I64},
    A::Vector{F64},
    ğ�•Š::Vector{F64},
    mesh::AbstractMesh,
    fmesh::AbstractMesh,
    Ï‡â‚€::F64,
    bsymm::Bool
    )
    Î· = get_x("eta")
    nmesh = length(mesh)

    iÏ‰â‚™ = mesh.mesh .+ im * Î·
    G = zeros(C64, nmesh)
    if bsymm == false
        _A = A .* ğ�•Š .* Ï‡â‚€ .* fmesh.mesh[P]
        for i in eachindex(mesh)
            G[i] = sum( @. _A / (iÏ‰â‚™[i] - fmesh.mesh[P]) )
        end
    #
    else
        _A = A .* ğ�•Š .* Ï‡â‚€ .* fmesh.mesh[P] .* 0.5
        for i in eachindex(mesh)
            Gâ‚Š = sum( @. _A / (iÏ‰â‚™[i] - fmesh.mesh[P]) )
            Gâ‚‹ = sum( @. _A / (iÏ‰â‚™[i] + fmesh.mesh[P]) )
            G[i] = Gâ‚Š - Gâ‚‹
        end
    #
    end

    return G
end

"""
    calc_chi2(Gâ‚™::Vector{F64}, Gáµ¥::Vector{F64})

Try to calculate the goodness-of-fit function (i.e, Ï‡Â²), which measures
the distance between input and regenerated correlators.

### Arguments
* Gâ‚™ -> Reconstructed Green's function.
* Gáµ¥ -> Original Green's function.

### Returns
* chi2 -> Goodness-of-fit functional, Ï‡Â².

See also: [`calc_green`](@ref).
"""
function calc_chi2(Gâ‚™::Vector{F64}, Gáµ¥::Vector{F64})
    Î”G = Gâ‚™ - Gáµ¥
    return dot(Î”G, Î”G)
end

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

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

If `offdiag = true`, the regions defined by `exclude` are for A(Ï‰) < 0,
the corresponding signs of poles are negative. While the other regions
that are not included in `exclude` are for A(Ï‰) > 0, and the corresponding
signs of poles are positive.

### Arguments
* S     -> A StochPXSolver struct.
* fmesh -> Very dense mesh for the poles.

### Returns
* allow -> Allowable indices.

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

    allow = I64[]
    unallow = 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)
        else
            push!(unallow, -i)
        end
    end

    # If it is offdiagonal, then the spectral function can be negative.
    # Now `allow` is for A(Ï‰) > 0 and `unallow` is for A(Ï‰) < 0. We have
    # to distinguish them.
    if offdiag
        append!(allow, unallow)
    end

    return allow
end

"""
    try_move_s(
        t::I64,
        MC::StochPXMC,
        SE::StochPXElement,
        SC::StochPXContext
    )

Change the position of one randomly selected pole.

### Arguments
* t -> Counter of attempts.
* MC -> A StochPXMC struct.
* SE -> A StochPXElement struct.
* SC -> A StochPXContext struct.

### Returns
N/A

See also: [`try_move_p`](@ref).
"""
function try_move_s(
    t::I64,
    MC::StochPXMC,
    SE::StochPXElement,
    SC::StochPXContext
    )
    # Get parameters
    ngrid = length(SC.Gáµ¥) # get_b("ngrid")
    nfine = get_x("nfine")
    npole = get_x("npole")
    move_window = nfine Ã· 100

    # It is used to save the change of Green's function
    Î´G = zeros(F64, ngrid)
    Gâ‚™ = zeros(F64, ngrid)

    # Try to go through each pole
    for _ = 1:npole

        # Select one pole randomly
        s = rand(MC.rng, 1:npole)

        # Try to change position of the s pole
        Aâ‚› = SE.A[s]
        ğ�•Šâ‚› = SE.ğ�•Š[s]
        #
        Î´P = rand(MC.rng, 1:move_window)
        #
        Pâ‚� = SE.P[s]
        Pâ‚‚ = Pâ‚�
        if rand(MC.rng) > 0.5
            Pâ‚‚ = Pâ‚� + Î´P
        else
            Pâ‚‚ = Pâ‚� - Î´P
        end
        #
        if ğ�•Šâ‚› > 0.0
            !(+Pâ‚‚ in SC.allow) && continue
        else
            !(-Pâ‚‚ in SC.allow) && continue
        end

        # Calculate change of Green's function
        Î›â‚� = view(SC.Î›, :, Pâ‚�)
        Î›â‚‚ = view(SC.Î›, :, Pâ‚‚)
        @. Î´G = ğ�•Šâ‚› * Aâ‚› * (Î›â‚‚ - Î›â‚�)

        # Calculate new Green's function and goodness-of-fit function
        @. Gâ‚™ = Î´G + SC.Gáµ§
        Ï‡Â² = calc_chi2(Gâ‚™, SC.Gáµ¥)
        Î´Ï‡Â² = Ï‡Â² - SC.Ï‡Â²

        # Simulated annealing algorithm
        MC.Stry = MC.Stry + 1
        if Î´Ï‡Â² < 0 || min(1.0, exp(-Î´Ï‡Â² * SC.Î˜)) > rand(MC.rng)
            # Update Monte Carlo configuration
            SE.P[s] = Pâ‚‚

            # Update reconstructed Green's function
            @. SC.Gáµ§ = Gâ‚™

            # Update goodness-of-fit function
            SC.Ï‡Â² = Ï‡Â²

            # Update Monte Carlo counter
            MC.Sacc = MC.Sacc + 1

            # Save optimal solution
            if SC.Ï‡Â² < SC.Ï‡Â²áµ¥[t]
                measure(t, SE, SC)
            end
        end

    end
end

"""
    try_move_p(
        t::I64,
        MC::StochPXMC,
        SE::StochPXElement,
        SC::StochPXContext
    )

Change the positions of two randomly selected poles.

### Arguments
* t -> Counter of attempts.
* MC -> A StochPXMC struct.
* SE -> A StochPXElement struct.
* SC -> A StochPXContext struct.

### Returns
N/A

See also: [`try_move_s`](@ref).
"""
function try_move_p(
    t::I64,
    MC::StochPXMC,
    SE::StochPXElement,
    SC::StochPXContext
    )
    # Get parameters
    ngrid = length(SC.Gáµ¥) # get_b("ngrid")
    npole = get_x("npole")

    # Sanity check
    if npole == 1
        return
    end

    # It is used to save the change of Green's function
    Î´G = zeros(F64, ngrid)
    Gâ‚™ = zeros(F64, ngrid)

    # Try to go through each pole
    for _ = 1:npole

        # Select two poles randomly
        # The two poles should not be the same.
        sâ‚� = 1
        sâ‚‚ = 1
        while sâ‚� == sâ‚‚
            sâ‚� = rand(MC.rng, 1:npole)
            sâ‚‚ = rand(MC.rng, 1:npole)
        end

        # Try to change position of the sâ‚� pole
        Aâ‚� = SE.A[sâ‚�]
        ğ�•Šâ‚� = SE.ğ�•Š[sâ‚�]
        Pâ‚� = SE.P[sâ‚�]
        Pâ‚ƒ = Pâ‚�
        while Pâ‚ƒ == Pâ‚� || sign(Pâ‚ƒ) != sign(ğ�•Šâ‚�)
            Pâ‚ƒ = rand(MC.rng, SC.allow)
        end
        Pâ‚ƒ = abs(Pâ‚ƒ)
        #
        # Try to change position of the sâ‚‚ pole
        Aâ‚‚ = SE.A[sâ‚‚]
        ğ�•Šâ‚‚ = SE.ğ�•Š[sâ‚‚]
        Pâ‚‚ = SE.P[sâ‚‚]
        Pâ‚„ = Pâ‚‚
        while Pâ‚„ == Pâ‚‚ || sign(Pâ‚„) != sign(ğ�•Šâ‚‚)
            Pâ‚„ = rand(MC.rng, SC.allow)
        end
        Pâ‚„ = abs(Pâ‚„)

        # Calculate change of Green's function
        Î›â‚� = view(SC.Î›, :, Pâ‚�)
        Î›â‚‚ = view(SC.Î›, :, Pâ‚‚)
        Î›â‚ƒ = view(SC.Î›, :, Pâ‚ƒ)
        Î›â‚„ = view(SC.Î›, :, Pâ‚„)
        @. Î´G = ğ�•Šâ‚� * Aâ‚� * (Î›â‚ƒ - Î›â‚�) + ğ�•Šâ‚‚ * Aâ‚‚ * (Î›â‚„ - Î›â‚‚)

        # Calculate new Green's function and goodness-of-fit function
        @. Gâ‚™ = Î´G + SC.Gáµ§
        Ï‡Â² = calc_chi2(Gâ‚™, SC.Gáµ¥)
        Î´Ï‡Â² = Ï‡Â² - SC.Ï‡Â²

        # Simulated annealing algorithm
        MC.Ptry = MC.Ptry + 1
        if Î´Ï‡Â² < 0 || min(1.0, exp(-Î´Ï‡Â² * SC.Î˜)) > rand(MC.rng)
            # Update Monte Carlo configuration
            SE.P[sâ‚�] = Pâ‚ƒ
            SE.P[sâ‚‚] = Pâ‚„

            # Update reconstructed Green's function
            @. SC.Gáµ§ = Gâ‚™

            # Update goodness-of-fit function
            SC.Ï‡Â² = Ï‡Â²

            # Update Monte Carlo counter
            MC.Pacc = MC.Pacc + 1

            # Save optimal solution
            if SC.Ï‡Â² < SC.Ï‡Â²áµ¥[t]
                measure(t, SE, SC)
            end
        end

    end
end

"""
    try_move_a(
        t::I64,
        MC::StochPXMC,
        SE::StochPXElement,
        SC::StochPXContext
    )

Change the amplitudes of two randomly selected poles.

### Arguments
* t -> Counter of attempts.
* MC -> A StochPXMC struct.
* SE -> A StochPXElement struct.
* SC -> A StochPXContext struct.

### Returns
N/A

See also: [`try_move_x`](@ref).
"""
function try_move_a(
    t::I64,
    MC::StochPXMC,
    SE::StochPXElement,
    SC::StochPXContext
    )
    # Get parameters
    ngrid = length(SC.Gáµ¥) # get_b("ngrid")
    npole = get_x("npole")

    # Sanity check
    if npole == 1
        return
    end

    # It is used to save the change of Green's function
    Î´G = zeros(F64, ngrid)
    Gâ‚™ = zeros(F64, ngrid)

    # Try to go through each pole
    for _ = 1:npole

        # Select two poles randomly
        # The two poles should not be the same.
        sâ‚� = 1
        sâ‚‚ = 1
        while sâ‚� == sâ‚‚
            sâ‚� = rand(MC.rng, 1:npole)
            sâ‚‚ = rand(MC.rng, 1:npole)
        end

        # Try to change amplitudes of the two poles, but their sum is kept.
        Pâ‚� = SE.P[sâ‚�]
        Pâ‚‚ = SE.P[sâ‚‚]
        Aâ‚� = SE.A[sâ‚�]
        Aâ‚‚ = SE.A[sâ‚‚]
        Aâ‚ƒ = 0.0
        Aâ‚„ = 0.0
        ğ�•Šâ‚� = SE.ğ�•Š[sâ‚�]
        ğ�•Šâ‚‚ = SE.ğ�•Š[sâ‚‚]

        if ğ�•Šâ‚� == ğ�•Šâ‚‚
            while true
                Î´A = rand(MC.rng) * (Aâ‚� + Aâ‚‚) - Aâ‚�
                Aâ‚ƒ = Aâ‚� + Î´A
                Aâ‚„ = Aâ‚‚ - Î´A

                if 1.0 > Aâ‚ƒ > 0.0 && 1.0 > Aâ‚„ > 0.0
                    break
                end
            end
        else
            while true
                _Î´A = rand(MC.rng) * (Aâ‚� + Aâ‚‚) - Aâ‚�
                Î´A = rand(MC.rng) > 0.5 ? _Î´A * (+1.0) : _Î´A * (-1.0)
                Aâ‚ƒ = (ğ�•Šâ‚� * Aâ‚� + Î´A) / ğ�•Šâ‚�
                Aâ‚„ = (ğ�•Šâ‚‚ * Aâ‚‚ - Î´A) / ğ�•Šâ‚‚

                if 1.0 > Aâ‚ƒ > 0.0 && 1.0 > Aâ‚„ > 0.0
                    break
                end
            end
        end

        # Calculate change of Green's function
        Î›â‚� = view(SC.Î›, :, Pâ‚�)
        Î›â‚‚ = view(SC.Î›, :, Pâ‚‚)
        @. Î´G = ğ�•Šâ‚� * (Aâ‚ƒ - Aâ‚�) * Î›â‚� + ğ�•Šâ‚‚ * (Aâ‚„ - Aâ‚‚) * Î›â‚‚

        # Calculate new Green's function and goodness-of-fit function
        @. Gâ‚™ = Î´G + SC.Gáµ§
        Ï‡Â² = calc_chi2(Gâ‚™, SC.Gáµ¥)
        Î´Ï‡Â² = Ï‡Â² - SC.Ï‡Â²

        # Simulated annealing algorithm
        MC.Atry = MC.Atry + 1
        if Î´Ï‡Â² < 0 || min(1.0, exp(-Î´Ï‡Â² * SC.Î˜)) > rand(MC.rng)
            # Update Monte Carlo configuration
            SE.A[sâ‚�] = Aâ‚ƒ
            SE.A[sâ‚‚] = Aâ‚„

            # Update reconstructed Green's function
            @. SC.Gáµ§ = Gâ‚™

            # Update goodness-of-fit function
            SC.Ï‡Â² = Ï‡Â²

            # Update Monte Carlo counter
            MC.Aacc = MC.Aacc + 1

            # Save optimal solution
            if SC.Ï‡Â² < SC.Ï‡Â²áµ¥[t]
                measure(t, SE, SC)
            end
        end

    end
end

"""
    try_move_x(
        t::I64,
        MC::StochPXMC,
        SE::StochPXElement,
        SC::StochPXContext
    )

Exchange the amplitudes of two randomly selected poles.

### Arguments
* t -> Counter of attempts.
* MC -> A StochPXMC struct.
* SE -> A StochPXElement struct.
* SC -> A StochPXContext struct.

### Returns
N/A

See also: [`try_move_a`](@ref).
"""
function try_move_x(t::I64, MC::StochPXMC, SE::StochPXElement, SC::StochPXContext)
    # Get parameters
    ngrid = length(SC.Gáµ¥) # get_b("ngrid")
    offdiag = get_b("offdiag")
    npole = get_x("npole")

    # Sanity check
    if offdiag
        if npole â‰¤ 3
            return
        end
    else
        if npole == 1
            return
        end
    end

    # It is used to save the change of Green's function
    Î´G = zeros(F64, ngrid)
    Gâ‚™ = zeros(F64, ngrid)

    # Try to go through each pole
    for _ = 1:npole

        # Select two poles randomly
        # The positions of the two poles are different,
        # but their signs should be the same.
        sâ‚� = 1
        sâ‚‚ = 1
        while ( sâ‚� == sâ‚‚ ) || ( SE.ğ�•Š[sâ‚�] != SE.ğ�•Š[sâ‚‚] )
            sâ‚� = rand(MC.rng, 1:npole)
            sâ‚‚ = rand(MC.rng, 1:npole)
        end

        # Try to swap amplitudes of the two poles, but their sum is kept.
        Pâ‚� = SE.P[sâ‚�]
        Pâ‚‚ = SE.P[sâ‚‚]
        Aâ‚� = SE.A[sâ‚�]
        Aâ‚‚ = SE.A[sâ‚‚]
        Aâ‚ƒ = Aâ‚‚
        Aâ‚„ = Aâ‚�
        ğ�•Šâ‚� = SE.ğ�•Š[sâ‚�]
        ğ�•Šâ‚‚ = SE.ğ�•Š[sâ‚‚]

        # Calculate change of Green's function
        Î›â‚� = view(SC.Î›, :, Pâ‚�)
        Î›â‚‚ = view(SC.Î›, :, Pâ‚‚)
        @. Î´G = ğ�•Šâ‚� * (Aâ‚ƒ - Aâ‚�) * Î›â‚� + ğ�•Šâ‚‚ * (Aâ‚„ - Aâ‚‚) * Î›â‚‚

        # Calculate new Green's function and goodness-of-fit function
        @. Gâ‚™ = Î´G + SC.Gáµ§
        Ï‡Â² = calc_chi2(Gâ‚™, SC.Gáµ¥)
        Î´Ï‡Â² = Ï‡Â² - SC.Ï‡Â²

        # Simulated annealing algorithm
        MC.Xtry = MC.Xtry + 1
        if Î´Ï‡Â² < 0 || min(1.0, exp(-Î´Ï‡Â² * SC.Î˜)) > rand(MC.rng)
            # Update Monte Carlo configuration
            SE.A[sâ‚�] = Aâ‚ƒ
            SE.A[sâ‚‚] = Aâ‚„

            # Update reconstructed Green's function
            @. SC.Gáµ§ = Gâ‚™

            # Update goodness-of-fit function
            SC.Ï‡Â² = Ï‡Â²

            # Update Monte Carlo counter
            MC.Xacc = MC.Xacc + 1

            # Save optimal solution
            if SC.Ï‡Â² < SC.Ï‡Â²áµ¥[t]
                measure(t, SE, SC)
            end
        end

    end
end
