Solvers

AbstractSolver

An abstract type representing the solver for analytic continuation problem. It is used to build the internal type system. All the other solvers are its subtypes.

It has the following subtypes:

• MaxEntSolver
• BarRatSolver
• NevanACSolver
• StochACSolver
• StochSKSolver
• StochOMSolver
• StochPXSolver

AbstractMC

An abstract type representing the Monte Carlo engine. It is used to build the internal type system.

MaxEntSolver

It represents the analytic continuation solver that implements the maximum entropy method.

This solver is highly recommended.

Features

• Both off-diagonal and diagonal Green's functions are supported.
• Both fermionic and bosonic correlators are supported.
• The spectral weights can be negative.
• Good robustness with respect to noisy correlators.
• Numerically stable.
• Fast.

Known Limitations

• It always tends to generate broad and smooth spectra.
• The fine features in the spectra could be smeared out.
• Sometimes it would generate keen-edged and unphysical peaks.
• When the spectra are sharp or discrete, this solver usually fails.

Implementation

See src/maxent.jl.

Tests

See test/basic.

MaxEntContext

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

Members

• Gᵥ -> Input data for correlator.
• σ² -> Actually 1.0 / σ².
• grid -> Grid for input data.
• mesh -> Mesh for output spectrum.
• model -> Default model function.
• kernel -> Default kernel function.
• Vₛ -> Matrix from singular value decomposition.
• W₂ -> Precomputed array.
• W₃ -> Precomputed array.
• Bₘ -> Precomputed array.

solve(S::MaxEntSolver, rd::RawData)

Solve the analytic continuation problem by the maximum entropy method. It is the driver for the MaxEnt 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 StochAC, StochSK, and StochOM solvers.

It seems that the MaxEnt solver is hard to create δ-like spectra.

Arguments

• S -> A MaxEntSolver 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(ω).

init(S::MaxEntSolver, rd::RawData)

Initialize the MaxEnt solver and return a MaxEntContext struct.

Arguments

• S -> A MaxEntSolver struct.
• rd -> A RawData struct, containing raw data for input correlator.

Returns

• mec -> A MaxEntContext struct.

run(mec::MaxEntContext)

Perform maximum entropy simulation with different algorithms. Now it supports the historic, classic, bryan, and chi2kink algorithms.

Arguments

• mec -> A MaxEntContext struct.

Returns

• svec -> A vector of dictionaries. It contains the intermediate solutions.
• sol -> Dictionary. It contains the final solution.

last(mec::MaxEntContext, svec::Vector, sol::Dict)

Postprocess the results generated during the maximum entropy simulations. Here sol is the final solution for the analytic continuation problem, while svec contains all the intermediate results (it is a vector of dictionary actually).

Arguments

• mec -> A MaxEntContext struct.
• svec -> See above explanations.
• sol -> See above explanations.

Returns

• G -> Retarded Green's function, G(ω).

historic(mec::MaxEntContext)

Apply the historic algorithm to solve the analytic continuation problem. It choose α in a way that χ² ≈ N.

For the historic algorithm, alpha is usually 10⁶, and ratio is 10.0. It is compatible with the Bayesian Reconstruction entropy.

Arguments

• mec -> A MaxEntContext struct.

Returns

• svec -> A vector of dictionaries. It contains the intermediate solutions.
• sol -> Dictionary. It contains the final solution.

See also: MaxEntContext.

classic(mec::MaxEntContext)

Apply the classic algorithm to solve the analytic continuation problem.

Classic algorithm uses Bayes statistics to approximately determine the most probable value of α. We always start at a large value of α, where the optimization yields basically the default model, therefore u_vec is only a few steps away from 0 (= default model). And then we gradually decrease α, step by step moving away from the default model towards data fitting. Using u_vec as start for the next (smaller) α brings a great speedup into this procedure.

For the classic algorithm, alpha is usually 10⁶, and ratio is 10.0. It is incompatible with the Bayesian Reconstruction entropy.

Arguments

• mec -> A MaxEntContext struct.

Returns

• svec -> A vector of dictionaries. It contains the intermediate solutions.
• sol -> Dictionary. It contains the final solution.

See also: MaxEntContext.

bryan(mec::MaxEntContext)

Apply the bryan algorithm to solve the analytic continuation problem.

Bryan's maxent calculates an average of spectral functions, weighted by their Bayesian probability.

For the bryan algorithm, alpha is usually 500, and ratio is 1.1. It is incompatible with the Bayesian Reconstruction entropy.

Arguments

• mec -> A MaxEntContext struct.

Returns

• svec -> A vector of dictionaries. It contains the intermediate solutions.
• sol -> Dictionary. It contains the final solution.

See also: MaxEntContext.

chi2kink(mec::MaxEntContext)

Apply the chi2kink algorithm to solve the analytic continuation problem.

We start with an optimization at a large value of α, where we should get only the default model. And then, α is decreased step-by-step, until the minimal value of α is reached. Then, we fit a function

ϕ(x; a, b, c, d) = a + b / [1 + exp(-d*(x-c))],

from which the optimal α is determined by

x_opt = c - fit_position / d,

and

alpha_opt = 10^x_opt.

For the chi2kink algorithm, alpha is usually 10⁹, ratio is 10.0, the number of alpha parameters is 12. It is compatible with the Bayesian Reconstruction entropy.

Arguments

• mec -> A MaxEntContext struct.

Returns

• svec -> A vector of dictionaries. It contains the intermediate solutions.
• sol -> Dictionary. It contains the final solution.

See also: MaxEntContext.

optimizer(
    mec::MaxEntContext,
    α::F64,
    us::Vector{F64},
    use_bayes::Bool
)

