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

#
# Note:
#
# The following codes for the NevanAC solver are mostly adapted from
#
#     https://github.com/SpM-lab/Nevanlinna.jl
#
# See
#
#     Nevanlinna.jl: A Julia implementation of Nevanlinna analytic continuation
#     Kosuke Nogaki, Jiani Fei, Emanuel Gull, Hiroshi Shinaoka
#     SciPost Phys. Codebases 19 (2023)
#
# for more details. And we thank Dr. Shuang Liang for her help.
#

#=
### *Customized Structs* : *NevanAC Solver*
=#

"""
    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.
"""
mutable struct NevanACContext
    Gáµ¥   :: Vector{APC}
    grid :: AbstractGrid
    mesh :: AbstractMesh
    Î¦    :: Vector{APC}
    ğ�’œ    :: Array{APC,3}
    â„‹    :: Array{APC,2}
    ğ�‘�ğ�‘�   :: Vector{C64}
    hmin :: I64
    hopt :: I64
end

#=
### *Global Drivers*
=#

"""
    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(Ï‰).
"""
function solve(S::NevanACSolver, rd::RawData)
    println("[ NevanAC ]")
    #
    nac = init(S, rd)
    run(nac)
    Aout, Gout = last(nac)
    #
    return nac.mesh.mesh, Aout, Gout
end

"""
    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.
"""
function init(S::NevanACSolver, rd::RawData)
    # Setup numerical precision. Note that the NAC method is extremely
    # sensitive to the float point precision.
    setprecision(128)

    # Convert the input data to APC, i.e., Complex{BigFloat}.
    Ï‰â‚™ = APC.(rd._grid * im)
    Gâ‚™ = APC.(rd.value)

    # Evaluate the optimal value for the size of input data.
    # Here we just apply the Pick criterion.
    ngrid = calc_noptim(Ï‰â‚™, Gâ‚™)

    # Prepera input data
    Gáµ¥ = calc_mobius(-Gâ‚™[1:ngrid])
    reverse!(Gáµ¥)
    println("Postprocess input data: ", length(Gáµ¥), " points")

    # Prepare grid for input data
    grid = make_grid(rd, T = APF)
    resize!(grid, ngrid)
    reverse!(grid)
    println("Build grid for input data: ", length(grid), " points")

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

    # Precompute key quantities to accelerate the computation
    Î¦, ğ�’œ, â„‹, ğ�‘�ğ�‘� = precompute(grid, mesh, Gáµ¥)
    println("Precompute key matrices")

    # Actually, the NevanACContext struct already contains enough
    # information to build the Nevanlinna interpolant and get the
    # spectrum, but Hardy basis optimization is needed to smooth
    # the results further.
    return NevanACContext(Gáµ¥, grid, mesh, Î¦, ğ�’œ, â„‹, ğ�‘�ğ�‘�, 1, 1)
end

"""
    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
"""
function run(nac::NevanACContext)
    hardy = get_n("hardy")
    #
    if hardy
        println("Activate Hardy basis optimization")

        # Determine the minimal Hardy order (`hmin`), update `â„‹` and `ğ�‘�ğ�‘�`.
        calc_hmin!(nac)

        # Determine the optimal Hardy order (`hopt`), update `â„‹` and `ğ�‘�ğ�‘�`.
        calc_hopt!(nac)
    end
end

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

    # Calculate full response function on real axis and write them
    # Note that _G is actually ğ�‘�G, so there is a `-` symbol for the
    # return value.
    _G = C64.(calc_green(nac.ğ�’œ, nac.â„‹, nac.ğ�‘�ğ�‘�))
    fwrite && write_complete(nac.mesh, -_G)

    # Calculate and write the spectral function
    Aout = F64.(imag.(_G) ./ Ï€)
    fwrite && write_spectrum(nac.mesh, Aout)

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

    return Aout, -_G
end

#=
### *Service Functions*
=#

#=
*Remarks* :

**Mobius transformation**

```math
\begin{equation}
z \mapsto \frac{z - i}{z + i}.
\end{equation}
```

It maps `z` from ``\overline{\mathcal{C}^{+}}`` to ``\overline{\mathcal{D}}``.
See `calc_mobius()`.

