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

#=
### *Customized Structs* : *MaxEnt Solver*
=#

"""
    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.
"""
mutable struct MaxEntContext
    Gáµ¥     :: Vector{F64}
    ÏƒÂ²     :: Vector{F64}
    grid   :: AbstractGrid
    mesh   :: AbstractMesh
    model  :: Vector{F64}
    kernel :: Array{F64,2}
    hess   :: Array{F64,2}
    Vâ‚›     :: Array{F64,2}
    Wâ‚‚     :: Array{F64,2}
    Wâ‚ƒ     :: Array{F64,3}
    Bâ‚˜     :: Vector{F64}
end

#=
### *Global Drivers*
=#

"""
    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(Ï‰).
"""
function solve(S::MaxEntSolver, rd::RawData)
    println("[ MaxEnt ]")
    #
    mec = init(S, rd)
    darr, sol = run(mec)
    gout = last(mec, darr, sol)
    #
    return mec.mesh.mesh, sol[:A], gout
end

"""
    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.
"""
function init(S::MaxEntSolver, rd::RawData)
    # Prepera input data
    G = make_data(rd)
    Gáµ¥ = G.value
    ÏƒÂ² = 1.0 ./ G.covar
    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")

    # Prepare default model function
    model = make_model(mesh)
    println("Build default model: ", get_b("mtype"))

    # Prepare kernel function
    kernel = make_kernel(mesh, grid)
    println("Build default kernel: ", get_b("ktype"))

    # Prepare some essential intermediate variables
    Vâ‚›, Wâ‚‚, Wâ‚ƒ, Bâ‚˜, hess = precompute(Gáµ¥, ÏƒÂ², mesh, model, kernel)
    println("Precompute key coefficients")

    return MaxEntContext(Gáµ¥, ÏƒÂ², grid, mesh, model,
                         kernel, hess, Vâ‚›, Wâ‚‚, Wâ‚ƒ, Bâ‚˜)
end

"""
    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.
"""
function run(mec::MaxEntContext)
    stype = get_m("stype")
    method = get_m("method")

    # Note that the Bayesian Reconstruction entropy is compatible with
    # all the four algorithms so far.
    if stype == "br"
        println("Bayesian Reconstruction entropy is used!")
    else
        println("Shannonâ€“Jaynes entropy is used!")
    end

    @cswitch method begin
        @case "historic"
            return historic(mec)
            break

        @case "classic"
            return classic(mec)
            break

        @case "bryan"
            return bryan(mec)
            break

        @case "chi2kink"
            return chi2kink(mec)
            break
    end
end

"""
    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(Ï‰).
"""
function last(mec::MaxEntContext, svec::Vector, sol::Dict)
    # By default, we should write the analytic continuation results
    # into the external files.
    _fwrite = get_b("fwrite")
    fwrite = isa(_fwrite, Missing) || _fwrite ? true : false

    # Write the spectral function
    fwrite && write_spectrum(mec.mesh, sol[:A])

    # Write the model function
    fwrite && write_model(mec.mesh, mec.model)

    # Write Î±-Ï‡Â² data
    Î±_vec = map(x -> x[:Î±], svec)
    Ï‡_vec = map(x -> x[:Ï‡Â²], svec)
    fwrite && write_misfit(Î±_vec, Ï‡_vec)

    # Write P[Î±|A] for bryan algorithm
    if haskey(svec[end], :prob)
        p_vec = map(x -> x[:prob], svec)
        fwrite && write_probability(Î±_vec, p_vec)
    end

    # Regenerate the input data and write them
    #
    # Perhaps the spectral function (sol[:A]) is blurred.
    # But we need the pure spectral function (sol[:Araw]).
    Aout = haskey(sol, :Araw) ? sol[:Araw] : sol[:A]
    G = reprod(mec.mesh, mec.kernel, Aout)
    fwrite && write_backward(mec.grid, G)

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

    return _G
end

#=
### *Core Algorithms*
=#