Optimization of maxent functional for a given value of α. Since a priori the best value of α is unknown, this function has to be called several times in order to find a good value.

α means a weight factor of the entropy. us is a vector in singular space. It is used as a starting value for the optimization. For the very first optimization, done at large α, we use zeros, which corresponds to the default model. Then we use the result of the previous optimization as a starting value. use_bayes determines whether to use the Bayesian inference parameters for α.

This function will return a dictionary object that holds the results of the optimization, e.g. spectral function, χ² deviation.

Arguments

• mec -> A MaxEntContext struct.
• α -> See above explanations.
• us -> See above explanations.
• use_bayes -> See above explanations.

Returns

• dict -> A dictionary, the solution to analytic continuation problem.

precompute(
    Gᵥ::Vector{F64},
    σ²::Vector{F64},
    am::AbstractMesh,
    D::Vector{F64},
    K::Matrix{F64}
)

Precompute some key coefficients. Here Gᵥ and σ² are input data, am is the mesh for spectrum, D is the default model, and K is the kernel function.

Arguments

• Gᵥ -> Input correlator.
• σ² -> Error bar for input correlator.
• am -> See above explanations.
• D -> See above explanations.
• K -> See above explanations.

Returns

• V -> An orthogonal matrix from singular value decomposition of kernel.
• W₂ -> The Wₘₗ matrix.
• W₃ -> The Wₘₗᵢ tensor.
• Bₘ -> The Bₘ vector.
• hess -> The Hessian matrix.

f_and_J(u::Vector{F64}, mec::MaxEntContext, α::F64)

This function evaluates the function whose root we want to find. Here u is a singular space vector that parametrizes the spectral function, and α is a (positive) weight factor of the entropy.

It returns f, value of the function whose zero we want to find, and J, jacobian at the current position.

Arguments

See above explanations.

Returns

See above explanations.

See also: f_and_J_od.

f_and_J_od(u::Vector{F64}, mec::MaxEntContext, α::F64)

This function evaluates the function whose root we want to find. Here u is a singular space vector that parametrizes the spectral function, and α is a (positive) weight factor of the entropy.

It returns f, value of the function whose zero we want to find, and J, jacobian at the current position.

This function is similar to f_and_J, but for offdiagonal elements.

Arguments

See above explanations.

Returns

See above explanations.

See also: f_and_J.

svd_to_real(mec::MaxEntContext, u::Vector{F64})

Go from singular value space to real space. It will transform the singular space vector u into real-frequency space (to get the spectral function) by A(ω) = D(ω) eⱽᵘ, where D(ω) is the default model, V is the matrix from the singular value decomposition. The argument u means a singular space vector that parametrizes the spectral function.

Arguments

See above explanations.

Returns

See above explanations.

See also: svd_to_real_od.

svd_to_real_od(mec::MaxEntContext, u::Vector{F64})

Go from singular value space to real space. It will transform the singular space vector u into real-frequency space in the case of an offdiagonal element. It will return the spectral function.

Arguments

• mec -> A MaxEntContext struct.
• u -> A singular space vector that parametrizes the spectral function.

Returns

See above explanations.

See also: svd_to_real.

calc_entropy(mec::MaxEntContext, A::Vector{F64}, u::Vector{F64})

It computes entropy for positive definite spectral function. Here the arguments A means spectral function and u means a singular space vector that parametrizes the spectral function.

Arguments

See above explanations.

Returns

• S -> Entropy.

See also: calc_entropy_od.

calc_entropy_od(mec::MaxEntContext, A::Vector{F64})

It compute positive-negative entropy for spectral function with norm 0. Here the argument A means spectral function.

Arguments

See above explanations.

Returns

• S -> Entropy.

See also: calc_entropy.

calc_bayes(
    mec::MaxEntContext,
    A::Vector{F64},
    S::F64,
    χ²::F64,
    α::F64
)

It calculates Bayesian convergence criterion (ng, tr, and conv) for classic maxent (maximum of probablility distribution) and then Bayesian a-posteriori probability (log_prob) for α after optimization of A.

Here, A is the spectral function, S the entropy, χ² the deviation, and α weight factor of the entropy.

Arguments

See above explanations.

Returns

• ng -> -2.0αS.
• tr -> Tr(Λ / (αI + Λ)).
• conv -> Ratio between ng and tr.
• prob -> Pr[α | \bar{G}].

See also: calc_bayes_od.

calc_bayes_od(
    mec::MaxEntContext,
    A::Vector{F64},
    S::F64,
    χ²::F64,
    α::F64
)

It calculates Bayesian convergence criterion (ng, tr, and conv) for classic maxent (maximum of probablility distribution) and then Bayesian a-posteriori probability (log_prob) for α after optimization of A.

Here, A is the spectral function, S the entropy, χ² the deviation, and α weight factor of the entropy.

It is just a offdiagonal version of calc_bayes().

Arguments

See above explanations.

Returns

• ng -> -2.0αS.
• tr -> Tr(Λ / (αI + Λ)).
• conv -> Ratio between ng and tr.
• prob -> Pr[α | \bar{G}].

See also: calc_bayes.

calc_chi2(mec::MaxEntContext, A::Vector{F64})

It computes the χ²-deviation of the spectral function A.

Arguments

• mec -> A MaxEntContext struct.
• A -> Spectral function.

Returns

• χ² -> Goodness-of-fit functional.

BarRatSolver

It represents the analytic continuation solver that implements the barycentric rational function approximation method.

This solver is highly recommended.

Features

• Both off-diagonal and diagonal Green's functions are supported.
• Both fermionic and bosonic correlators are supported.
• Both continuous and discrete spectra are supported.
• The spectral weights can be negative.
• Moderate robustness with respect to noisy correlators.
• Numerically stable.
• It can provide analytic expressions to approximate the correlators.
• Extremely fast.