---

**Inverse Mobius transformation**

```math
\begin{equation}
z \mapsto i \frac{1 + z}{1 - z}.
\end{equation}
```

It maps `z` from ``\overline{\mathcal{D}}`` to ``\overline{\mathcal{C}^{+}}``.
See `calc_inv_mobius()`.

---

**Pick matrix**

```math
\begin{equation}
\mathcal{P} =
\left[
    \frac{1-\lambda_i \lambda^*_j}{1-h(Y_i)h(Y_j)^*}
\right]_{i,j}.
\end{equation}
```

Here, ``Y_i`` is the *i*th Matsubara frequency, ``C_i`` is the value of
``\mathcal{NG}`` at ``Y_i``, and ``\lambda_i`` is the value of ``\theta``
at ``Y_i``:

```math
\begin{equation}
\lambda_i = \theta(Y_i) = \frac{C_i - i}{C_i + i},
\end{equation}
```

```math
\begin{equation}
h(Y_i) = \frac{Y_i - i}{Y_i + i}.
\end{equation}
```

See `calc_pick()`.

---

**Î¦ vector**

At first, we have

```math
\begin{equation}
\phi_{\alpha} = \theta_{\alpha}(Y_{\alpha}).
\end{equation}
```

So

```math
\begin{equation}
\phi_1 = \theta_1 (Y_1),
\end{equation}
```

```math
\begin{equation}
\phi_{\beta} =
\frac{-d_{\beta}\theta(Y_{\beta}) + b_{\beta}}{c_{\beta}\theta(Y_{\beta}) - \alpha_{\beta}}
\end{equation}
```

where

```math
\begin{equation}
\begin{pmatrix}
a_{\beta} & b_{\beta} \\
c_{\beta} & d_{\beta}
\end{pmatrix}
= \prod^{\beta-1}_{\alpha=1}
\begin{pmatrix}
\frac{Y_{\beta}-Y_{\alpha}}{Y_{\beta}-Y^*_{\alpha}} & \phi_{\alpha} \\
\phi^*_{\alpha}\frac{Y_{\beta} - Y_{\alpha}}{Y_{\beta} - Y^*_{\alpha}} & 1
\end{pmatrix},
\end{equation}
```

See `calc_phis()`.

---

**Contractive function Î¸(z)**

```math
\begin{equation}
\theta(z)[z; Î¸_{M+1}(z)] =
\frac{a(z)\theta_{M+1}(z) + b(z)}{c(z)\theta_{M+1}(z) + d(z)}
\end{equation}
```

where

```math
\begin{equation}
\begin{pmatrix}
a(z) & b(z) \\
c(z) & d(z)
\end{pmatrix}
= \prod^{M}_{j=1}
\begin{pmatrix}
\frac{z - Y_j}{z - Y^*_j} & \phi_j \\
\phi^*_j \frac{z - Y_j}{z - Y^*_j} & 1
\end{pmatrix}
\end{equation}
```

See `calc_abcd()` and `calc_theta()`.

---

**Hardy basis**

```math
\begin{equation}
f^k(z) = \frac{1}{\sqrt{\pi}(z + i)}
    \left( \frac{z - i}{z + i} \right)^k.
\end{equation}
```

See `calc_hbasis()` and `calc_hmatrix()`.

---

**Expanding ``\theta_{M+1}`` using Hardy basis**

```math
\theta_{M+1} = \sum^{H}_{k=0} \left[a_k f^k(z) + b_k f^k(z)^*\right]
```

See `calc_theta()`.

---

**Smooth norm**

```math
\begin{equation}
F[A_{\theta_{M+1}}(\omega)] =
    \left|
        1 - \int A_{\theta_{M+1}}(\omega) d\omega
    \right|^2 +
    \alpha \int \left[A^{''}_{\theta_{M+1}}(\omega)\right]^2 d\omega,
\end{equation}
```
where the first term enforces proper normalization while the second
term promotes smoothness by minimizing second derivatives. Here Î±
is a regulation parameter.

See `smooth_norm()`.
=#