"""
    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`](@ref).
"""
function historic(mec::MaxEntContext)
    function root_fun(_Î±, _u)
        res = optimizer(mec, _Î±, _u, use_bayes)
        @. _u = res[:u]
        return length(mec.ÏƒÂ²) / res[:Ï‡Â²] - 1.0
    end

    println("Apply historic algorithm to determine optimized Î±")

    use_bayes = false
    alpha = get_m("alpha")
    ratio = get_m("ratio")
    n_svd = length(mec.Bâ‚˜)

    u_vec = zeros(F64, n_svd)
    s_vec = []

    conv = 0.0
    while conv < 1.0
        sol = optimizer(mec, alpha, u_vec, use_bayes)
        push!(s_vec, sol)
        alpha = alpha / ratio
        conv = length(mec.ÏƒÂ²) / sol[:Ï‡Â²]
    end

    u_vec = s_vec[end-1][:u]
    alpha = s_vec[end][:Î±]
    Î±_opt = secant(root_fun, alpha, u_vec)

    sol = optimizer(mec, Î±_opt, u_vec, use_bayes)
    println("Optimized Î± : $Î±_opt log10(Î±) : $(log10(Î±_opt))")

    return s_vec, sol
end

"""
    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`](@ref).
"""
function classic(mec::MaxEntContext)
    function root_fun(_Î±, _u)
        res = optimizer(mec, _Î±, _u, use_bayes)
        @. _u = res[:u]
        return res[:conv] - 1.0
    end

    println("Apply classic algorithm to determine optimized Î±")

    use_bayes = true
    alpha = get_m("alpha")
    ratio = get_m("ratio")
    n_svd = length(mec.Bâ‚˜)

    u_vec = zeros(F64, n_svd)
    s_vec = []

    conv = 0.0
    while conv < 1.0
        sol = optimizer(mec, alpha, u_vec, use_bayes)
        push!(s_vec, sol)
        alpha = alpha / ratio
        @. u_vec = sol[:u]
        conv = sol[:conv]
    end

    c_vec = [x[:conv] for x in s_vec]
    Î±_vec = [x[:Î±] for x in s_vec]
    exp_opt = log10(Î±_vec[end] / Î±_vec[end-1])
    exp_opt = exp_opt / log10(c_vec[end] / c_vec[end-1])
    exp_opt = log10(Î±_vec[end-1]) - log10(c_vec[end-1]) * exp_opt

    # Starting from the predicted value of Î±, and starting optimization
    # at the solution for the next-lowest Î±, we find the optimal Î± by
    # secant root finding method.
    u_vec = s_vec[end-1][:u]
    alpha = 10.0 ^ exp_opt
    Î±_opt = secant(root_fun, alpha, u_vec)

    sol = optimizer(mec, Î±_opt, u_vec, use_bayes)
    println("Optimized Î± : $Î±_opt log10(Î±) : $(log10(Î±_opt))")

    return s_vec, sol
end

"""
    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`](@ref).
"""
function bryan(mec::MaxEntContext)
    println("Apply bryan algorithm to determine optimized Î±")

    use_bayes = true
    alpha = get_m("alpha")
    ratio = get_m("ratio")
    n_svd = length(mec.Bâ‚˜)
    nmesh = length(mec.mesh)

    u_vec = zeros(F64, n_svd)
    s_vec = []

    maxprob = 0.0
    while true
        sol = optimizer(mec, alpha, u_vec, use_bayes)
        push!(s_vec, sol)
        alpha = alpha / ratio
        @. u_vec = sol[:u]
        prob = sol[:prob]
        if prob > maxprob
            maxprob = prob
        elseif prob < 0.01 * maxprob
            break
        end
    end

    Î±_vec = map(x->x[:Î±], s_vec)
    p_vec = map(x->x[:prob], s_vec)
    p_vec = -p_vec ./ trapz(Î±_vec, p_vec)
    A_vec = map(x->x[:A], s_vec)

    nprob = length(p_vec)
    A_opt = zeros(F64, nmesh)
    spectra = zeros(F64, nmesh, nprob)
    for i = 1:nprob
        spectra[:,i] = A_vec[i] * p_vec[i]
    end
    for j = 1:nmesh
        A_opt[j] = -trapz(Î±_vec, spectra[j,:])
    end

    sol = Dict(:A => A_opt)

    return s_vec, sol