Known Limitations

• The imaginary-time correlation correlators are not supported.
• Sometimes it would generate keen-edged and unphysical peaks.

Implementation

See src/rfa.jl.

Tests

See test/rfa.

BarRatContext

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

Members

• Gᵥ -> Input data for correlator.
• grid -> Grid for input data.
• mesh -> Mesh for output spectrum.
• 𝒫 -> Prony approximation for the input data.
• ℬ -> Barycentric rational function approximation for the input data.
• ℬP -> It means the positions of the poles.
• ℬA -> It means the weights / amplitudes of the poles.

BarycentricFunction

Mutable struct. Barycentric representation of a rational function.

Members

• nodes -> Nodes of the rational function, $z_i$.
• values -> Values of the rational function, $r(z_i)$.
• weights -> Weights of the rational function, $w_i$.
• wtimesf -> Weighted values of the rational function, $w_i f_i$.

PronyApproximation

Mutable struct. Prony approximation to a complex-valued Matsubara function.

Members

• 𝑁ₚ -> Number of nodes for Prony approximation.
• ωₚ -> Non-negative Matsubara frequency.
• 𝐺ₚ -> Complex values at ωₚ.
• Γₚ -> Nodes for Prony approximation, $γ_i$.
• Ωₚ -> Weights for Prony approximation, $w_i$.

solve(S::BarRatSolver, rd::RawData)

Solve the analytic continuation problem by the barycentric rational function method. This is the driver for the BarRat solver.

This solver suits Matsubara Green's functions. It supports both bosonic and fermionic systems, irrespective of diagonal or non-diagonal functions. It is extremely efficient. But sometimes the resulting spectral functions could violate the sum-rules.

Similar to the StochPX and NevanAC solvers, it always returns A(ω).

Arguments

• S -> A BarRatSolver 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(ω).

init(S::BarRatSolver, rd::RawData)

Initialize the BarRat solver and return a BarRatContext struct.

Arguments

• S -> A BarRatSolver struct.
• rd -> A RawData struct, containing raw data for input correlator.

Returns

• mec -> A BarRatContext struct.

run(brc::BarRatContext)

At first, it will try to construct a Prony approximation for the input Matsubara data. Then the Prony approximation is used to build smooth data set (data denoising). Finally, the barycentric rational function for this data set is constructed. The member ℬ of the BarRatContext struct (brc) should be updated in this function.

Arguments

• brc -> A BarRatContext struct.

Returns

N/A

last(brc::BarRatContext)

It will process and write the calculated results by the BarRat solver, including correlator at real axis, final spectral function, reproduced correlator. The information about Prony approximation and barycentric rational function approximation will be written as well.

Arguments

• brc -> A BarRatContext struct.

Returns

• Aout -> Spectral function, A(ω).
• Gout -> Retarded Green's function, G(ω).

aaa(z, y)

Adaptively compute a barycentric rational interpolant.

Arguments

• z::AbstractVector{<:Number} -> Interpolation nodes.
• y::AbstractVector{<:Number} -> Values at nodes.
• max_degree::Integer=150 -> Maximum numerator/denominator degree to use.
• float_type::Type=F64 -> Floating point type to use for the computation.
• tol::Real=1000*eps(float_type) -> Tolerance for stopping.
• lookahead::Integer=10 -> Number of iterations to determines stagnation.
• stats::Bool=false -> Return convergence statistics.

Returns

• r::BarycentricFunction -> The rational interpolant.
• stats::NamedTuple -> Convergence statistics, if keyword stats = true.

Examples

julia> z = 1im * range(-10, 10, 500);

julia> y = @. exp(z);

julia> r = aaa(z, y);

julia> bc_degree(r) # both numerator and denominator
12

julia> first(nodes(r), 4)
4-element Vector{ComplexF64}:
 0.0 - 6.272545090180361im
 0.0 + 9.43887775551102im
 0.0 - 1.1022044088176353im
 0.0 + 4.909819639278557im

julia> r(1im * π / 2)
-2.637151617496356e-15 + 1.0000000000000002im

poles!(brc::BarRatContext)

Convert the barycentric rational function approximation to the classic pole representation. Note that this feature is only suitable for the atype = "delta" case. In such case, the barycenteric algorithm can find the accurate positions for the poles via the bc_poles() function. But it seems that the weights for these poles are wrong. In this function, we just use the BFGS method to solve this optimization problem to get the correct weights for the poles. And then the positions and weights of these poles will be stored in brc, a BarRatContext struct.

Arguments

• brc -> A BarRatContext struct.

Returns

N/A

bc_nodes(r::BarycentricFunction)

Returns the nodes of the rational interpolant r as a vector. Actually, they are $z_i$.

bc_values(r::BarycentricFunction)

Returns the nodal values of the rational interpolant r as a vector. They are $r(z_i) = f_i$.

bc_weights(r::BarycentricFunction)

Returns the weights of the rational interpolant r as a vector. Actually, they are $w_i$.

bc_degree(r::BarycentricFunction)

Returns the degree of the numerator and denominator of the rational r.

bc_poles(r::BarycentricFunction)

Return the poles of the rational function r.

Arguments

• r -> A BarycentricFunction struct.

Returns

• pole -> List of poles.

prony_data(ω₁::Vector{F64}, 𝐺₁::Vector{C64})

Prepare essential data for the later Prony approximation. It will return the number of nodes 𝑁ₚ, frequency mesh ωₚ, and Green's function data 𝐺ₚ at this mesh.

Arguments

• ω₁ -> Non-negative Matsubara frequency (raw values).
• 𝐺₁ -> Complex values at ωₚ (raw values), 𝐺₁ = G(ω₁).