"""
    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`](@ref).
"""
function precompute(
    grid::AbstractGrid,
    mesh::AbstractMesh,
    Gáµ¥::Vector{APC}
    )
    # Evaluate Ï• and `abcd` matrices
    Î¦ = calc_phis(grid, Gáµ¥)
    ğ�’œ = calc_abcd(grid, mesh, Î¦)

    # Allocate memory for evaluating Î¸
    # The initial Hardy order is just 1.
    â„‹ = calc_hmatrix(mesh, 1)
    ğ�‘�ğ�‘� = zeros(C64, 2)

    return Î¦, ğ�’œ, â„‹, ğ�‘�ğ�‘�
end

"""
    calc_mobius(z::Vector{APC})

A direct Mobius transformation.

### Arguments
* z -> Complex vector.

### Returns
* val -> Ï†(z), Mobius transformation of z.
"""
function calc_mobius(z::Vector{APC})
    return @. (z - im) / (z + im)
end

"""
    calc_inv_mobius(z::Vector{APC})

An inverse Mobius transformation.

### Arguments
* z -> Complex vector.

### Returns
* val -> Ï†â�»Â¹(z), inverse Mobius transformation of z.
"""
function calc_inv_mobius(z::Vector{APC})
    return @. im * (one(APC) + z) / (one(APC) - z)
end

"""
    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.
"""
function calc_pick(k::I64, â„�::Vector{APC}, Î»::Vector{APC})
    pick = zeros(APC, k, k)

    # Calculate the Pick matrix
    for j = 1:k
        for i = 1:k
            num = one(APC) - Î»[i] * conj(Î»[j])
            den = one(APC) - â„�[i] * conj(â„�[j])
            pick[i,j] = num / den
        end
        pick[j,j] += APC(1e-250)
    end

    # Cholesky decomposition
    return issuccess(cholesky(pick, check = false))
end

"""
    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.
"""
function calc_phis(grid::AbstractGrid, Gáµ¥::Vector{APC})
    ngrid = length(grid)

    # Allocate memory
    Î¦ = zeros(APC, ngrid)
    ğ�’œ = zeros(APC, 2, 2, ngrid)
    âˆ� = zeros(APC, 2, 2)
    ğ�‘” = grid.Ï‰ * im

    # Initialize the `abcd` matrix
    for i = 1:ngrid
        ğ�’œ[:,:,i] .= Matrix{APC}(I, 2, 2)
    end

    # Evaluate Î¦ using recursive algorithm
    Î¦[1] = Gáµ¥[1]
    for j = 1:ngrid-1
        for k = j+1:ngrid
            âˆ�[1,1] = ( ğ�‘”[k] - ğ�‘”[j] ) / ( ğ�‘”[k] - conj(ğ�‘”[j]) )
            âˆ�[1,2] = Î¦[j]
            âˆ�[2,1] = conj(Î¦[j]) * âˆ�[1,1]
            âˆ�[2,2] = one(APC)
            view(ğ�’œ,:,:,k) .= view(ğ�’œ,:,:,k) * âˆ�
        end
        num = ğ�’œ[1,2,j+1] - ğ�’œ[2,2,j+1] * Gáµ¥[j+1]
        den = ğ�’œ[2,1,j+1] * Gáµ¥[j+1] - ğ�’œ[1,1,j+1]
        Î¦[j+1] = num / den
    end

    return Î¦
end

"""
    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.
"""
function calc_abcd(grid::AbstractGrid, mesh::AbstractMesh, Î¦::Vector{APC})
    eta::APF = get_n("eta")

    ngrid = length(grid)
    nmesh = length(mesh)

    ğ�‘” = grid.Ï‰ * im
    ğ�‘š = mesh.mesh .+ im * eta

    ğ�’œ = zeros(APC, 2, 2, nmesh)
    âˆ� = zeros(APC, 2, 2)

    for i = 1:nmesh
        result = Matrix{APC}(I, 2, 2)
        ğ�‘§ = ğ�‘š[i]
        for j = 1:ngrid
            âˆ�[1,1] = ( ğ�‘§ - ğ�‘”[j] ) / ( ğ�‘§ - conj(ğ�‘”[j]) )
            âˆ�[1,2] = Î¦[j]
            âˆ�[2,1] = conj(Î¦[j]) * âˆ�[1,1]
            âˆ�[2,2] = one(APC)
            result *= âˆ�
        end

        ğ�’œ[:,:,i] .= result
    end

    return ğ�’œ