end

"""
    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`](@ref).
"""
function chi2kink(mec::MaxEntContext)
    function fitfun(x, p)
        return @. p[1] + p[2] / (1.0 + exp(-p[4] * (x - p[3])))
    end

    println("Apply chi2kink algorithm to determine optimized Î±")

    use_bayes = false
    alpha = get_m("alpha")
    ratio = get_m("ratio")
    nalph = get_m("nalph")
    Î±_end = alpha / (ratio^nalph)
    n_svd = length(mec.Bâ‚˜)

    u_vec = zeros(F64, n_svd)
    s_vec = []
    Ï‡_vec = []
    Î±_vec = []

    while true
        sol = optimizer(mec, alpha, u_vec, use_bayes)
        push!(s_vec, sol)
        push!(Î±_vec, alpha)
        push!(Ï‡_vec, sol[:Ï‡Â²])
        @. u_vec = sol[:u]
        alpha = alpha / ratio
        if alpha < Î±_end
            break
        end
    end

    good = isfinite.(Ï‡_vec)
    guess = [0.0, 5.0, 2.0, 0.0]
    fit = curve_fit(fitfun, log10.(Î±_vec[good]), log10.(Ï‡_vec[good]), guess)
    _, _, c, d = fit.param

    # `fit_pos` is a control parameter for under/overfitting.
    # Good values are usually between 2 and 2.5. Smaller values usually
    # lead to underfitting, which is sometimes desirable. Larger values
    # lead to overfitting, which should be avoided.
    fit_pos = 2.5
    Î±_opt = c - fit_pos / d
    close = argmin( abs.( log10.(Î±_vec) .- Î±_opt ) )
    u_vec = s_vec[close][:u]
    Î±_opt = 10.0 ^ Î±_opt

    sol = optimizer(mec, Î±_opt, u_vec, use_bayes)
    println("Optimized Î± : $Î±_opt log10(Î±) : $(log10(Î±_opt))")

    return s_vec, sol
end

"""
    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.
"""
function optimizer(
    mec::MaxEntContext,
    Î±::F64,
    us::Vector{F64},
    use_bayes::Bool
    )
    blur = get_m("blur")
    offdiag = get_b("offdiag")

    if offdiag
        solution, call = newton(f_and_J_od, us, mec, Î±)
        u = copy(solution)
        A = svd_to_real_od(mec, solution)
        S = calc_entropy_od(mec, A)
    else
        solution, call = newton(f_and_J, us, mec, Î±)
        u = copy(solution)
        A = svd_to_real(mec, solution)
        S = calc_entropy(mec, A, u)
    end

    Ï‡Â² = calc_chi2(mec, A)
    norm = trapz(mec.mesh, A)

    dict = Dict{Symbol,Any}(
        :u => u,
        :Î± => Î±,
        :S => S,
        :Ï‡Â² => Ï‡Â²,
        :norm => norm,
        :Q => Î± * S - 0.5 * Ï‡Â²,
        :Araw => deepcopy(A),
    )

    if use_bayes
        if offdiag
            ng, tr, conv, prob = calc_bayes_od(mec, A, S, Ï‡Â², Î±)
        else
            ng, tr, conv, prob = calc_bayes(mec, A, S, Ï‡Â², Î±)
        end
        dict[:ngood] = ng
        dict[:trace] = tr
        dict[:conv] = conv
        dict[:prob] = prob
    end

    if blur > 0.0
        make_blur(mec.mesh, A, blur)
    end
    dict[:A] = A

    @printf("log10(Î±) = %8.4f ", log10(Î±))
    @printf("Ï‡Â² = %8.4e ", Ï‡Â²)
    @printf("S = %8.4e ", S)
    @printf("call = %4i ", call)
    @printf("norm = %8.4f\n", norm)

    return dict
end

#=
### *Service Functions*
=#

#=
*Remarks* :

Try to calculate some key variables by using the Einstein summation trick.

```math
\begin{equation}
B_m = \sum^{N}_{n = 1} \frac{1}{\sigma^2_n} \xi_m U_{nm} G_n,
\end{equation}
```