Returns

See above explanations.

prony_svd(𝑁ₚ::I64, 𝐺ₚ::Vector{C64})

Perform singular value decomposition for the matrix ℋ that is constructed from 𝐺ₚ. It will return the singular values S and orthogonal matrix V.

Arguments

• 𝑁ₚ -> Number of nodes.
• 𝐺ₚ -> Truncated Green's function data.

Returns

See above explanations.

See also: prony_data.

prony_idx(S::Vector{F64}, ε::F64)

The diagonal matrix (singular values) S is used to test whether the threshold ε is reasonable and figure out the index for extracting v from V.

Arguments

• S -> Singular values of ℋ.
• ε -> Threshold provided by the users.

Returns

• idx -> Index for S[idx] < ε.

See also: prony_v and prony_svd.

prony_idx(S::Vector{F64})

The diagonal matrix (singular values) S is used to figure out the index for extracting v from V. This function is try to evaluate the maximum index for the exponentially decaying region of S.

Arguments

• S -> Singular values of ℋ.

Returns

• idx -> Index for extracting v from V.

See also: prony_v.

prony_v(V::Adjoint{C64, Matrix{C64}}, idx::I64)

Extract suitable vector v from orthogonal matrix V according to the threshold ε.

Arguments

• V -> Orthogonal matrix from singular value decomposition of ℋ.
• idx -> Index for extracting v from V.

Returns

• v -> Vector v extracted from V according to idx.

See also: prony_svd.

prony_gamma(v::Vector{C64}, Λ::F64)

Try to calculate Γₚ. Actually, Γₚ are eigenvalues of a matrix constructed by v. Λ is a cutoff for Γₚ. Only those Γₚ that are smaller than Λ are kept.

Arguments

• v -> A vector extracted from V.
• Λ -> A cutoff for Γₚ.

Returns

• Γₚ -> Roots of a polynominal with coefficients given in v.

See also: prony_v.

prony_omega(𝐺ₚ::Vector{C64}, Γₚ::Vector{C64})

Try to calculate Ωₚ.

Arguments

• 𝐺ₚ -> Complex values at ωₚ.
• Γₚ -> Nodes for Prony approximation, $γ_i$.

Returns

• Ωₚ -> Weights for Prony approximation, $w_i$.

NevanACSolver

It represents the analytic continuation solver that implements the Nevanlinna analytical continuation (It doesn't support the analytic confinuations for matrix-valued Green's functions).

This solver is not recommended.

This solver is not well tested. If you would like to use the Nevanlinna analytical continuation method, perhaps the official Nevanlinna.jl is a better choice.

Features

• Only diagonal Green's functions are supported.
• Only fermionic correlators are supported.
• Both continuous and discrete spectra are supported.
• It is quite accurate if the input correlators are noise-free.

Known Limitations

• It needs Matsubara data.
• It is quite slow if Hardy basis optimization is activated.
• It doesn't support bosonic correlators directly.
• It doesn't support off-diagonal Green's functions.
• It doesn't suit the noisy quantum Monte Carlo data.
• It is numerically unstable when the input correlators are noisy.

Implementation

See src/nac.jl.

Tests

See test/nac.

NevanACContext

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

Members

• Gᵥ -> Input data for correlator.
• grid -> Grid for input data.
• mesh -> Mesh for output spectrum.
• Φ -> Φ vector in Schur algorithm.
• 𝒜 -> Coefficients matrix abcd in Schur algorithm.
• ℋ -> Hardy matrix for Hardy basis optimization.
• 𝑎𝑏 -> Coefficients matrix for expanding θ with Hardy basis.
• hmin -> Minimal value of the order of Hardy basis functions.
• hopt -> Optimal value of the order of Hardy basis functions.

solve(S::NevanACSolver, rd::RawData)

Solve the analytic continuation problem by the Nevanlinna analytical continuation method. It is the driver for the NevanAC solver.

This solver suits Matsubara Green's functions for fermionic systems. It can not be used directly to treat the bosonic correlators. It will return A(ω) all the time, similar to the StochPX and BarRat solvers.

This solver is numerically unstable. Sometimes it is hard to get converged solution, especially when the noise is medium.

Arguments

• S -> A NevanACSolver 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(ω).

init(S::NevanACSolver, rd::RawData)

Initialize the NevanAC solver and return a NevanACContext struct.

Arguments

• S -> A NevanACSolver struct.
• rd -> A RawData struct, containing raw data for input correlator.

Returns

• mec -> A NevanACContext struct.

run(nac::NevanACContext)

Perform Hardy basis optimization to smooth the spectrum. the members ℋ, 𝑎𝑏, hmin, and hopt of the NevanACContext struct (nac) should be updated in this function.

Arguments

• nac -> A NevanACContext struct.

Returns

N/A

last(nac::NevanACContext)

Postprocess the results generated during the Nevanlinna analytical continuation simulations.

Arguments

• nac -> A NevanACContext struct.

Returns

• Aout -> Spectral function, A(ω).
• Gout -> Retarded Green's function, G(ω).

precompute(
    grid::AbstractGrid,
    mesh::AbstractMesh,
    Gᵥ::Vector{APC}
)

Precompute some key quantities, such as Φ, 𝒜, ℋ, and 𝑎𝑏. Note that Φ and 𝒜 won't be changed any more. But ℋ and 𝑎𝑏 should be updated by the Hardy basis optimization to get a smooth spectrum. Here Gᵥ is input data, grid is the grid for input data, and mesh is the mesh for output spectrum.

Arguments

See above explanations.

Returns