end

"""
    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.
"""
function calc_hbasis(z::APC, k::I64)
    w = ( z - im ) / ( z + im )
    return 1.0 / ( sqrt(Ï€) * (z + im) ) * w^k
end

"""
    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.
"""
function calc_hmatrix(mesh::AbstractMesh, H::I64)
    # Build real axis
    eta::APF = get_n("eta")
    ğ�‘š = mesh.mesh .+ eta * im

    # Allocate memory for the Hardy matrix
    nmesh = length(mesh)
    â„‹ = zeros(APC, nmesh, 2*H)

    # Build the Hardy matrix
    for k = 1:H
        â„‹[:,2*k-1] .= calc_hbasis.(ğ�‘š,k-1)
        â„‹[:,2*k]   .= conj(â„‹[:,2*k-1])
    end

    return â„‹
end

"""
    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.
"""
function calc_theta(ğ�’œ::Array{APC,3}, â„‹::Array{APC,2}, ğ�‘�ğ�‘�::Vector{C64})
    # Well, we should calculate Î¸â‚˜â‚Šâ‚� at first.
    Î¸â‚˜â‚Šâ‚� = â„‹ * ğ�‘�ğ�‘�

    # Then we evaluate Î¸ according Eq. (7)
    num = ğ�’œ[1,1,:] .* Î¸â‚˜â‚Šâ‚� .+ ğ�’œ[1,2,:]
    den = ğ�’œ[2,1,:] .* Î¸â‚˜â‚Šâ‚� .+ ğ�’œ[2,2,:]
    Î¸ = num ./ den

    return Î¸
end

"""
    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(Ï‰).
"""
function calc_green(ğ�’œ::Array{APC,3}, â„‹::Array{APC,2}, ğ�‘�ğ�‘�::Vector{C64})
    Î¸ = calc_theta(ğ�’œ, â„‹, ğ�‘�ğ�‘�)
    gout = calc_inv_mobius(Î¸)

    return gout
end

"""
    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.
"""
function calc_noptim(Ï‰â‚™::Vector{APC}, Gâ‚™::Vector{APC})
    # Get size of input data
    ngrid = length(Ï‰â‚™)

    # Check whether the Pick criterion is applied
    pick = get_n("pick")
    if !pick
        return ngrid
    end

    # Apply invertible Mobius transformation. We actually work at
    # the \bar{ğ�’Ÿ} space.
    ğ�“� = calc_mobius(Ï‰â‚™)
    ğ�’¢ = calc_mobius(-Gâ‚™)

    # Find the optimal value of k until the Pick criterion is violated
    k = 0
    success = true
    while success && k â‰¤ ngrid
        k += 1
        success = calc_pick(k, ğ�“�, ğ�’¢)
    end

    # Return the optimal value for the size of input data
    if !success
        println("The size of input data is optimized to $(k-1)")
        return k - 1
    else
        println("The size of input data is optimized to $(ngrid)")
        return ngrid
    end
end

"""
    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
"""
function calc_hmin!(nac::NevanACContext)
    hmax = get_n("hmax")

    h = 1
    while h â‰¤ hmax
        println("H (Order of Hardy basis) -> $h")

        # Prepare initial â„‹ and ğ�‘�ğ�‘�
        â„‹ = calc_hmatrix(nac.mesh, h)
        ğ�‘�ğ�‘� = zeros(C64, 2*h)

        # Hardy basis optimization
        causality, optim = hardy_optimize!(nac, â„‹, ğ�‘�ğ�‘�, h)

        # Check whether the causality is preserved and the
        # optimization is successful.
        if causality && optim
            nac.hmin = h
            break
        else
            h = h + 1
        end
    end
end