```math
\begin{equation}
W_{ml} = \sum_{pn} \frac{1}{\sigma^2_n}
    U_{nm}\xi_m U_{np} \xi_p V_{lp} \Delta_l D_l,
\end{equation}
```

```math
\begin{equation}
W_{mli} = W_{ml} V_{li}.
\end{equation}
```

Note that these variables do not depend on the spectral function
``A(\omega)``, so they could be computed at advance to improve the
computational efficiency.

---

The `hessian matrix` is also calculated here.

```math
\begin{equation}
L = \frac{1}{2} \chi^2,
\end{equation}
```

```math
\begin{equation}
\frac{\partial^2 L}{\partial A_i \partial A_j} =
\sum_n \frac{K_{ni} \Delta_i K_{nj} \Delta_j}{\sigma^2_n}.
\end{equation}
```
=#

"""
    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.
"""
function precompute(
    Gáµ¥::Vector{F64},
    ÏƒÂ²::Vector{F64},
    am::AbstractMesh,
    D::Vector{F64},
    K::Matrix{F64}
    )
    # Create singular value space
    U, V, S = make_singular_space(K)

    # Evaluate sizes of the arrays
    nmesh = length(am)
    n_svd = length(S)

    # Allocate memories
    Wâ‚‚ = zeros(F64, n_svd, nmesh)
    Wâ‚ƒ = zeros(F64, n_svd, n_svd, nmesh)
    Bâ‚˜ = zeros(F64, n_svd)
    hess = zeros(F64, nmesh, nmesh)

    # Get weight of the mesh, Î”Ï‰â‚—.
    Î” = am.weight

    # Compute Wâ‚˜â‚—
    @einsum Wâ‚‚[m,l] = ÏƒÂ²[k] * U[k,m] * S[m] * U[k,n] * S[n] * V[l,n] * Î”[l] * D[l]

    # Compute Wâ‚˜â‚—áµ¢
    @einsum Wâ‚ƒ[m,k,l] = Wâ‚‚[m,l] * V[l,k]

    # Compute Bâ‚˜
    @einsum Bâ‚˜[m] = S[m] * U[k,m] * ÏƒÂ²[k] * Gáµ¥[k]

    # Compute the Hessian matrix
    @einsum hess[i,j] = Î”[i] * Î”[j] * K[k,i] * K[k,j] * ÏƒÂ²[k]

    return V, Wâ‚‚, Wâ‚ƒ, Bâ‚˜, hess
end

#=
*Remarks* :

For Shannon-Jaynes entropy,

```math
\begin{equation}
w_l = \exp \left(\sum_m V_{lm} u_m\right).
\end{equation}
```

```math
\begin{equation}
f_m = \alpha u_m + \sum_l W_{ml} w_l - B_m.
\end{equation}
```

```math
\begin{equation}
J_{mi} = \alpha \delta_{mi} + \sum_l W_{mli} w_l.
\end{equation}
```

---

For Bayesian Reconstruction entropy,

```math
\begin{equation}
w_l = \frac{1}{1 - D_l \sum_m V_{lm} u_m}.
\end{equation}
```

```math
\begin{equation}
f_m = \alpha u_m + \sum_l W_{ml} w_l - B_m.
\end{equation}
```

```math
\begin{equation}
J_{mi} = \alpha \delta_{mi} + \sum_l W_{mli} D_l w_l w_l.
\end{equation}
```
=#

"""
    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`](@ref).
"""
function f_and_J(u::Vector{F64}, mec::MaxEntContext, Î±::F64)
    stype = get_m("stype")

    n_svd = length(mec.Bâ‚˜)
    J = diagm([Î± for i = 1:n_svd])

    # For Shannonâ€“Jaynes entropy
    if stype == "sj"
        w = exp.(mec.Vâ‚› * u)
        #
        for j = 1:n_svd
            for i = 1:n_svd
                J[i,j] = J[i,j] + dot(mec.Wâ‚ƒ[i,j,:], w)
            end
        end
        #
        f = Î± * u + mec.Wâ‚‚ * w - mec.Bâ‚˜
    # For Bayesian Reconstruction entropy
    else
        w = mec.Vâ‚› * u
        wâ‚� = 1.0 ./ (1.0 .- mec.model .* w)
        wâ‚‚ = wâ‚� .* wâ‚� .* mec.model
        #
        for j = 1:n_svd
            for i = 1:n_svd
                J[i,j] = J[i,j] + dot(mec.Wâ‚ƒ[i,j,:], wâ‚‚)
            end
        end
        #
        f = Î± * u + mec.Wâ‚‚ * wâ‚� - mec.Bâ‚˜
    end

    return f, J