• Φ -> Φ vector in Schur algorithm.
• 𝒜 -> Coefficients matrix abcd in Schur algorithm.
• ℋ -> Hardy matrix for Hardy basis optimization.
• 𝑎𝑏 -> Coefficients matrix for expanding θ with Hardy basis.

See also: NevanACContext.

calc_mobius(z::Vector{APC})

A direct Mobius transformation.

Arguments

• z -> Complex vector.

Returns

• val -> φ(z), Mobius transformation of z.

calc_inv_mobius(z::Vector{APC})

An inverse Mobius transformation.

Arguments

• z -> Complex vector.

Returns

• val -> φ⁻¹(z), inverse Mobius transformation of z.

calc_pick(k::I64, ℎ::Vector{APC}, λ::Vector{APC})

Try to calculate the Pick matrix, anc check whether it is a positive semidefinite matrix. See Eq. (5) in Fei's NAC paper.

Arguments

• k -> Size of the Pick matrix.
• ℎ -> Vector ℎ. It is actually 𝑧.
• λ -> Vector λ. It is actually 𝒢(𝑧).

Returns

• success -> Test that a factorization of the Pick matrix succeeded.

calc_phis(grid::AbstractGrid, Gᵥ::Vector{APC})

Try to calculate the Φ vector, which is used to calculate the 𝒜 matrix. Note that Φ should not be changed anymore once it has been established.

Arguments

• grid -> Grid in imaginary axis for input Green's function.
• Gᵥ -> Input Green's function.

Returns

• Φ -> Φ vector in Schur algorithm.

calc_abcd(grid::AbstractGrid, mesh::AbstractMesh, Φ::Vector{APC})

Try to calculate the coefficients matrix abcd (here it is called 𝒜), which is then used to calculate θ. See Eq. (8) in Fei's NAC paper.

Arguments

• grid -> Grid in imaginary axis for input Green's function.
• mesh -> Real frequency mesh.
• Φ -> Φ vector calculated by calc_phis().

Returns

• 𝒜 -> Coefficients matrix abcd in Schur algorithm.

calc_hbasis(z::APC, k::I64)

Try to calculate the Hardy basis $f^k(z)$.

Arguments

• z -> A complex variable.
• k -> Current order for the Hardy basis.

Returns

See above explanations.

calc_hmatrix(mesh::AbstractMesh, H::I64)

Try to calculate $[f^k(z), f^k(z)^*]$ for 0 ≤ 𝑘 ≤ 𝐻-1, which is called the hardy matrix (ℋ) and is used to evaluate θₘ₊₁.

Arguments

• mesh -> Real frequency mesh.
• H -> Maximum order for the Hardy basis.

Returns

• ℋ -> Hardy matrix for Hardy basis optimization.

calc_theta(𝒜::Array{APC,3}, ℋ::Array{APC,2}, 𝑎𝑏::Vector{C64})

Try to calculate the contractive function θ(z). 𝒜 is the coefficients matrix abcd, ℋ is the Hardy matrix, and 𝑎𝑏 are complex coefficients for expanding θₘ₊₁. See Eq. (7) in Fei's NAC paper.

Arguments

• 𝒜 -> Matrix 𝑎𝑏𝑐𝑑.
• ℋ -> Hardy matrix.
• 𝑎𝑏 -> Expansion coefficients 𝑎 and 𝑏 for the contractive function θ.

Returns

See above explanations.

calc_green(𝒜::Array{APC,3}, ℋ::Array{APC,2}, 𝑎𝑏::Vector{C64})

Firstly we try to calculate θ. Then θ is back transformed to a Nevanlinna interpolant via the inverse Mobius transform. Here, 𝒜 (abcd matrix), ℋ (Hardy matrix), and 𝑎𝑏 are used to evaluate θ.

Arguments

• 𝒜 -> Matrix 𝑎𝑏𝑐𝑑.
• ℋ -> Hardy matrix.
• 𝑎𝑏 -> Expansion coefficients 𝑎 and 𝑏 for the contractive function θ.

Returns

Gout -> Retarded Green's function, G(ω).

calc_noptim(ωₙ::Vector{APC}, Gₙ::Vector{APC})

Evaluate the optimal value for the size of input data (how may frequency points are actually used in the analytic continuation simulations) via the Pick criterion.

Arguments

• ωₙ -> Matsubara frequency points (the 𝑖 unit is not included).
• Gₙ -> Matsubara Green's function.

Returns

• ngrid -> Optimal number for the size of input data.

calc_hmin!(nac::NevanACContext)

Try to perform Hardy basis optimization. Such that the Hardy matrix ℋ and the corresponding coefficients 𝑎𝑏 are updated. They are used to calculate θ, which is then back transformed to generate smooth G (i.e., the spectrum) at real axis.

This function will determine the minimal value of H (hmin). Of course, ℋ and 𝑎𝑏 in NevanACContext struct are also changed.

Arguments

• nac -> A NevanACContext struct.

Returns

N/A

calc_hopt!(nac::NevanACContext)

Try to perform Hardy basis optimization. Such that the Hardy matrix ℋ and the corresponding coefficients 𝑎𝑏 are updated. They are used to calculate θ, which is then back transformed to generate smooth G (i.e., the spectrum) at real axis.

This function will determine the optimal value of H (hopt). Of course, ℋ and 𝑎𝑏 in NevanACContext struct are also changed. Here, H means order of the Hardy basis.

Arguments

• nac -> A NevanACContext struct.

Returns

N/A

hardy_optimize!(
    nac::NevanACContext,
    ℋ::Array{APC,2},
    𝑎𝑏::Vector{C64},
    H::I64
)

For given Hardy matrix ℋ, try to update the expanding coefficients 𝑎𝑏 by minimizing the smooth norm.

Arguments