"""
    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
"""
function calc_hopt!(nac::NevanACContext)
    hmax = get_n("hmax")

    for h = nac.hmin + 1:hmax
        println("H (Order of Hardy basis) -> $h")

        # Prepare initial â„‹ and ğ�‘�ğ�‘�
        â„‹ = calc_hmatrix(nac.mesh, h)
        ğ�‘�ğ�‘�  = copy(nac.ğ�‘�ğ�‘�)
        push!(ğ�‘�ğ�‘�, zero(C64))
        push!(ğ�‘�ğ�‘�, zero(C64))
        @assert size(â„‹)[2] == length(ğ�‘�ğ�‘�)

        # Hardy basis optimization
        causality, optim = hardy_optimize!(nac, â„‹, ğ�‘�ğ�‘�, h)

        # Check whether the causality is preserved and the
        # optimization is successful.
        if !(causality && optim)
            break
        end
    end
end

"""
    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.
"""
function hardy_optimize!(
    nac::NevanACContext,
    â„‹::Array{APC,2},
    ğ�‘�ğ�‘�::Vector{C64},
    H::I64
    )
    # Function call to the smooth norm.
    function ğ�‘“(x::Vector{C64})
        return smooth_norm(nac, â„‹, x)
    end

    # Function call to the gradient of the smooth norm.
    #
    # Here we adopt the Zygote package, which implements an automatic
    # differentiation algorithm, to evaluate the gradient of the smooth
    # norm. Of course, we can turn to the finite difference algorithm,
    # which is less efficient.
    function ğ��½!(J::Vector{C64}, x::Vector{C64})
        #J .= Zygote.gradient(ğ�‘“, x)[1]

        # Finite difference algorithm
        J .= gradient_via_fd(ğ�‘“, x)
    end

    # Perform numerical optimization by the BFGS algorithm.
    # If it is failed, please turn to the Optim.jl package.
    # A simplified version is implemented in math.jl.
    res = optimize(ğ�‘“, ğ��½!, ğ�‘�ğ�‘�, max_iter = 500)

    # Check whether the BFGS algorithm is converged
    if !converged(res)
        @info("Sorry, faild to optimize the smooth norm!")
    end

    # Check causality of the solution
    causality = check_causality(â„‹, res.minimizer)

    # Update â„‹ and the corresponding ğ�‘�ğ�‘�
    if causality && (converged(res))
        nac.hopt = H
        nac.ğ�‘�ğ�‘� = res.minimizer
        nac.â„‹ = â„‹
    end

    return causality, converged(res)
end

"""
    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.
"""
function smooth_norm(nac::NevanACContext, â„‹::Array{APC,2}, ğ�‘�ğ�‘�::Vector{C64})
    # Get regulation parameter
    Î± = get_n("alpha")

    # Generate output spectrum
    _G = calc_green(nac.ğ�’œ, â„‹, ğ�‘�ğ�‘�)
    A = F64.(imag.(_G) ./ Ï€)

    # Normalization term
    ğ�‘“â‚� = trapz(nac.mesh, A)

    # Smoothness term
    sd = second_derivative(nac.mesh.mesh, A)
    x_sd = nac.mesh.mesh[2:end-1]
    ğ�‘“â‚‚ = trapz(x_sd, abs.(sd) .^ 2)

    # Assemble the final smooth norm
    ğ��¹ = abs(1.0 - ğ�‘“â‚�)^2 + Î± * ğ�‘“â‚‚

    return F64(ğ��¹)
end

"""
    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
"""
function check_pick(wn::Vector{APC}, gw::Vector{APC}, Nopt::I64)
    freq = calc_mobius(wn[1:Nopt])
    val = calc_mobius(-gw[1:Nopt])

    success = calc_pick(Nopt, val, freq)
    #
    if success
        println("Pick matrix is positive semi-definite.")
    else
        println("Pick matrix is non positive semi-definite matrix in Schur method.")
    end
end

"""
    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.
"""
function check_causality(â„‹::Array{APC,2}, ğ�‘�ğ�‘�::Vector{C64})
    Î¸â‚˜â‚Šâ‚� = â„‹ * ğ�‘�ğ�‘�

    max_theta = findmax(abs.(Î¸â‚˜â‚Šâ‚�))[1]

    if max_theta <= 1.0
        println("max_theta = ", max_theta)
        println("Hardy optimization was success.")
        causality = true
    else
        println("max_theta = ", max_theta)
        println("Hardy optimization was failure.")
        causality = false
    end

    return causality
end