end

#=
*Remarks* :

For Shannon-Jaynes entropy,

```math
\begin{equation}
w_l = \exp \left(\sum_m V_{lm} u_m\right).
\end{equation}
```

```math
\begin{equation}
f_m = \alpha u_m +
      \sum_l W_{ml}\left(w_l - \frac{1}{w_l}\right) - B_m.
\end{equation}
```

```math
\begin{equation}
J_{mi} = \alpha \delta_{mi} +
         \sum_{l} W_{mli} \left(w_l + \frac{1}{w_l}\right).
\end{equation}
```

---

For Bayesian Reconstruction entropy,

```math
\begin{equation}
w^+_l = \frac{1}{ 1 - D_l \sum_m V_{lm} u_m}.
\end{equation}
```

```math
\begin{equation}
w^-_l = \frac{1}{ 1 + D_l \sum_m V_{lm} u_m}.
\end{equation}
```

```math
\begin{equation}
f_m = \alpha u_m + \sum_l W_{ml} (w^+_l - w^-_l) - B_m.
\end{equation}
```

```math
\begin{equation}
J_{mi} = \alpha \delta_{mi} + \sum_l W_{mli} D_l (w^+_l w^+_l + w^-_l w^-_l).
\end{equation}
```
=#

"""
    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`](@ref).
"""
function f_and_J_od(u::Vector{F64}, mec::MaxEntContext, Î±::F64)
    stype = get_m("stype")

    n_svd = length(mec.Bâ‚˜)
    J = diagm([Î± for i = 1:n_svd])

    # For Shannonâ€“Jaynes entropy
    if stype == "sj"
        w = exp.(mec.Vâ‚› * u)
        #
        aâ�º = 1.0 .* w
        aâ�» = 1.0 ./ w
        aâ‚� = aâ�º - aâ�»
        aâ‚‚ = aâ�º + aâ�»
        #
        for j = 1:n_svd
            for i = 1:n_svd
                J[i,j] = J[i,j] + dot(mec.Wâ‚ƒ[i,j,:], aâ‚‚)
            end
        end
        #
        f = Î± * u + mec.Wâ‚‚ * aâ‚� - mec.Bâ‚˜
    # For Bayesian Reconstruction entropy
    else
        w = mec.Vâ‚› * u
        #
        aâ�º = 1.0 ./ (1.0 .- mec.model .* w)
        aâ�» = 1.0 ./ (1.0 .+ mec.model .* w)
        aâ‚� = aâ�º - aâ�»
        aâ‚‚ = (aâ�º .* aâ�º + aâ�» .* aâ�») .* mec.model
        #
        for j = 1:n_svd
            for i = 1:n_svd
                J[i,j] = J[i,j] + dot(mec.Wâ‚ƒ[i,j,:], aâ‚‚)
            end
        end
        #
        f = Î± * u + mec.Wâ‚‚ * aâ‚� - mec.Bâ‚˜
    end

    return f, J
end

#=
*Remarks* :

For Shannon-Jaynes entropy,

```math
\begin{equation}
A_l = D_l \exp \left(\sum_m V_{lm} u_m\right).
\end{equation}
```

---

For Bayesian Reconstruction entropy,

```math
\begin{equation}
A_l = \frac{D_l}{ 1 - D_l \sum_m V_{lm} u_m}.
\end{equation}
```
=#

"""
    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`](@ref).
"""
function svd_to_real(mec::MaxEntContext, u::Vector{F64})
    stype = get_m("stype")
    #
    # For Shannonâ€“Jaynes entropy
    if stype == "sj"
        w = exp.(mec.Vâ‚› * u)
        return mec.model .* w
    # For Bayesian Reconstruction entropy
    else
        w = mec.Vâ‚› * u
        return mec.model ./ (1.0 .- mec.model .* w)
    end