• nac -> A NevanACContext struct.
• ℋ -> Hardy matrix, which contains the Hardy basis.
• 𝑎𝑏 -> Expansion coefficients 𝑎 and 𝑏 for the contractive function θ.
• H -> Current order of the Hardy basis.

Returns

• causality -> Test whether the solution is causal.
• converged -> Check whether the optimization is converged.

smooth_norm(nac::NevanACContext, ℋ::Array{APC,2}, 𝑎𝑏::Vector{C64})

Establish the smooth norm, which is used to improve the smoothness of the output spectrum. See Fei's paper for more details.

Arguments

• nac -> A NevanACContext struct.
• ℋ -> Hardy matrix, which contains the Hardy basis.
• 𝑎𝑏 -> Expansion coefficients 𝑎 and 𝑏 for the contractive function θ.

Returns

• 𝐹 -> Value of smooth norm.

check_pick(wn::Vector{APC}, gw::Vector{APC}, Nopt::I64)

Check whether the input data are valid (the Pick criterion is satisfied). Here, wn is the Matsubara frequency, gw is the Matsubara function, and Nopt is the optimized number of Matsubara data points.

Arguments

See above explanations.

Returns

N/A

check_causality(ℋ::Array{APC,2}, 𝑎𝑏::Vector{C64})

Check causality of the Hardy coefficients 𝑎𝑏.

Arguments

• ℋ -> Hardy matrix for Hardy basis optimization.
• 𝑎𝑏 -> Coefficients matrix for expanding θ with Hardy basis.

Returns

• causality -> Test whether the Hardy coefficients are causal.

StochACSolver