end

#=
*Remarks* :

For Shannon-Jaynes entropy,

```math
\begin{equation}
A_l = D_l \exp \left(\sum_m V_{lm} u_m\right) -
      D_l \exp \left(-\sum_m V_{lm} u_m\right).
\end{equation}
```

---

For Bayesian Reconstruction entropy,

```math
\begin{equation}
w^+_l = \frac{1}{ 1 - D_l \sum_m V_{lm} u_m}.
\end{equation}
```

```math
\begin{equation}
w^-_l = \frac{1}{ 1 + D_l \sum_m V_{lm} u_m}.
\end{equation}
```

```math
\begin{equation}
A_l = D_l (w^+_l - w^-_l).
\end{equation}
```
=#

"""
    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`](@ref).
"""
function svd_to_real_od(mec::MaxEntContext, u::Vector{F64})
    stype = get_m("stype")
    #
    # For Shannonâ€“Jaynes entropy
    if stype == "sj"
        w = exp.(mec.Vâ‚› * u)
        wâ�º = w
        wâ�» = 1.0 ./ w
        return mec.model .* (wâ�º .- wâ�»)
    # For Bayesian Reconstruction entropy
    else
        w = mec.Vâ‚› * u
        wâ�º = 1.0 ./ (1.0 .- mec.model .* w)
        wâ�» = 1.0 ./ (1.0 .+ mec.model .* w)
        return mec.model .* (wâ�º .- wâ�»)
    end
end

#=
*Remarks* :

Shannonâ€“Jaynes entropy

```math
\begin{equation}
S[A] = \int^{\infty}_0 d\omega
\left[
    A - m -
    A \log{\left(\frac{A}{m}\right)}
\right],
\end{equation}
```

```math
\begin{equation}
S[A^{+},A^{-}] = \int^{+\infty}_0 d\omega
\left[
    \sqrt{A^2 + 4m^2} - 2m -
    A\log{\left(\frac{\sqrt{A^2 + 4m^2} + A}{2m}\right)}
\right].
\end{equation}
```

---

Bayesian Reconstruction entropy

```math
\begin{equation}
S[A] = \int^{\infty}_0 d\omega
\left[
    1 - \frac{A}{m} + \log{\left(\frac{A}{m}\right)}
\right],
\end{equation}
```

```math
\begin{equation}
S[A^{+},A^{-}] = \int^{+\infty}_0 d\omega
\left[
    2 - \frac{\sqrt{A^2 + m^2} + m}{m} +
    \log{\left(\frac{\sqrt{A^2 + m^2} + m}{2m}\right)}
\right].
\end{equation}
```
=#

"""
    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`](@ref).
"""
function calc_entropy(mec::MaxEntContext, A::Vector{F64}, u::Vector{F64})
    stype = get_m("stype")
    #
    # For Shannonâ€“Jaynes entropy
    if stype == "sj"
        f = A - mec.model - A .* (mec.Vâ‚› * u)
    # For Bayesian Reconstruction entropy
    else
        ð�‘… = A ./ mec.model
        #
        if any(x -> x < 0.0, ð�‘…)
            @info "Negative spectrum occurs!"
            @info "The results might be questionable."
            @info "Perhaps you should switch to the Shannonâ€“Jaynes entropy."
            f = 1.0 .- ð�‘… + log.(abs.(ð�‘…))
        else
            f = 1.0 .- ð�‘… + log.(ð�‘…)
        end
    end
    #
    return trapz(mec.mesh, f)
end

"""
    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`](@ref).
"""
function calc_entropy_od(mec::MaxEntContext, A::Vector{F64})
    stype = get_m("stype")
    #
    # For Shannonâ€“Jaynes entropy
    if stype == "sj"
        root = sqrt.(A .^ 2.0 + 4.0 .* mec.model .* mec.model)
        f = root - 2.0 .* mec.model
        f = f - A .* log.((root + A) ./ (2.0 .* mec.model))
    # For Bayesian Reconstruction entropy
    else
        root = sqrt.(A .^ 2.0 + mec.model .^ 2.0) + mec.model
        f = 2.0 .- (root ./ mec.model) + log.(root ./ (2.0 .* mec.model))
    end
    #
    return trapz(mec.mesh, f)
end

#=
*Remarks* :

**Posterior distribution of ``\alpha``**

Because
```math
\begin{equation}
\text{Pr}[\alpha | \bar{G}] =
\text{Pr}[\alpha] \frac{e^Q}{Z_L Z_S}
\frac{(2\pi)^{N/2}}{\sqrt{\det{[\alpha I + \Lambda]}}},
\end{equation}
```

so
```math
\begin{equation}
\log \text{Pr}[\alpha | \bar{G}] =
\text{constant} + \log \text{Pr} [\alpha] +
\frac{1}{2} \text{Tr} \log \left[\frac{\alpha I}{\alpha I + A}\right] +
\alpha S - \frac{1}{2}\chi^2.
\end{equation}
```

The defining equation for the `classic MaxEnt` equation reads:

```math
-2\alpha S = \text{Tr} \left(\frac{\Lambda}{\alpha I + \Lambda}\right).
```

The summation on the right hand side of the above equation is defined to
be ``N_g``, the number of good measurements:

```math
\begin{equation}
N_g = \text{Tr} \left(\frac{\Lambda}{\alpha I + \Lambda}\right).
\end{equation}
```

If ``\lambda_i`` are the eigenvalues of ``\Lambda``, then

```math
\begin{equation}
N_g = \sum_i \frac{\lambda_i}{\alpha + \lambda_i}.
\end{equation}
```

**``\Lambda`` matrix**

For Shannon-Jaynes entropy,

```math
\begin{equation}
\frac{\partial^2 S[A]}{\partial A_i \partial A_j} =
-\delta_{ij}\frac{\Delta_i}{A_i} =
-\delta_{ij}\frac{\sqrt{\Delta_i \Delta_j}}{\sqrt{A_i A_j}},
\end{equation}
```

```math
\begin{equation}
\Lambda_{ij} = \sqrt{\frac{A_i}{\Delta_i}}
               \frac{\partial^2 L}{\partial A_i \partial A_j}
               \sqrt{\frac{A_j}{\Delta_j}}.
\end{equation}
```

```math
\begin{equation}
\frac{\partial^2 S[A^+,A^-]}{\partial A_i \partial A_j} =
-\delta_{ij}\frac{\Delta_i}{\sqrt{A_i^2 + 4m_i^2}} =
-\delta_{ij}\frac{\sqrt[4]{\Delta^2_i}\sqrt[4]{\Delta^2_j}}
                 {\sqrt[4]{A_i^2 + 4m_i^2} \sqrt[4]{A_j^2 + 4m_j^2}},
\end{equation}
```

```math
\begin{equation}
\Lambda_{ij} = \sqrt[4]{\frac{A^2_i + 4m_i^2}{\Delta^2_i}}
               \frac{\partial^2 L}{\partial A_i \partial A_j}
               \sqrt[4]{\frac{A^2_j + 4m_j^2}{\Delta^2_j}}.
\end{equation}
```

---

For Bayesian Reconstruction entropy,

```math
\begin{equation}
\frac{\partial^2 S[A]}{\partial A_i \partial A_j} =
-\delta_{ij} \frac{\Delta_i}{A^2_i} =
-\delta_{ij} \frac{\sqrt{\Delta_i \Delta_j}}{A_i A_j},
\end{equation}
```

```math
\begin{equation}
\Lambda_{ij} = \frac{A_i}{\sqrt{\Delta_i}}
               \frac{\partial^2 L}{\partial A_i \partial A_j}
               \frac{A_j}{\sqrt{\Delta_j}}.
\end{equation}
```

```math
\begin{equation}
\frac{\partial^2 S[A^+,A^-]}{\partial A_i \partial A_j} =
-\delta_{ij} X_{ij} Y_{ij},
\end{equation}
```