It represents the analytic continuation solver that implements the stochastic analytic continuation method (K. S. D. Beach's version).

This solver is moderately recommended. It is an alternative of the StochSK solver.

Features

• Only diagonal Green's functions are supported.
• Both fermionic and bosonic correlators are supported.
• Both continuous and discrete spectra are supported.
• Both imaginary-time and Matsubara data are supported.
• Good robustness with respect to noisy correlators.
• It supports the constrained sampling algorithm.
• Numerically stable.
• It is parallelized.

Known Limitations

• It is quite slow.
• It is less accurate than the BarRat solver at most cases.
• It doesn't support off-diagonal Green's functions.
• It doesn't support negative spectral weights.
• If there are multiple δ-like peaks, it is not good.

Implementation

See src/sac.jl.

Tests

See test/basic.

StochACMC

Mutable struct. It is used within the StochAC solver. It includes random number generator and some counters.

Members

• rng -> Random number generator.
• Macc -> Counter for move operation (accepted).
• Mtry -> Counter for move operation (tried).
• Sacc -> Counter for swap operation (accepted).
• Stry -> Counter for swap operation (tried).

See also: StochACSolver.

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.

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.

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(ω).

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.

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.

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.

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.

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(ω).

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

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

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

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.

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.

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.

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.

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, DynamicMesh.

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.

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.

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.

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.

calc_alpha()

Generate a list for the α parameters.

Arguments

N/A

Returns

• αₗ -> List of the α parameters.

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.

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.

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.

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.

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.

StochSKSolver

It represents the analytic continuation solver that implements the stochastic analytic continuation method (A. W. Sandvik's version).

This solver is moderately recommended. It is an alternative of the StochAC solver.

Features

• Only diagonal Green's functions are supported.
• Both fermionic and bosonic correlators are supported.
• Both continuous and discrete spectra are supported.
• Both imaginary-time and Matsubara data are supported.
• Good robustness with respect to noisy correlators.
• It supports the constrained sampling algorithm.
• Numerically stable.
• It is parallelized.

Known Limitations

• It is quite slow.
• It is less accurate than the BarRat solver at most cases.
• It doesn't support off-diagonal Green's functions.
• It doesn't support negative spectral weights.
• If there are multiple δ-like peaks, it is not good.

Implementation

See src/san.jl.

Tests

See test/basic.

StochSKMC

Mutable struct. It is used within the StochSK solver. It includes random number generator and some counters.

Members

• rng -> Random number generator.
• Sacc -> Counter for single-updated operation (accepted).
• Stry -> Counter for single-updated operation (tried).
• Pacc -> Counter for pair-updated operation (accepted).
• Ptry -> Counter for pair-updated operation (tried).
• Qacc -> Counter for quadruple-updated operation (accepted).
• Qtry -> Counter for quadruple-updated operation (tried).

See also: StochSKSolver.

StochSKElement

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

In the present implementation of StochSK solver, the amplitudes of the δ functions are fixed. But in principles, they could be sampled in the Monte Carlo procedure.

Members

• P -> It means the positions of the δ functions.
• A -> It means the weights / amplitudes of the δ functions.
• W -> It denotes the window that is used to shift the δ functions.

StochSKContext

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

Members

• Gᵥ -> Input data for correlator.
• Gᵧ -> Generated correlator.
• σ¹ -> Actually 1.0 / σ¹.
• allow -> Allowable indices.
• grid -> Imaginary axis grid for input data.
• mesh -> Real frequency mesh for output spectrum.
• kernel -> Default kernel function.
• Aout -> Calculated spectral function.
• χ² -> Current goodness-of-fit function.
• χ²min -> Mininum goodness-of-fit function.
• χ²vec -> Vector of goodness-of-fit function.
• Θ -> Current Θ parameter.
• Θvec -> Vector of Θ parameter.

solve(S::StochSKSolver, rd::RawData)

Solve the analytic continuation problem by using the stochastic analytic continuation algorithm (A. W. Sandvik's version). It is the driver for the StochSK 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 StochOM solvers.

Now the StochSK solver supports both continuous and δ-like spectra.

Arguments

• S -> A StochSKSolver 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(ω).

init(S::StochSKSolver, rd::RawData)

Initialize the StochSK solver and return the StochSKMC, StochSKElement, and StochSKContext structs. Please don't call this function directly.

Arguments

• S -> A StochSKSolver struct.
• rd -> A RawData struct, containing raw data for input correlator.

Returns

• MC -> A StochSKMC struct.
• SE -> A StochSKElement struct.
• SC -> A StochSKContext struct.

run(MC::StochSKMC, SE::StochSKElement, SC::StochSKContext)

Perform stochastic analytic continuation simulation, sequential version.

Arguments

• MC -> A StochSKMC struct.
• SE -> A StochSKElement struct.
• SC -> A StochSKContext struct.

Returns

• Aout -> Spectral function, A(ω).
• χ²vec -> Θ-dependent χ², χ²(Θ).
• Θvec -> List of Θ parameters.

prun(
    S::StochSKSolver,
    p1::Dict{String,Vector{Any}},
    p2::Dict{String,Vector{Any}},
    MC::StochSKMC,
    SE::StochSKElement,
    SC::StochSKContext
)

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

Arguments

• S -> A StochSKSolver struct.
• p1 -> A copy of PBASE.
• p2 -> A copy of PStochSK.
• MC -> A StochSKMC struct.
• SE -> A StochSKElement struct.
• SC -> A StochSKContext struct.

Returns

• Aout -> Spectral function, A(ω).
• χ²vec -> Θ-dependent χ², χ²(Θ).
• Θvec -> List of Θ parameters.

average(step::F64, SC::StochSKContext)

Postprocess the results generated during the stochastic analytic continuation simulations. It will generate the spectral functions.

Arguments

• step -> Number of Monte Carlo samplings.
• SC -> A StochSKContext struct.

Returns

• Aout -> Spectral function, A(ω).
• χ²vec -> Θ-dependent χ², χ²(Θ).
• Θvec -> List of Θ parameters.

last(
    SC::StochSKContext,
    Asum::Vector{F64},
    χ²vec::Vector{F64},
    Θvec::Vector{F64}
)

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

Arguments

• SC -> A StochSKContext struct.
• Asum -> Spectral function, A(ω).
• χ²vec -> Θ-dependent χ².
• Θvec -> List of Θ parameters.

Returns

• G -> Retarded Green's function, G(ω).

warmup(MC::StochSKMC, SE::StochSKElement, SC::StochSKContext)

Warmup the Monte Carlo engine to acheieve thermalized equilibrium. Then it will try to figure out the optimized Θ and the corresponding Monte Carlo field configuration.

Arguments

• MC -> A StochSKMC struct.
• SE -> A StochSKElement struct.
• SC -> A StochSKContext struct.

Returns

N/A

sample(MC::StochSKMC, SE::StochSKElement, SC::StochSKContext)

Perform Monte Carlo sweeps and sample the field configurations.

Arguments

• MC -> A StochSKMC struct.
• SE -> A StochSKElement struct.
• SC -> A StochSKContext struct.

Returns

N/A

measure(SE::StochSKElement, SC::StochSKContext)

Accumulate the final spectral functions A(ω).

Arguments

• SE -> A StochSKElement struct.
• SC -> A StochSKContext struct.

Returns

N/A

See also: nearest.

shuffle(MC::StochSKMC, SE::StochSKElement, SC::StochSKContext)

Try to shuffle the Monte Carlo field configuration via the Metropolis algorithm. Then the window for shifting the δ functions is adjusted.

Arguments

• MC -> A StochSKMC struct.
• SE -> A StochSKElement struct.
• SC -> A StochSKContext struct.

Returns

N/A

init_iodata(S::StochSKSolver, rd::RawData)

Preprocess the input data (rd).

Arguments

• S -> A StochSKSolver struct.
• rd -> A RawData struct, which contains essential input data.

Returns

• Gᵥ -> Input correlator.
• σ¹ -> 1.0 / σ¹.

See also: RawData.

init_mc(S::StochSKSolver)

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

Arguments

• S -> A StochSKSolver struct.

Returns

• MC -> A StochSKMC struct.

See also: StochSKMC.

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

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

Arguments

• S -> A StochSKSolver struct.
• rng -> Random number generator.
• allow -> Allowed positions for the δ peaks.

Returns

• SE -> A StochSKElement struct.

See also: StochSKElement.

init_context(
    SE::StochSKElement,
    Gᵥ::Vector{F64},
    σ¹::Vector{F64},
    allow::Vector{I64},
    grid::AbstractGrid,
    mesh::AbstractMesh,
    fmesh::AbstractMesh
)

Try to create a StochSKContext struct, which contains some key variables, including grid, mesh, input correlator and the corresponding standard deviation, kernel matrix, spectral function, and goodness-of-fit function.

Arguments

• SE -> A StochSKElement 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 StochSKContext struct.

calc_fmesh(S::StochSKSolver)

Try to calculate very fine (dense) linear mesh in [wmin, wmax], which is used internally to build the kernel function. Note that the stochastic analytic continuation method (A. W. Sandvik's version) does not support the self-adaptive mesh.

Arguments

• S -> A StochSKSolver struct.

Returns

• fmesh -> A very fine, uniform mesh in [wmin, wmax].

See also: LinearMesh, DynamicMesh.

calc_correlator(SE::StochSKElement, kernel::Array{F64,2})

Try to calculate correlator with the kernel function and the Monte Carlo field configuration. This correlator will then be used to evaluate the goodness-of-fit function χ².

Arguments

• SE -> A StochSKElement struct.
• kernel -> The fermionic or bosonic kernel.

Returns

• G -> Reconstructed correlator.

See also: calc_goodness.

calc_goodness(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 correlators.
• Gᵥ -> Input (original) correlators.

Returns

• χ² -> Goodness-of-fit function.

See also: calc_correlator.

calc_theta(𝒜::Array{APC,3}, ℋ::Array{APC,2}, 𝑎𝑏::Vector{C64})

Try to calculate the contractive function θ(z). 𝒜 is the coefficients matrix abcd, ℋ is the Hardy matrix, and 𝑎𝑏 are complex coefficients for expanding θₘ₊₁. See Eq. (7) in Fei's NAC paper.

Arguments

• 𝒜 -> Matrix 𝑎𝑏𝑐𝑑.
• ℋ -> Hardy matrix.
• 𝑎𝑏 -> Expansion coefficients 𝑎 and 𝑏 for the contractive function θ.

Returns

See above explanations.

calc_theta(len::I64, SC::StochSKContext)

Try to locate the optimal Θ and χ². This function implements the chi2min and chi2kink algorithms. Note that the chi2min algorithm is preferred.

Arguments

• len -> Length of vector Θ.
• SC -> A StochSKContext struct.

Returns

• c -> Selected index for optimal Θ.

constraints(S::StochSKSolver, 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, fmesh should be a fine linear mesh.

Arguments

• S -> A StochSKSolver struct.
• fmesh -> Very dense mesh for the δ peaks.

Returns

• allow -> Allowable indices.

See also: StochSKSolver.

try_move_s(MC::StochSKMC, SE::StochSKElement, SC::StochSKContext)

Try to update the Monte Carlo field configurations via the Metropolis algorithm. In each update, only single δ function is shifted.

Arguments

• MC -> A StochSKMC struct.
• SE -> A StochSKElement struct.
• SC -> A StochSKContext struct.

Returns

N/A

See also: try_move_p.

try_move_p(MC::StochSKMC, SE::StochSKElement, SC::StochSKContext)

Try to update the Monte Carlo field configurations via the Metropolis algorithm. In each update, only a pair of δ functions are shifted.

Arguments

• MC -> A StochSKMC struct.
• SE -> A StochSKElement struct.
• SC -> A StochSKContext struct.

Returns

N/A

See also: try_move_s.

try_move_q(MC::StochSKMC, SE::StochSKElement, SC::StochSKContext)

Try to update the Monte Carlo field configurations via the Metropolis algorithm. In each update, four different δ functions are shifted.

Arguments

• MC -> A StochSKMC struct.
• SE -> A StochSKElement struct.
• SC -> A StochSKContext struct.

Returns

N/A

See also: try_move_s.

StochOMSolver

It represents the analytic continuation solver that implements the stochastic optimization method.

This solver is less recommended.

Features

• Only diagonal Green's functions are supported.
• Both fermionic and bosonic correlators are supported.
• Only continuous spectra are supported.
• Good robustness with respect to noisy correlators.
• It supports the constrained sampling algorithm.
• Numerically stable.
• It is parallelized.

Known Limitations

• It is extremely slow.
• It is less accurate than the BarRat solver at most cases.
• It doesn't support off-diagonal Green's functions.
• It doesn't support negative spectral weights.
• If there are multiple δ-like peaks, it is not good.

Implementation

See src/som.jl.

Tests

See test/basic.

StochOMMC

Mutable struct. It is used within the StochOM solver. It includes random number generator and some counters.

Because the StochOM solver supports many Monte Carlo updates, so Macc and Mtry are vectors.

Members

• rng -> Random number generator.
• Macc -> Counter for move operation (accepted).
• Mtry -> Counter for move operation (tried).

See also: StochOMSolver.

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.

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.

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.

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(ω).

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.

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(ω).

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(ω).

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(ω).

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(ω).

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

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, GreenData.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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

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

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

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

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

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

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

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 δξ.

StochPXSolver

It represents the analytic continuation solver that implements the stochastic pole expansion method.

This solver is highly recommended.

Features

• Both off-diagonal and diagonal Green's functions are supported.
• Both fermionic and bosonic correlators are supported.
• Both continuous and discrete spectra are supported.
• The spectral weights can be negative.
• Good robustness with respect to noisy correlators.
• It supports the constrained sampling algorithm.
• Numerically stable.
• It can provide analytic expressions to approximate the correlators.
• It is parallelized.

Known Limitations

• It is quite slow.
• It is less accurate than the BarRat solver at most cases.
• It requires Matsubara data as input.
• If the spectra are discrete, a priori knowledge is essential.
• If the spectra exhibit negative weights, a priori knowledge is essential.

Implementation

See src/spx.jl.

Tests

See test/pole and test/lqcd.

StochPXMC

Mutable struct. It is used within the StochPX solver. It includes random number generator and some counters.

Members

• rng -> Random number generator.
• Sacc -> Counter for position-updated (type 1) operation (accepted).
• Stry -> Counter for position-updated (type 1) operation (tried).
• Pacc -> Counter for position-updated (type 2) operation (accepted).
• Ptry -> Counter for position-updated (type 2) operation (tried).
• Aacc -> Counter for amplitude-updated operation (accepted).
• Atry -> Counter for amplitude-updated operation (tried).
• Xacc -> Counter for exchange operation (accepted).
• Xtry -> Counter for exchange operation (tried).

See also: StochPXSolver.

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.

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.

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(ω).

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.

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ωₙ).

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ωₙ).

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ωₙ).

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

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

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

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, GreenData.

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.

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.

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.

reset_mc(MC::StochPXMC)

Reset the counters in StochPXMC struct.

Arguments

• MC -> A StochPXMC struct.

Returns

N/A

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

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

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, DynamicMesh.

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.

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.

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.

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ωₙ).

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ωₙ).

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(ω).

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(ω).

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.

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.

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.

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.

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.

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.