```math
\begin{equation}
X_{ij} = \frac{\sqrt{2\Delta_i} \sqrt{2\Delta_j}}{
        \left(\sqrt{A^2_i + m^2_i} + m_i + A_i\right)
        \left(\sqrt{A^2_j + m^2_j} + m_j + A_j\right)
    },
\end{equation}
```

```math
\begin{equation}
Y_{ij} =
    \frac{
        \sqrt{A_i + \sqrt{A^2_i + m^2_i}}
        \sqrt{A_j + \sqrt{A^2_j + m^2_j}}
    }{
        \sqrt[4]{A^2_i + m^2_i}
        \sqrt[4]{A^2_j + m^2_j}
    },
\end{equation}
```

```math
\begin{equation}
\Lambda_{ij} = Z_i \frac{\partial^2 L}{\partial A_i \partial A_j} Z_j
\end{equation}
```

```math
\begin{equation}
Z_i = \frac{\left(\sqrt{A^2_i + m^2_i} + m_i + A_i\right)}{\sqrt{2\Delta_i}}
      \frac{\sqrt[4]{A^2_i + m^2_i}}{\sqrt{A_i + \sqrt{A^2_i + m^2_i}}}
\end{equation}
```

```math
\begin{equation}
Z_j = \frac{\left(\sqrt{A^2_j + m^2_j} + m_j + A_j\right)}{\sqrt{2\Delta_j}}
      \frac{\sqrt[4]{A^2_j + m^2_j}}{\sqrt{A_j + \sqrt{A^2_j + m^2_j}}}
\end{equation}
```

**Reference:**

[1] G. J. Kraberger, *et al.*, Phys. Rev. B **96**, 155128 (2017).

[2] M. Jarrell, *et al.*, Phys. Rep. **269**, 133 (1996).
=#

"""
    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`](@ref).
"""
function calc_bayes(
    mec::MaxEntContext,
    A::Vector{F64},
    S::F64,
    Ï‡Â²::F64,
    Î±::F64
    )
    stype = get_m("stype")
    mesh = mec.mesh

    if stype == "sj"
        T = sqrt.(A ./ mesh.weight)
    else
        T = A ./ sqrt.(mesh.weight)
    end
    Î› = (T * T') .* mec.hess

    Î» = eigvals(Hermitian(Î›))
    ng = -2.0 * Î± * S
    tr = sum(Î» ./ (Î± .+ Î»))
    conv = tr / ng

    eig_sum = sum(log.(Î± ./ (Î± .+ Î»)))
    log_prob = Î± * S - 0.5 * Ï‡Â² + log(Î±) + 0.5 * eig_sum

    return ng, tr, conv, exp(log_prob)
end

"""
    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`](@ref).
"""
function calc_bayes_od(
    mec::MaxEntContext,
    A::Vector{F64},
    S::F64,
    Ï‡Â²::F64,
    Î±::F64
    )
    stype = get_m("stype")
    mesh = mec.mesh

    if stype == "sj"
        R = (A .^ 2.0 + 4.0 * mec.model .^ 2.0) ./ (mesh.weight .^ 2.0)
        T = R .^ 0.25
    else
        R = sqrt.(A .^ 2.0 + mec.model .^ 2.0)
        X = (R .+ mec.model .+ A) ./ sqrt.(2.0 * mesh.weight)
        Y = sqrt.(R) ./ sqrt.(A .+ R)
        T = X .* Y
    end
    Î› = (T * T') .* mec.hess

    Î» = eigvals(Hermitian(Î›))
    ng = -2.0 * Î± * S
    tr = sum(Î» ./ (Î± .+ Î»))
    conv = tr / ng

    eig_sum = sum(log.(Î± ./ (Î± .+ Î»)))
    log_prob = Î± * S - 0.5 * Ï‡Â² + log(Î±) + 0.5 * eig_sum

    return ng, tr, conv, exp(log_prob)
end

"""
    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.
"""
function calc_chi2(mec::MaxEntContext, A::Vector{F64})
    Gâ‚™ = reprod(mec.mesh, mec.kernel, A)
    Ï‡Â² = sum(mec.ÏƒÂ² .* ((mec.Gáµ¥ - Gâ‚™) .^ 2.0))
    return Ï‡Â²
end
