text:
using LinearAlgebra
using Printf
using Parameters
using Profile
using LinearMaps
using BenchmarkTools
using OrderedCollections
using PyCall

struct FCIProblem
    no::Int  # number of orbitals
    na::Int  # number of alpha
    nb::Int  # number of beta
    dima::Int
    dimb::Int
    dim::Int
    converged::Bool
    restarted::Bool
    iteration::Int
    algorithm::String   #  options: direct/davidson
    n_roots::Int
end

function FCIProblem(no, na, nb)
    na <= no || throw(DimensionMismatch)
    nb <= no || throw(DimensionMismatch)
    dima = calc_nchk(no,na)
    dimb = calc_nchk(no,nb)
    return FCIProblem(no, na, nb, dima, dimb, dima*dimb, false, false, 1, "direct", 1)
end

function display(p::FCIProblem)
    @printf(" FCIProblem:: #Orbs = %-3i #α = %-2i #β = %-2i Dimension: %-9i\n",p.no,p.na,p.nb,p.dim)
    #@printf(" FCIProblem::  NOrbs: %2i NAlpha: %2i NBeta: %2i Dimension: %-9i\n",p.no,p.na,p.nb,p.dim)
end

function compute_spin_diag_terms_full!(H, P::FCIProblem, Hmat)
    #={{{=#

    print(" Compute same spin terms.\n")
    size(Hmat,1) == P.dim || throw(DimensionMismatch())

    Hdiag_a = FCI.precompute_spin_diag_terms(H,P,P.na)
    Hdiag_b = FCI.precompute_spin_diag_terms(H,P,P.nb)
    Hmat .+= kron(Matrix(1.0I, P.dimb, P.dimb), Hdiag_a)
    Hmat .+= kron(Hdiag_b, Matrix(1.0I, P.dima, P.dima))

end
#=}}}=#

"""
    build_H_matrix(ints, P::FCIProblem)

Build the Hamiltonian defined by `ints` in the Slater Determinant Basis
in the sector of Fock space specified by `P`
"""
function build_H_matrix(ints, P::FCIProblem)

    Hmat = zeros(P.dim, P.dim)

    Hdiag_a = precompute_spin_diag_terms(ints,P,P.na)
    Hdiag_b = precompute_spin_diag_terms(ints,P,P.nb)
    #
    #   Create ci_strings
    ket_a = DeterminantString(P.no, P.na)
    ket_b = DeterminantString(P.no, P.nb)
    bra_a = DeterminantString(P.no, P.na)
    bra_b = DeterminantString(P.no, P.nb)
    #
    #   Add spin diagonal components
    Hmat += kron(Matrix(1.0I, P.dimb, P.dimb), Hdiag_a)
    Hmat += kron(Hdiag_b, Matrix(1.0I, P.dima, P.dima))
    #
    #   Add opposite spin term (todo: make this reasonably efficient)
    Hmat += compute_ab_terms_full(ints, P)

    Hmat = .5*(Hmat+Hmat')

    return Hmat
end

"""
    compute_fock_diagonal!(H, P::FCIProblem, e_mf::Float64)
"""
function compute_fock_diagonal(P::FCIProblem, orb_energies::Vector{Float64}, e_mf::Float64)
    ket_a = DeterminantString(P.no, P.na)
    ket_b = DeterminantString(P.no, P.nb)

    meanfield_e = e_mf
    meanfield_e -=sum(orb_energies[ket_a.config])
    meanfield_e -=sum(orb_energies[ket_b.config])

    fdiag = zeros(ket_a.max, ket_b.max)

    reset!(ket_a)
    reset!(ket_b)
    for Ib in 1:ket_b.max
        eb = sum(orb_energies[ket_b.config]) + meanfield_e
        for Ia in 1:ket_a.max
            fdiag[Ia,Ib] = eb + sum(orb_energies[ket_a.config])
            incr!(ket_a)
        end
        incr!(ket_b)
    end
    fdiag = reshape(fdiag,ket_a.max*ket_b.max)
    #display(fdiag)
    return fdiag
end

"""
    compute_ab_terms_full!(H, P::FCIProblem, Hmat)
"""
function compute_ab_terms_full!(H, P::FCIProblem, Hmat)
    #={{{=#

    print(" Compute opposite spin terms.\n")
    @assert(size(Hmat,1) == P.dim)

    #v = transpose(vin)

    #   Create local references to ci_strings
    ket_a = DeterminantString(P.no, P.na)
    ket_b = DeterminantString(P.no, P.nb)
    bra_a = DeterminantString(P.no, P.na)
    bra_b = DeterminantString(P.no, P.nb)

    ket_a_lookup = fill_ca_lookup(ket_a)
    ket_b_lookup = fill_ca_lookup(ket_b)

    reset!(ket_b)

    for Kb in 1:ket_b.max

        reset!(ket_a)
        for Ka in 1:ket_a.max
            K = Ka + (Kb-1) * ket_a.max

            #  <pq|rs> p'q'sr  --> (pr|qs) (a,b)
            for r in 1:ket_a.no
                for p in 1:ket_a.no
                    sign_a, La = ket_a_lookup[Ka][p+(r-1)*ket_a.no]
                    if La == 0
                        continue
                    end

                    Lb = 1
                    sign_b = 1
                    L = 1
                    for s in 1:ket_b.no
                        for q in 1:ket_b.no
                            sign_b, Lb = ket_b_lookup[Kb][q+(s-1)*ket_b.no]

                            if Lb == 0
                                continue
                            end

                            L = La + (Lb-1) * bra_a.max

                            Hmat[K,L] += H.h2[p,r,q,s] * sign_a * sign_b

                        end
                    end
                end
            end
            incr!(ket_a)

        end
        incr!(ket_b)
    end
    return
end
#=}}}=#

"""
    compute_ab_terms_full(H, P::FCIProblem)
"""
function compute_ab_terms_full(H, P::FCIProblem)
    #={{{=#

    #print(" Compute opposite spin terms. Shape of v: ", size(v), "\n")

    #v = transpose(vin)

    Hmat = zeros(Float64, P.dim, P.dim)

    #   Create local references to ci_strings
    ket_a = DeterminantString(P.no, P.na)
    ket_b = DeterminantString(P.no, P.nb)
    bra_a = DeterminantString(P.no, P.na)
    bra_b = DeterminantString(P.no, P.nb)

    ket_a_lookup = fill_ca_lookup(ket_a)
    ket_b_lookup = fill_ca_lookup(ket_b)
    ket_a_lookup2 = fill_ca_lookup2(ket_a)
    ket_b_lookup2 = fill_ca_lookup2(ket_b)

    a_max = bra_a.max
    reset!(ket_b)

    for Kb in 1:ket_b.max

        reset!(ket_a)
        for Ka in 1:ket_a.max
            K = Ka + (Kb-1) * ket_a.max

            #  <pq|rs> p'q'sr  --> (pr|qs) (a,b)
            for r in 1:ket_a.no
                for p in 1:ket_a.no
                    #sign_a, La = ket_a_lookup[Ka][p+(r-1)*ket_a.no]
                    La = ket_a_lookup2[p,r,Ka]
                    if La == 0
                        continue
                    end
                    sign_a = sign(La)
                    La = abs(La)

                    Lb = 1
                    sign_b = 1
                    L = 1
                    for s in 1:ket_b.no
                        for q in 1:ket_b.no
                            Lb = ket_b_lookup2[q,s,Kb]
                            if Lb == 0
                                continue
                            end
                            sign_b = sign(Lb)
                            Lb = abs(Lb)
                            #sign_b, Lb = ket_b_lookup[Kb][q+(s-1)*ket_b.no]

                            if Lb == 0
                                continue
                            end

                            L = La + (Lb-1) * a_max

                            Hmat[K,L] += H.h2[p,r,q,s] * sign_a * sign_b
                            continue
                        end
                    end
                end
            end
            incr!(ket_a)

        end
        incr!(ket_b)
    end
    #sig = transpose(sig)
    return Hmat
end
#=}}}=#

# Helper functions for Olsen's agorithm
function _gather!(FJb::Vector{Float64}, occ::Vector{Int}, vir::Vector{Int}, vkl::Array{Float64,2}, Ib::Int,ket_b_lookup)
#={{{=#
    i::Int = 1
    j::Int = 1
    Jb::Int = 1
    sgn::Float64 = 1.0
    @inbounds @simd for j in occ
        for i in vir
            Jb = ket_b_lookup[i,j,Ib]
            sgn = sign(Jb)
            Jb = abs(Jb)
            FJb[Jb] = FJb[Jb] + vkl[j,i]*sgn
        end
    end
end
#=}}}=#

function _mult!(Ckl::Array{Float64,3}, FJb::Array{Float64,1}, VI::Array{Float64,2})
    #={{{=#
    VI .= 0
    nI = size(Ckl)[1]
    n_roots::Int = size(Ckl)[3]
    ket_max = size(FJb)[1]
    tmp = 0.0
    for si in 1:n_roots
        @views V = VI[:,si]
        for Jb in 1:ket_max
            tmp = FJb[Jb]
            if abs(tmp) > 1e-14
                @inbounds @simd for I in 1:nI
                    VI[I,si] += tmp*Ckl[I,Jb,si]
                end
                #@views LinearAlgebra.axpy!(tmp, Ckl[:,Jb,si], VI[:,si])
                #@inbounds VI[:,si] .+= tmp .* Ckl[:,Jb,si]
                #@inbounds @views @. VI[:,si] += tmp * Ckl[:,Jb,si]
            end
        end
    end
end
#=}}}=#

function _scatter!(sig::Array{Float64,3}, VI::Array{Float64,2}, L::Vector{Int}, R::Vector{Int}, Ib::Int)
    #={{{=#
    n_roots = size(sig)[3]
    @inbounds @simd for si in 1:n_roots
        for Li in 1:length(L)
            sig[R[Li],Ib,si] += VI[Li,si]
            #@inbounds sig[R[Li],Ib,si] += VI[Li,si]
        end
    end
end
#=}}}=#

function _getCkl!(Ckl::Array{Float64,3}, v,L::Vector{Int})
    #={{{=#
    nI = length(L)
    n_roots = size(v)[3]
    ket_max = size(v)[2]
    @inbounds @simd for si in 1:n_roots
        for Jb in 1:ket_max
            for Li in 1:nI
                Ckl[Li,Jb,si] = v[abs(L[Li]), Jb, si] * sign(L[Li])
            end
        end
    end
end
#=}}}=#

"""
    compute_ab_terms2(v, H, P::FCIProblem,
                          ket_a_lookup, ket_b_lookup)
"""
function compute_ab_terms2(v, H, P::FCIProblem,
                          ket_a_lookup, ket_b_lookup)
    #={{{=#

    #print(" Compute opposite spin terms. Shape of v: ", size(v), "\n")
    @assert(size(v,1)*size(v,2) == P.dim)

    #   Create local references to ci_strings
    ket_a = DeterminantString(P.no, P.na)
    ket_b = DeterminantString(P.no, P.nb)
    bra_a = DeterminantString(P.no, P.na)
    bra_b = DeterminantString(P.no, P.nb)

    no = ket_a.no

    a_max::Int = bra_a.max
    reset!(ket_b)

    #
    #   sig3(Ia,Ib,s) = <Ia|k'l|Ja> <Ib|i'j|Jb> V(ij,kl) C(Ja,Jb,s)
    n_roots::Int = size(v,3)
    #v = reshape(v,ket_a.max, ket_b.max, n_roots)
    sig = zeros(Float64, ket_a.max, ket_b.max, n_roots)
    FJb_scr1 = zeros(Float64, ket_b.max)
    Ckl_scr1 = zeros(Float64, binomial(ket_a.no-1,ket_a.ne-1), size(v)[2], size(v)[3])
    Ckl_scr2 = zeros(Float64, binomial(ket_a.no-2,ket_a.ne-1), size(v)[2], size(v)[3])
    Ckl = Array{Float64,3}
    virt = zeros(Int,ket_b.no-ket_b.ne)
    diff_ref = Set(collect(1:ket_b.no))
    FJb = copy(FJb_scr1)
    for k in 1:ket_a.no,  l in 1:ket_a.no
        #@printf(" %4i, %4i\n",k,l)
        L = Vector{Int}()
        R = Vector{Int}()
        #sgn = Vector{Int}()
        #for (Iidx,I) in enumerate(ket_a_lookup[k,l,:])
        #    if I[2] != 0
        #        push!(R,Iidx)
        #        push!(sgn,I[1])
        #        push!(L,I[2])
        #    end
        #end
        for (Iidx,I) in enumerate(ket_a_lookup[k,l,:])
            if I != 0
                push!(R,Iidx)
                push!(L,I)
            end
        end
        VI = zeros(Float64, length(L),n_roots)
        #Ckl = zeros(Float64, size(v)[2], length(L), size(v)[3])
        if k==l
            Ckl = deepcopy(Ckl_scr1)
        else
            Ckl = deepcopy(Ckl_scr2)
        end
        #nI = length(L)
        #for si in 1:n_roots
        #    @simd for Jb in 1:ket_b.max
        #        for Li in 1:nI
        #            @inbounds Ckl[Li,Jb,si] = v[abs(L[Li]), Jb, si] * sign(L[Li])
        #        end
        #    end
        #end
        _getCkl!(Ckl, v, L)

        vkl = H.h2[:,:,l,k]
        reset!(ket_b)
        for Ib in 1:ket_b.max
            fill!(FJb,0.0)
            Jb = 1
            sgn = 1
            zero_num = 0

            no = ket_b.no
            ne = ket_b.ne
            nv = no-ne
            scr1 = 0.0
            get_unoccupied!(virt, ket_b)
            #virt = setdiff(diff_ref, ket_b.config)
            #@btime $_gather!($FJb, $ket_b.config, $virt, $vkl, $Ib, $ket_b_lookup)
            _gather!(FJb, ket_b.config, virt, vkl, Ib, ket_b_lookup)
            #
            # diagonal part
            for j in ket_b.config
                FJb[Ib] += H.h2[j,j,l,k]
            end

            if 1==0
                _mult_old!(Ckl, FJb, VI)
            end
            if 1==0
                @tensor begin
                    VI[I,s] = FJb[J] * Ckl[J,I,s]
                end
            end
            if 1==1
                _mult!(Ckl, FJb, VI)
            end

            _scatter!(sig,VI,L,R,Ib)
            #@btime $_scatter!($sig,$VI,$L,$R,$Ib)

            incr!(ket_b)
        end
    end

    return sig

end
#=}}}=#

function _ss_sum!(sig::Array{Float64,3}, v::Array{Float64,3}, F::Vector{Float64},Ia::Int)
    nKb     = size(v)[1]
    n_roots = size(v)[2]
    nJa     = size(v)[3]

    for Ja in 1:nJa
        if abs(F[Ja]) > 1e-14
            @inbounds @simd for si in 1:n_roots
                for Kb in 1:nKb
                    sig[Kb,si,Ia] += F[Ja]*v[Kb,si,Ja]
                end
            end
        end
    end
end

function _ss_sum_Ia!(sig::Array{Float64,3}, v::Array{Float64,3}, F::Vector{Float64},Ia::Int)
    nJa = size(v)[3]
    nKb = size(v)[1]
    n_roots = size(v)[2]

    for Ja in 1:nJa
        if abs(F[Ja]) > 1e-14
            @inbounds @simd for si in 1:n_roots
                for Kb in 1:nKb
                    sig[Kb,si,Ia] += F[Ja]*v[Kb,si,Ja]
                end
            end
        end
    end
end

"""
    compute_aa_terms2(v, H, P::FCIProblem,
                          ket_a_lookup)
"""
function compute_ss_terms2(v, H, P::FCIProblem, ket_a_lookup, ket_b_lookup)
    #={{{=#

    #print(" Compute opposite spin terms. Shape of v: ", size(v), "\n")
    @assert(size(v,1)*size(v,2) == P.dim)

    #v = transpose(vin)

    #   Create local references to ci_strings
    ket_a = DeterminantString(P.no, P.na)
    ket_b = DeterminantString(P.no, P.nb)
    bra_a = DeterminantString(P.no, P.na)
    bra_b = DeterminantString(P.no, P.nb)

    no = ket_a.no
    n_roots = P.n_roots

    a_max::Int = bra_a.max
    reset!(ket_a)

    #
    #   sig1(Ia,Ib,s) = <Ib|i'j|Jb> V(ij,kl) C(Ja,Jb,s)
    n_roots::Int = size(v,3)
    #v = reshape(v,ket_a.max, ket_b.max, n_roots)
    sig = zeros(Float64, ket_a.max, ket_b.max, n_roots)
    size(sig) == size(v) || throw(DimensionError())

    h1eff = deepcopy(H.h1)
    @tensor begin
        h1eff[p,q] -= .5 * H.h2[p,j,j,q]
    end

    #bb
    sig = permutedims(sig,[1,3,2])
    v = permutedims(v,[1,3,2])

    ket = ket_b
    reset!(ket)
    F = zeros(ket_b.max)
    for I in 1:ket.max
        F .= 0
        for k in 1:ket.no, l in 1:ket.no
            K = ket_b_lookup[k,l,I]
            if K == 0
                continue
            end
            sign_kl = sign(K)
            K = abs(K)

            @inbounds F[K] += sign_kl * h1eff[k,l]
            for i in 1:ket.no, j in 1:ket.no
                J = ket_b_lookup[i,j,K]
                if J == 0
                    continue
                end
                sign_ij = sign(J)
                J = abs(J)
                if sign_kl == sign_ij
                    @inbounds F[J] += .5 * H.h2[i,j,k,l]
                else
                    @inbounds F[J] -= .5 * H.h2[i,j,k,l]
                end
            end
        end
        _ss_sum!(sig,v,F,I)
    end
    #sig = permutedims(sig,[1,3,2])
    #v = permutedims(v,[1,3,2])

    #aa
    sig = permutedims(sig,[3,2,1])
    v = permutedims(v,[3,2,1])

    ket = ket_a
    reset!(ket)
    F = zeros(ket_a.max)
    bra = deepcopy(ket)
    for I in 1:ket.max
        F .= 0
        for k in 1:ket.no, l in 1:ket.no
            K = ket_a_lookup[k,l,I]
            if K == 0
                continue
            end
            sign_kl = sign(K)
            K = abs(K)

            @inbounds F[K] += sign_kl * h1eff[k,l]
            for i in 1:ket.no, j in 1:ket.no
                J = ket_a_lookup[i,j,K]
                if J == 0
                    continue
                end
                sign_ij = sign(J)
                J = abs(J)
                if sign_kl == sign_ij
                    @inbounds F[J] += .5 * H.h2[i,j,k,l]
                else
                    @inbounds F[J] -= .5 * H.h2[i,j,k,l]
                end

            end
        end

        _ss_sum_Ia!(sig,v,F,I)

    end

    sig = permutedims(sig,[3,1,2])
    v = permutedims(v,[3,1,2])

    return sig

end
#=}}}=#

"""
    compute_ab_terms(v, H, P::FCIProblem)
"""
function compute_ab_terms(v, H, P::FCIProblem)
    #={{{=#

    #print(" Compute opposite spin terms. Shape of v: ", size(v), "\n")
    @assert(size(v,1) == P.dim)

    #v = transpose(vin)

    sig = deepcopy(v)
    sig .= 0

    #   Create local references to ci_strings
    ket_a = DeterminantString(P.no, P.na)
    ket_b = DeterminantString(P.no, P.nb)
    bra_a = DeterminantString(P.no, P.na)
    bra_b = DeterminantString(P.no, P.nb)

    ket_a_lookup = fill_ca_lookup2(ket_a)
    ket_b_lookup = fill_ca_lookup2(ket_b)

    a_max = bra_a.max
    reset!(ket_b)

    n_roots = size(sig,2)
    scr = zeros(1,ket_a.max*ket_b.max)
    for Kb in 1:ket_b.max

        reset!(ket_a)
        for Ka in 1:ket_a.max
            K = Ka + (Kb-1) * ket_a.max

            #  <pq|rs> p'q'sr  --> (pr|qs) (a,b)
            for r in 1:ket_a.no
                for p in 1:ket_a.no
                    #sign_a, La = ket_a_lookup[Ka][p+(r-1)*ket_a.no]
                    La = ket_a_lookup[p,r,Ka]
                    sign_a = sign(La)
                    La = abs(La)
                    if La == 0
                        continue
                    end

                    Lb = 1
                    sign_b = 1
                    L = 1
                    for s in 1:ket_b.no
                        for q in 1:ket_b.no
                            Lb = ket_b_lookup[q,s,Kb]
                            sign_b = sign(Lb)
                            Lb = abs(Lb)

                            if Lb == 0
                                continue
                            end

                            L = La + (Lb-1) * a_max

                            #sig[K,:] += H.h2[p,r,q,s] * v[L,:]
                            #sig[K,:] .+= H.h2[p,r,q,s] * sign_a * sign_b * v[L,:]
                            for si in 1:n_roots
                                sig[K,si] += H.h2[p,r,q,s] * sign_a * sign_b * v[L,si]
                                #@views sig[K,si] .+= H.h2[p,r,q,s] * sign_a * sign_b * v[L,si]
                            end
                            continue
                        end
                    end
                end
            end
            incr!(ket_a)

        end
        incr!(ket_b)
    end
    #sig = transpose(sig)
    return sig
end
#=}}}=#

"""
    compute_ab_terms(v, H, P::FCIProblem, ket_a_lookup, ket_b_lookup)
"""
function compute_ab_terms(v, H, P::FCIProblem, ket_a_lookup, ket_b_lookup)
    #={{{=#

    #print(" Compute opposite spin terms. Shape of v: ", size(v), "\n")
    size(v,1) == P.dim || throw(DimensionError())

    #v = transpose(vin)

    sig = 0*v

    #   Create local references to ci_strings
    ket_a = DeterminantString(P.no, P.na)
    ket_b = DeterminantString(P.no, P.nb)
    bra_a = DeterminantString(P.no, P.na)
    bra_b = DeterminantString(P.no, P.nb)

    a_max = bra_a.max
    reset!(ket_b)

    n_roots = size(sig,2)
    scr = zeros(1,ket_a.max*ket_b.max)
    for Kb in 1:ket_b.max

        reset!(ket_a)
        for Ka in 1:ket_a.max
            K = Ka + (Kb-1) * ket_a.max

            #  <pq|rs> p'q'sr  --> (pr|qs) (a,b)
            for r in 1:ket_a.no
                for p in 1:ket_a.no
                    sign_a, La = ket_a_lookup[Ka][p+(r-1)*ket_a.no]
                    if La == 0
                        continue
                    end

                    Lb = 1
                    sign_b = 1
                    L = 1
                    for s in 1:ket_b.no
                        for q in 1:ket_b.no
                            sign_b, Lb = ket_b_lookup[Kb][q+(s-1)*ket_b.no]

                            if Lb == 0
                                continue
                            end

                            L = La + (Lb-1) * a_max

                            #sig[K,:] += H.h2[p,r,q,s] * v[L,:]
                            #sig[K,:] += H.h2[p,r,q,s] * sign_a * sign_b * v[L,:]
                            for si in 1:n_roots
                                sig[K,si] += H.h2[p,r,q,s] * sign_a * sign_b * v[L,si]
                                #@views sig[K,si] .+= H.h2[p,r,q,s] * sign_a * sign_b * v[L,si]
                            end
                            continue
                        end
                    end
                end
            end
            incr!(ket_a)

        end
        incr!(ket_b)
    end
    #sig = transpose(sig)
    return sig
end
#=}}}=#

"""
    precompute_spin_diag_terms(H, P::FCIProblem, e)
"""
function precompute_spin_diag_terms(H, P::FCIProblem, e)
    #={{{=#

    #   Create local references to ci_strings
    ket = DeterminantString(P.no, e)
    bra = DeterminantString(P.no, e)

    ket_ca_lookup = fill_ca_lookup(ket)

    Hout = zeros(ket.max,ket.max)

    reset!(ket)

    for K in 1:ket.max

        #  hpq p'q
        for p in 1:ket.no
            for q in 1:ket.no
                bra = deepcopy(ket)
                apply_annihilation!(bra,q)
                if bra.sign == 0
                    continue
                end
                apply_creation!(bra,p)
                if bra.sign == 0
                    continue
                end

                L = calc_linear_index(bra)

                term = H.h1[q,p]
                Hout[K,L] += term * bra.sign
            end
        end

        #  <pq|rs> p'q'sr -> (pr|qs)
        for r in 1:ket.no
            for s in r+1:ket.no
                for p in 1:ket.no
                    for q in p+1:ket.no

                        bra = deepcopy(ket)

                        apply_annihilation!(bra,r)
                        if bra.sign == 0
                            continue
                        end
                        apply_annihilation!(bra,s)
                        if bra.sign == 0
                            continue
                        end
                        apply_creation!(bra,q)
                        if bra.sign == 0
                            continue
                        end
                        apply_creation!(bra,p)
                        if bra.sign == 0
                            continue
                        end
                        L = calc_linear_index(bra)
                        Ipqrs = H.h2[p,r,q,s]-H.h2[p,s,q,r]
                        Hout[K,L] += bra.sign*Ipqrs
#                        if bra.sign == -1
#                            Hout[K,L] -= Ipqrs
#                        elseif bra.sign == +1
#                            Hout[K,L] += Ipqrs
#                        else
#                            throw(Exception())
#                        end
                    end
                end
            end
        end
        incr!(ket)
    end
    return Hout
end
#=}}}=#

"""
    get_map(ham, prb::FCIProblem, HdiagA, HdiagB)

Assumes you've already computed the spin diagonal components
"""
function get_map(ham, prb::FCIProblem, HdiagA, HdiagB)
    #=
    Get LinearMap with takes a vector and returns action of H on that vector
    =#
    #={{{=#
    ket_a = DeterminantString(prb.no, prb.na)
    ket_b = DeterminantString(prb.no, prb.nb)

    lookup_a = fill_ca_lookup2(ket_a)
    lookup_b = fill_ca_lookup2(ket_b)
    iters = 0
    function mymatvec(v)
        iters += 1
        #@printf(" Iter: %4i\n", iters)
        nr = 0
        if length(size(v)) == 1
            nr = 1
            v = reshape(v,ket_a.max*ket_b.max, nr)
        else
            nr = size(v)[2]
        end
        v = reshape(v, ket_a.max, ket_b.max, nr)
        sig = compute_ab_terms2(v, ham, prb, lookup_a, lookup_b)
        @tensor begin
            sig[I,J,s] += HdiagA[I,K] * v[K,J,s]
            sig[I,J,s] += HdiagB[J,K] * v[I,K,s]
        end

        v = reshape(v, ket_a.max*ket_b.max, nr)
        sig = reshape(sig, ket_a.max*ket_b.max, nr)
        return sig
    end
    return LinearMap(mymatvec, prb.dim, prb.dim; issymmetric=true, ismutating=false, ishermitian=true)
end
#=}}}=#

"""
    get_map(ham, prb::FCIProblem)

Get LinearMap with takes a vector and returns action of H on that vector
"""
function get_map(ham, prb::FCIProblem)
    #={{{=#
    ket_a = DeterminantString(prb.no, prb.na)
    ket_b = DeterminantString(prb.no, prb.nb)

    #@btime lookup_a = $fill_ca_lookup2($ket_a)
    lookup_a = fill_ca_lookup2(ket_a)
    lookup_b = fill_ca_lookup2(ket_b)
    iters = 0
    function mymatvec(v)
        iters += 1
        #@printf(" Iter: %4i\n", iters)
        nr = 0
        if length(size(v)) == 1
            nr = 1
            v = reshape(v,ket_a.max*ket_b.max, nr)
        else
            nr = size(v)[2]
        end
        v = reshape(v, ket_a.max, ket_b.max, nr)
        sig = compute_ab_terms2(v, ham, prb, lookup_a, lookup_b)
        sig += compute_ss_terms2(v, ham, prb, lookup_a, lookup_b)

        v = reshape(v, ket_a.max*ket_b.max, nr)
        sig = reshape(sig, ket_a.max*ket_b.max, nr)
        return sig
    end
    return LinearMap(mymatvec, prb.dim, prb.dim; issymmetric=true, ismutating=false, ishermitian=true)
end
#=}}}=#

"""
    run_fci(ints, prb::FCIProblem)

input: ints is a struct containing 0, 2, and 4 dimensional tensors
- `h0`: energy shift
- `h1`: 1 electron integrals
- `h2`: 2 electron integrals (chemists notation)
- `prb`: FCIProblem just defines the current CI problem (i.e., fock sector)

ints is simply an InCoreInts object from FermiCG

"""
function run_fci(ints, problem::FCIProblem; v0=nothing, nroots=1, tol=1e-6,
                precompute_ss = false)

    if precompute_ss
        print(" Compute spin_diagonal terms\n")
        @time Hdiag_a = StringCI.precompute_spin_diag_terms(ints,problem,problem.na)
        @time Hdiag_b = StringCI.precompute_spin_diag_terms(ints,problem,problem.nb)
        print(" done\n")

        Hmap = StringCI.get_map(ints, problem, Hdiag_a, Hdiag_b)
    else
        Hmap = StringCI.get_map(ints, problem)
    end

    e = 0
    v = Array{Float64,2}
    if v0 == nothing
        @time e,v = eigs(Hmap, nev = nroots, which=:SR, tol=tol)
        e = real(e)
        for ei in e
            @printf(" Energy: %12.8f\n",ei+ints.h0)
        end
    else
        @time e,v = eigs(Hmap, v0=v0[:,1], nev = nroots, which=:SR, tol=tol)
        e = real(e)
        for ei in e
            @printf(" Energy: %12.8f\n",ei+ints.h0)
        end
    end
    return e,v
end

"""
    build_s2(prb::FCIProblem)
- `prb`: FCIProblem just defines the current CI problem (i.e., fock sector)
"""
function build_S2_matrix(P::FCIProblem)

    #={{{=#
    S2 = zeros(P.dim, P.dim)

    #   Create ci_strings
    ket_a = DeterminantString(P.no, P.na)
    ket_b = DeterminantString(P.no, P.nb)
    bra_a = DeterminantString(P.no, P.na)
    bra_b = DeterminantString(P.no, P.nb)

    #   lookup the ket space
    ket_a_lookup = fill_ca_lookup(ket_a)
    ket_b_lookup = fill_ca_lookup(ket_b)

    reset!(ket_b)
    for Kb in 1:ket_b.max

        reset!(ket_a)
        for Ka in 1:ket_a.max
            K = Ka + (Kb-1) * ket_a.max

	    # Sz.Sz

            for ai in ket_a.config
                for aj in ket_a.config
                    if ai != aj
                        S2[K,K] += 0.25
                    end
                end
            end

            for bi in ket_b.config
                for bj in ket_b.config
                    if bi != bj
                        S2[K,K] += 0.25
                    end
                end
            end

            for ai in ket_a.config
                for bj in ket_b.config
                    if ai != bj
                        S2[K,K] -= 0.50
                    end
                end
            end

	    # Sp.Sm
            for ai in ket_a.config
                if ai in ket_b.config
		    temp = 10
		else
                    S2[K,K] += 0.75
                end
            end

	    # Sm.Sp
            for bi in ket_b.config
                if bi in ket_a.config
		    temp = 10
		else
                    S2[K,K] += 0.75
                end
            end

            ket_a2 = DeterminantString(P.no, P.na)
            ket_b2 = DeterminantString(P.no, P.nb)

	    for ai in ket_a.config
	        for bj in ket_b.config
	            if ai ∉ ket_b.config
	                if bj ∉ ket_a.config

	        	    ket_a2 = deepcopy(ket_a)
	        	    ket_b2 = deepcopy(ket_b)

                            apply_annihilation!(ket_a2,ai)
	        	    ket_a2.sign != 0 || continue

                            apply_creation!(ket_b2,ai)
	        	    ket_b2.sign != 0 || continue

                            apply_creation!(ket_a2,bj)
	        	    ket_a2.sign != 0 || continue

                            apply_annihilation!(ket_b2,bj)
	        	    ket_b2.sign != 0 || continue

	        	    sign_a = ket_a2.sign
	        	    sign_b = ket_b2.sign

                            La = calc_linear_index(ket_a2)
                            Lb = calc_linear_index(ket_b2)

                            L = La + (Lb-1) * ket_a.max
			    #print("Init ",ket_a.config,"    ",ket_b.config,"\n")
			    #print("Final",ket_a2.config,"    ",ket_b2.config,"\n")
    	        	    #print(K,"  ",L,"\n")
                            S2[K,L] += 1 * sign_a * sign_b
	        	end
	            end
	        end
	    end

            incr!(ket_a)

        end
        incr!(ket_b)
    end
    return S2
end
#=}}}=#

"""
    svd_state(prb::FCIProblem)
Do an SVD of the FCI vector partitioned into clusters with (norbs1 | norbs2)
where the orbitals are assumed to be ordered for cluster 1| cluster 2 haveing norbs1 and
norbs2, respectively.

- `prb`: FCIProblem just defines the current CI problem (i.e., fock sector)
- `norbs1`:number of orbitals in left cluster
- `norbs2`:number of orbitals in right cluster
- `svd_thresh`: the threshold below which the states will be discarded
"""
function svd_state(v,P::FCIProblem,norbs1,norbs2,svd_thresh)

    #={{{=#

    @assert(norbs1+norbs2 ==P.no)

    schmidt_basis = OrderedDict()
    #vector = OrderedDict{Tuple{UInt8,UInt8},Float64}()
    vector = OrderedDict{Tuple{Int,Int},Any}()

    #schmidt_basis = Dict{Tuple,Matrix{Float64}}

    println("----------------------------------------")
    println("          SVD of state")
    println("----------------------------------------")

    # Create ci_strings
    ket_a = DeterminantString(P.no, P.na)
    ket_b = DeterminantString(P.no, P.nb)

    v = reshape(v,(ket_a.max, ket_b.max))
    @assert(size(v,1) == ket_a.max)
    @assert(size(v,2) == ket_b.max)

    fock_labels_a = Array{Int,1}(undef,ket_a.max)
    fock_labels_b = Array{Int,1}(undef,ket_b.max)

    # Get the fock space using the bisect method in python
    bisect = pyimport("bisect")
    for I in 1:ket_a.max
        fock_labels_a[I] = bisect.bisect(ket_a.config,norbs1)
        #print("nick: ", ket_a.config, " " , norbs1, " ", fock_labels_a[I], "\n")
        incr!(ket_a)
    end
    for I in 1:ket_b.max
        fock_labels_b[I] = bisect.bisect(ket_b.config,norbs1)
        incr!(ket_b)
    end
    for J in 1:ket_b.max
        for I in 1:ket_a.max
            fock = (fock_labels_a[I], fock_labels_b[J])

            #if fock in vector
            #    append!(vector[fock],v[I,J])
            #else
            #    vector[fock] = [v[I,J]]
            #end
            try
                append!(vector[tuple(fock_labels_a[I],fock_labels_b[J])],v[I,J])
            catch
                vector[tuple(fock_labels_a[I],fock_labels_b[J])] = [v[I,J]]
            end
        end
    end

    for (fock,fvec)  in vector

        println()
        @printf("Prepare Fock Space:  %iα, %iβ\n",fock[1] ,fock[2] )

        ket_a1 = DeterminantString(norbs1, fock[1])
        ket_b1 = DeterminantString(norbs1, fock[2])

        ket_a2 = DeterminantString(norbs2, P.na - fock[1])
        ket_b2 = DeterminantString(norbs2, P.nb - fock[2])

        temp_fvec = reshape(fvec,ket_b1.max*ket_b2.max,ket_a1.max*ket_a2.max)'
        #temp_fvec = reshape(fvec,ket_b1.max*ket_b2.max,ket_a1.max*ket_a2.max)'
        #st = "temp_fvec"*string(fock)
        #npzwrite(st, temp_fvec)

        #when swapping alpha2 and beta1 do we flip sign?
        sign = 1
        if (P.na-fock[1])%2==1 && fock[2]%2==1
            sign = -1
        end
        #println("sign",sign)
        @printf("   Dimensions: %5i x %-5i \n",ket_a1.max*ket_b1.max, ket_a2.max*ket_b2.max)

        norm_curr = fvec' * fvec
        @printf("   Norm: %12.8f\n",sqrt(norm_curr))
        #println(size(fvec))
        #display(fvec)

        fvec = sign *fvec

        #opposite to python with transpose on fvec
        #fvec2 = reshape(fvec',ket_b2.max,ket_b1.max,ket_a2.max,ket_a1.max)
        fvec2 = reshape(fvec,ket_a1.max,ket_a2.max,ket_b1.max,ket_b2.max)
        fvec3 = permutedims(fvec2, [ 1, 3, 2, 4])
        fvec4 = reshape(fvec3,ket_a1.max*ket_b1.max,ket_a2.max*ket_b2.max)

        # fvec4 is transpose of what we have in python code
        fvec5 = fvec4'

        F = svd(fvec5,full=true)

        nkeep = 0
        @printf("   %5s %12s\n","State","Weight")
        for (ni_idx,ni) in enumerate(F.S)
            if ni > svd_thresh
                nkeep += 1
                @printf("   %5i %12.8f\n",ni_idx,ni)
            else
                @printf("   %5i %12.8f (discarded)\n",ni_idx,ni)
            end
        end

        if nkeep > 0
            schmidt_basis[fock] = Matrix(F.U[:,1:nkeep])
            #st = "fin_vec"*string(fock)
            #npzwrite(st, F.U[:,1:nkeep])
        end

        #norm += norm_curr
    end

    return schmidt_basis
end
#=}}}=#

max_stars_repo_name: nmayhall-vt/FermiCG

max_issues_repo_name: nmayhall-vt/FermiCG

max_forks_repo_name: nmayhall-vt/FermiCG


text:
#
# basics.jl -
#
# Implement interface to Tcl interpreter, evaluation of scripts, callbacks...
#

#------------------------------------------------------------------------------
# Make TclStatus behaves as an integer when used in comparisons with another
# status or an integer.

for op in (:(==), :(!=), :(<), :(<=), :(>), :(>=))
    @eval begin
        Base.$op(a::TclStatus, b::Integer  ) = $op(a.code, b     )
        Base.$op(a::Integer,   b::TclStatus) = $op(a     , b.code)
        Base.$op(a::TclStatus, b::TclStatus) = $op(a.code, b.code)
    end
end
TclStatus(status::TclStatus) = status
Base.convert(::Type{T}, status::TclStatus) where {T<:Integer} =
    convert(T, status.code)
Base.convert(::Type{Integer}, status::TclStatus) = status.code

#------------------------------------------------------------------------------
# Automatically name objects.

const __counter = Dict{String,Int}()

"""
    autoname(pfx = "jl_auto")

yields a unique name with given prefix.  The result is a string of the form
`pfx#` where `#` is a unique number for that prefix.
"""
function autoname(pfx::AbstractString = "jl_auto")
    global __counter
    n = get(__counter, pfx, 0) + 1
    __counter[pfx] = n
    return pfx*string(n)
end

#------------------------------------------------------------------------------
# Management of Tcl interpreters.

"""
When Tcl package is imported, an initial interpreter is created which can be
retrieved by:

```julia
interp = Tcl.getinterp()
```

A new Tcl interpreter can also be created by the command:

```julia
interp = TclInterp()
```

The resulting object can be used as a function to evaluate a Tcl script, for
instance:

```julia
interp("set x 42")
```

which yields the result of the script (here the string `"42"`).  An alternative
syntax is:

```
interp("set", "x", 42)
```

which yields the value `42`.  See method [`Tcl.eval`](@ref) for more details
about script evaluation.

The object can also be used as an array to access global Tcl variables (the
variable name can be specified as a string or as a symbol):

```julia
interp["x"]          # yields value of variable "x"
interp[:tcl_version] # yields version of Tcl
interp[:x] = 33      # set the value of "x" and yields its value
```

The Tcl interpreter is initialized and will be deleted when no longer in use.
If Tk has been properly installed, then:

```julia
interp("package require Tk")
```

should load Tk extension and create the "." toplevel Tk window.  But see
[`tkstart`](@ref) method to load Tk.

"""
function TclInterp(permanent::Bool=false)
    if (ptr = Tcl_CreateInterp()) == C_NULL
        Tcl.error("unable to create Tcl interpreter")
    end
    if Tcl_Init(ptr) != TCL_OK
        Tcl_DeleteInterp(ptr)
        Tcl.error("unable to initialize Tcl interpreter")
    end
    obj = TclInterp(ptr)
    permanent || finalizer(__finalize, obj)
    return obj
end

function __finalize(interp::TclInterp)
    Tcl_DeleteInterp(interp.ptr)
end

(interp::TclInterp)(::Type{T}, args...; kwds...) where {T} =
    Tcl.eval(T, interp, args...; kwds...)

(interp::TclInterp)(args...; kwds...) = Tcl.eval(interp, args...; kwds...)

isdeleted(interp::TclInterp) = Tcl_InterpDeleted(interp.ptr)
isactive(interp::TclInterp) = Tcl_InterpActive(interp.ptr)

"""
```julia
Tcl.getinterp()
```

yields the initial Tcl interpreter which is used by default by many methods.
An argument can be provided:

```julia
Tcl.getinterp(w)
```

yields the Tcl interpreter for widget (or image) `w`.

"""
getinterp() = __initial_interpreter[]

# This constant wraps a NULL pointer in a TclInterp structure.
const __NO_INTERP = TclInterp(TclInterpPtr(0))

# Many things do not work properly (e.g., freeing a Tcl object yield a
# segmentation fault) if no interpreter has been created, so we always create
# an initial Tcl interpreter (this is done by the __init__ method).
const __initial_interpreter = Ref{TclInterp}(__NO_INTERP)

# Interpreter pointer for callbacks and objects which need a Tcl interpreter.
const __context = [__NO_INTERP for i in 1:Threads.nthreads()]

"""
```julia
__interp()
```

yields the contextual Tcl interpreter of the calling thread.  Beware that this
may be a NULL interpreter.  The contextual Tcl interpreter shall only be used
for error messages and for creating callbacks.

See also: [`__intptr`](@ref), [`__set_context`](@ref),
          [`__reset_context`](@ref).

"""
@inline __interp() =
    # one-line definition does not work as expected...
    (@inbounds interp = __context[Threads.threadid()];
     return interp)

"""
```julia
__intptr()
```

yields a reference to the contextual Tcl interpreter of the calling thread.
Beware that this may be a NULL pointer.  The contextual Tcl interpreter shall
only be used for error messages and for creating callbacks.

See also: [`__interp`](@ref), [`__set_context`](@ref),
          [`__reset_context`](@ref).

"""
@inline __intptr() = __interp().ptr

"""
```julia
__reset_context()
```

resets the context for the calling thread.

See also: [`__interp`](@ref), [`__set_context`](@ref),
          [`__contextual_error`](@ref).

"""
@inline __reset_context() = __set_context(__NO_INTERP)

"""
```julia
__set_context(interp)
```

sets the contextual Tcl interpreter for the calling thread.

See also: [`__interp`](@ref), [`__reset_context`](@ref),
          [`__contextual_error`](@ref).

"""
@inline __set_context(interp::TclInterp) =
    @inbounds __context[Threads.threadid()] = interp

"""
```julia
__contextual_error(msg)
```

throws a Tcl error exception with message retrieved from the contextual
Tcl interpreter of the calling thread, if any, or given by `msg`, otherwise.

See also: [`__set_context`](@ref),  [`__reset_context`](@ref).

"""
__contextual_error(msg::String) =
    (intptr = __intptr();
     (intptr == C_NULL ? msg : __objptr_to(String, Tcl_GetObjResult(intptr))))

#------------------------------------------------------------------------------
# Evaluation of Tcl scripts.

"""
```julia
Tcl.setresult([interp,] args...) -> nothing
```

set result stored in Tcl interpreter `interp` or in the initial interpreter if
this argument is omitted.

"""
setresult() = setresult(getinterp())
setresult(arg) = setresult(getinterp(), arg)
setresult(args...) = setresult(getinterp(), args...)

# To set Tcl interpreter result, we can call `Tcl_SetObjResult` for any object,
# or `Tcl_SetResult` but only for string results with no embedded nulls.  There
# may be a slight advantage for calling `Tcl_SetResult` with non-volatile
# strings as copies are avoided.  Julia strings are immutable but I am not sure
# that they are non-volatile, so I prefer to not try using `Tcl_SetResult` and
# rather use `Tcl_SetObjResult` for any object.
#
# `Tcl_SetObjResult` does manage the reference count of its object argument so
# it is OK to directly pass a temporary object.
setresult(interp::TclInterp) =
    Tcl_SetObjResult(interp.ptr, __objptr())

setresult(interp::TclInterp, arg) =
    Tcl_SetObjResult(interp.ptr, __objptr(arg))

setresult(interp::TclInterp, args...) =
    Tcl_SetObjResult(interp.ptr, __newlistobj(args))

@static if false
    # The code for strings (taking care of embedded nulls and using
    # `Tcl_SetResult` if possible) is written below for reference but not
    # compiled.
    function setresult(interp::TclInterp, str::AbstractString;
                       volatile::Bool = true)
        ptr, siz = Base.unsafe_convert(Ptr{Cchar}, str), sizeof(str)
        if Base.containsnul(ptr, siz)
            # String has embedded NULLs, wrap it into a temporary object.  It
            # is not necessary to dela with its reference count as
            # Tcl_SetObjResult takes care of that.
            Tcl_SetObjResult(interp.ptr, __newstringobj(ptr, siz))
        else
            Tcl_SetResult(interp.ptr, ptr,
                          (volatile ? TCL_VOLATILE : TCL_STATIC))
        end
    end
end

"""
```julia
Tcl.getresult([T=Any,][interp])
```

yields the current result stored in Tcl interpreter `interp` or in the initial
interpreter if this argument is omitted.  If optional argument `T` is omitted,
the type of the returned value reflects that of the internal representation of
the result stored in Tcl interpreter; otherwise, `T` can be used to specify
how Tcl result should be converted (see [`Tcl.getvar`](@ref) for details).

See also: [`Tcl.getvar`](@ref).

"""
getresult() = getresult(getinterp())

getresult(interp::TclInterp) = getresult(Any, interp)

getresult(::Type{T}) where {T} = getresult(T, getinterp())

getresult(::Type{T}, interp::TclInterp) where {T} = __getresult(T, interp)

# Tcl_GetStringResult calls Tcl_GetObjResult, so we only interface to this
# latter function.  Incrementing the reference count of the result is only
# needed if we want to keep a long-term reference to it,
# `__objptr_to(TclObj,...)` takes care of that).
__getresult(::Type{T}, interp::TclInterp) where {T} =
    __objptr_to(T, Tcl_GetObjResult(interp.ptr))

"""
```julia
Tcl.exec([T,][interp,], args...; kwds...)
```

makes a list of the arguments `args...` and keywords `kwds...` and evaluates it
as a Tcl script.  See [`Tcl.eval`](@ref) for details about optional result type
`T` and Tcl interpreter `interp`.

Any specified keyword, say `key=val`, is automatically converted in the pair of
arguments `-key val` in this list (note the hyphen before the keyword name).
To allow for option names that are Julia keywords, a leading underscore is
stripped, if any, in `key`.  All keywords appear at the end of the list in
unspecific order.

Apart from the accounting of keywords, the main difference with
[`Tcl.eval`](@ref) is that each input argument is interpreted as a different
*word* of the Tcl script.  Using [`Tcl.eval`](@ref), [`Tcl.exec`](@ref) is
equivalent to:

```julia
Tcl.eval([T,][interp,], Tcl.list(args...; kwds...))
```

Specify `T` as [`TclStatus`](@ref), if you want to avoid throwing errors and
call [`Tcl.getresult`](@ref) to retrieve the result.

See also: [`Tcl.eval`](@ref), [`Tcl.list`](@ref), [`Tcl.getresult`](@ref).

"""
exec(args...; kwds...) =
    exec(getinterp(), args...; kwds...)

exec(::Type{T}, args...; kwds...) where {T} =
    exec(T, getinterp(), args...; kwds...)

exec(interp::TclInterp, args...; kwds...) =
    exec(Any, interp, args...; kwds...)

function exec(::Type{T}, interp::TclInterp, args...; kwds...) where {T}
    length(args) ≥ 1 || Tcl.error("expecting at least one argument")
    listptr = C_NULL
    __set_context(interp)
    try
        listptr = Tcl_IncrRefCount(__newlistobj())
        @__build_list listptr args kwds
        return __eval(T, interp, TclObj{List}, listptr)
    finally
        listptr == C_NULL || Tcl_DecrRefCount(listptr)
         __reset_context()
    end
end

"""
```julia
Tcl.eval([T,][interp,], args...)
```

concatenates arguments `args...` into a list and evaluates it as a Tcl script
with interpreter `interp` (or in the initial interpreter if this argument is
omitted).  See [`Tcl.concat`](@ref) for details about how arguments are
concatenated into a list.

If optional argument `T` is omitted, the type of the returned value reflects
that of the internal representation of the result of the script; otherwise, `T`
specifies the type of the result (see [`getvar`](@ref) for details).  As a
special case, when `T` is `TclStatus`, `Tcl.eval` behaves like the Tcl `catch`
command: the script is evaluated and no exception get thrown in case of error
but a status such as `TCL_OK` or `TCL_ERROR` is always returned and
[`Tcl.getresult`](@ref) can be used to retrieve the value of the result (which
is an error message if the returned status is `TCL_ERROR`).

Use `Tcl.exec` if you want to consider each argument in `args...` as a distinct
command argument.

Specify `T` as `TclStatus`, if you want to avoid throwing errors and
`Tcl.getresult` to retrieve the result.

See also: [`Tcl.concat`](@ref), [`Tcl.exec`](@ref), [`getvar`](@ref).

"""
Tcl.eval(args...) =
    Tcl.eval(getinterp(), args...)

Tcl.eval(::Type{T}, args...) where {T} =
    Tcl.eval(T, getinterp(), args...)

Tcl.eval(interp::TclInterp, args...) =
    Tcl.eval(Any, interp, args...)

function Tcl.eval(::Type{T}, interp::TclInterp,
                  script::Union{AbstractString,Symbol}) where {T}
    return Tcl.eval(T, interp, string(script))
end

function Tcl.eval(::Type{T}, interp::TclInterp, script::String) where {T}
    # A single string (or equivalent) is argument is to be evaluated.
    __set_context(interp)
    try
        return __eval(T, interp, String, script)
    finally
        __reset_context()
    end
end

function Tcl.eval(::Type{T}, interp::TclInterp, script::TclObj{S}) where{T,S}
    # A managed Tcl object is directly evaluated.  This saves the concatenation
    # which produces the same thing. FIXME: check this for a string.
    __set_context(interp)
    try
        return __eval(T, interp, TclObj{S}, __objptr(script))
    finally
        __reset_context()
    end
end

function Tcl.eval(::Type{T}, interp::TclInterp, args...) where {T}
    # Zero or more arguments, or a single non-list and non-string-like argument
    # are to be evaluated.  Concatenate all arguments and evaluate the
    # resulting list.
    length(args) ≥ 1 || Tcl.error("missing script to evaluate")
    listptr = C_NULL
    __set_context(interp)
    try
        listptr = Tcl_IncrRefCount(__newlistobj())
        @__concat_args listptr args
        return __eval(T, interp, TclObj{List}, listptr)
    finally
        listptr == C_NULL || Tcl_DecrRefCount(listptr)
        __reset_context()
    end
end

# Called in a properly set context, this method should produce consistent error
# messages and result of given type T.
function __eval(::Type{T}, interp::TclInterp, ::Type{S}, script) where {T,S}
    if __eval(TclStatus, interp, S, script) != TCL_OK
        Tcl.error(interp)
    end
    return __getresult(T, interp)
end

# For a not well defined Tcl object type, we peek the real object type to
# decide which of `Tcl_EvalObjEx` or `Tcl_EvalObjv` to call.
function __eval(::Type{TclStatus}, interp::TclInterp,
                ::Type{TclObj}, objptr::TclObjPtr)
    return __eval(TclStatus, interp, TclObj{__objtype(objptr)}, objptr)
end

# For a script given as a string, `Tcl_EvalEx` is slightly faster than
# `Tcl_Eval` (says the doc.) and, more importantly, the script may contain
# embedded nulls.
__eval(::Type{TclStatus}, interp::TclInterp, ::Type{String}, script::String) =
    Tcl_EvalEx(interp.ptr, Base.unsafe_convert(Ptr{Cchar}, script),
               sizeof(script), TCL_EVAL_GLOBAL|TCL_EVAL_DIRECT)

# For non-list objects, we use `Tcl_EvalObjEx` to evaluate a single argument
# script.  `Tcl_EvalObjEx` does manage the reference count of its object
# argument.
function __eval(::Type{TclStatus}, interp::TclInterp,
                ::Type{TclObj{S}}, objptr::TclObjPtr) where S
    flags = TCL_EVAL_GLOBAL
    if Tcl_GetRefCount(objptr) < 1
        # For a temporary object there is no needs to compile the script.
        flags |= TCL_EVAL_DIRECT
    end
    return Tcl_EvalObjEx(interp.ptr, objptr, flags)
end

function __eval(::Type{TclStatus}, interp::TclInterp,
                ::Type{TclObj{List}}, listptr::TclObjPtr)
    flags = TCL_EVAL_GLOBAL
    status, objc, objv = Tcl_ListObjGetElements(interp.ptr, listptr)
    if status == TCL_OK
        status = Tcl_EvalObjv(interp.ptr, objc, objv, flags)
    end
    return status
end

#------------------------------------------------------------------------------
# Exceptions

"""
```julia
Tcl.error(arg)
```

throws a [`TclError`](@ref) exception, argument `arg` can be the error message
as a string or a Tcl interpreter (in which case the error message is assumed to
be the current result of the Tcl interpreter).

"""
Tcl.error(msg::AbstractString) = throw(TclError(string(msg)))
Tcl.error(interp::TclInterp) = Tcl.error(getresult(String, interp))

"""
```julia
geterrmsg(ex)
```

yields the error message associated with exception `ex`.

"""
geterrmsg(ex::Exception) = sprint(io -> showerror(io, ex))

#------------------------------------------------------------------------------
# Processing Tcl/Tk events.  The function `doevents` must be repeatedly
# called too process events when Tk is loaded.

const __timer = Ref{Timer}()

"""
```julia
Tcl.isrunning()
```

yields whether the processing of Tcl/Tk events is running.

"""
isrunning() = (isdefined(__timer, 1) && isopen(__timer[]))

"""
```julia
Tcl.resume(sec=0.05)
```

resumes or starts the processing of Tcl/Tk events with an interval of `sec`
seconds.  This manages to repeatedly call function `Tcl.doevents`.  The method
`Tcl.suspend` can be called to suspend the processing of events.

Calling `Tcl.resume` is mandatory when Tk extension is loaded.  Thus:

```julia
Tcl.eval(interp, "package require Tk")
Tcl.resume()
```

is the recommended way to load Tk package.  Alternatively:

```julia
Tcl.tkstart()
```

can be called to do that.

"""
resume(sec::Real=0.05) =
    (isrunning() || (__timer[] = Timer(doevents, 0.1; interval=sec)); nothing)

"""
```julia
Tcl.suspend()
```

suspends the processing of Tcl/Tk events for all interpreters.  The method
[`Tcl.resume`](@ref) can be called to resume the processing of events.

"""
suspend() =
    (isrunning() && close(__timer[]); nothing)

"""
```julia
Tcl.doevents(flags = TCL_DONT_WAIT|TCL_ALL_EVENTS)
```

processes Tcl/Tk events for all interpreters.  Normally this is automatically
called by the timer set by [`Tcl.resume`](@ref).

"""
doevents(::Timer) = doevents()

function doevents(flags::Integer = TCL_DONT_WAIT|TCL_ALL_EVENTS)
    while Tcl_DoOneEvent(flags)
    end
end

#------------------------------------------------------------------------------
# Implement callbacks.

# Dictionary of objects shared with Tcl to make sure they are not garbage
# collected until Tcl deletes their reference.
const __references = Dict{Any,Int}()

function preserve(obj)
    __references[obj] = get(__references, obj, 0) + 1
end

function release(obj)
    if haskey(__references, obj)
        if __references[obj] > 1
            __references[obj] -= 1
        else
            pop!(__references, obj)
        end
    end
    nothing
end

const __releaseobject_proc = Ref{Ptr{Cvoid}}() # will be set by __init__

__releaseobject(ptr::Ptr{Cvoid}) = release(unsafe_pointer_to_objref(ptr))

const __evalcommand_proc = Ref{Ptr{Cvoid}}() # will be set by __init__

function __evalcommand(fptr::ClientData, iptr::TclInterpPtr,
                       objc::Cint, objv::Ptr{TclObjPtr}) :: Cint
    func = unsafe_pointer_to_objref(fptr)
    interp = TclInterp(iptr) # weak reference
    __set_context(interp) # FIXME: not needed?
    try
        args = __buildvector(objc, objv)
        return __setcommandresult(interp, func(interp, args...))
    catch ex
        setresult(interp, "(callback error) " * geterrmsg(ex))
        return TCL_ERROR
    finally
        __reset_context() # FIXME: not needed?
    end
end

# If the function provides a return code, we do want to return it to the
# interpreter, otherwise TCL_OK is assumed.
__setcommandresult(interp::TclInterp, args...) =
    __setcommandresult(interp, TCL_OK, args...)

__setcommandresult(interp::TclInterp, args::Tuple) =
    __setcommandresult(interp, TCL_OK, args...)

__setcommandresult(interp::TclInterp, args::Tuple{TclStatus,Vararg}) =
    __setcommandresult(interp, args[1], args[2:end]...)

__setcommandresult(interp::TclInterp, status::TclStatus) =
    (Tcl_SetObjResult(interp.ptr, __newobj()); return status)

__setcommandresult(interp::TclInterp, status::TclStatus, arg) =
    (Tcl_SetObjResult(interp.ptr, __newobj(arg)); return status)

__setcommandresult(interp::TclInterp, status::TclStatus, args...) =
    (Tcl_SetObjResult(interp.ptr, __newobj(args)); return status)

# With precompilation, `__init__()` carries on initializations that must occur
# at runtime like `@cfunction` which returns a raw pointer.
function __init__()
    __initial_interpreter[] = TclInterp(true)
    __releaseobject_proc[] = @cfunction(__releaseobject, Cvoid, (Ptr{Cvoid},))
    __evalcommand_proc[] = @cfunction(__evalcommand, Cint,
                                      (ClientData, TclInterpPtr,
                                       Cint, Ptr{TclObjPtr}))
    __init_types()
end

"""
```julia
Tcl.CallBack([interp,] [name,] f) -> name
```

creates a callback command named `name` in Tcl interpreter `interp` (or in the
initial Tcl interpreter if this argument is omitted).  If `name` is missing
`Tcl.autoname("jl_callback")` is used to automatically define a name.  The
command name is returned as a string.  The Tcl command will call the Julia
function `f` as follows:

```julia
f(interp, name, args...) -> [status::TclStatus], vals...
```

where `interp` is the Tcl interpreter which calls the command, `name` is the
command name and `args...` are the arguments of the command.

The function can return any number of values.  The status retuned to the
interpreter by the command is assumed to be the first of these values if its
type is `TclStatus` (one of `TCL_OK`, `TCL_ERROR`, `TCL_RETURN`, `TCL_BREAK` or
`TCL_CONTINUE`); otherwise, `TCL_OK` is assumed.  The other values are stored
in the interpreter's result.

If the function throws any exception, the error message associated with the
exception is stored in the interpreter's result and `TCL_ERROR` is retuned to
the interpreter.

See also: [`Tcl.deletecommand`](@ref), [`Tcl.autoname`](@ref).

"""
Callback(func::Function) =
    Callback(getinterp(), func)

Callback(name::StringOrSymbol, func::Function) =
    Callback(getinterp(), name, func)

Callback(interp::TclInterp, func::Function) =
    __newcallback(interp.ptr,func)

Callback(interp::TclInterp, name::StringOrSymbol, func::Function) =
    __newcallback(interp.ptr, name, func)

__newcallback(func::Function) =
    __newcallback(__intptr(), func)

__newcallback(intptr::TclInterpPtr, func::Function) =
    __newcallback(intptr, autoname("jl_callback"), func)

__newcallback(intptr::TclInterpPtr, name::Symbol, func::Function) =
    __newcallback(intptr, string(name), func)

function __newcallback(intptr::TclInterpPtr, name::String, func::Function)
    local cmd :: TclObj{Function}
    if intptr != C_NULL
        preserve(func)
        token = Tcl_CreateObjCommand(intptr, name, __evalcommand_proc[],
                                     pointer_from_objref(func),
                                     __releaseobject_proc[])
        if token == C_NULL
            release(func)
            Tcl.error(__objptr_to(String, Tcl_GetObjResult(intptr)))
        end
        cmd = TclObj{Function}(__newobj(""))
        Tcl_GetCommandFullName(intptr, token, __objptr(cmd))
    else
        cmd = TclObj{Function}(__newobj(name))
    end
    return Callback(intptr, cmd, func)
end

TclObj(func::Function) = Callback(func)

AtomicType(::Type{<:Union{Function,Callback}}) = Atomic()

__objptr(func::Function) = __objptr(__newobj(func))

__newobj(func::Function) = __objptr(__newcallback(func))

__objptr(func::Callback) = __objptr(func.obj)

"""
```julia
Tcl.deletecommand([interp,] cmd)
```

deletes a command named `cmd` in Tcl interpreter `interp` (or in the initial
Tcl interpreter if this argument is omitted). `cmd` can also be a callback.

See also: [`Tcl.Callback`](@ref).

"""
deletecommand(cmd::Union{StringOrSymbol,Callback}) =
    # Since a callback just has a weak reference to the interpreter where it
    # was created, we must not use it as a default.
    deletecommand(getinterp(), cmd)

deletecommand(interp::TclInterp, cmd::Callback) =
    deletecommand(interp, cmd.name)

function deletecommand(interp::TclInterp, name::StringOrSymbol)
    if Tcl_DeleteCommand(interp.ptr, name) != TCL_OK
        Tcl.error(interp)
    end
    return nothing
end

ext: jl


text:
export
        FloatLeft,   FloatRight,
        Relative, Fixed, Absolute, HasAbsolute,
        Bottom, Right,
        HasOpacity, HasImage, HasPath,
        TextPath, LinearGrad, RadialGrad,

        RowFinalized,
        HasBorder,   BordersSame,  Clip, IsRoundBox,
        FixedHeight,
        DisplayBlock, DisplayInlineBlock, DisplayInline, DisplayNone, DisplayTable, DisplayFlex,
        LineBreakBefore, LineBreakAfter,
        IsHidden,
        AlignBase,   AlignMiddle,
        TextCenter,  TextRight,    TextJustify,
        TextItalic,  TextOblique,  TextBold,
        IsHScroll, IsVScroll,
        Percent,
        Marked

const FloatLeft, FloatRight, RowFinalized = 1, 2, 3
const Absolute, HasAbsolute, Relative, Fixed = 4, 5, 6, 7
const Bottom, Right = 8, 9
const HasOpacity    = 10
# position: static, relative, fixed, absolute
const HasBorder, BordersSame, Clip, IsRoundBox, FixedHeight    = 11, 12, 13, 14, 15
# https://developer.mozilla.org/en-US/docs/Web/CSS/display
const DisplayBlock, DisplayInlineBlock, DisplayInline,
      DisplayNone, DisplayTable, DisplayFlex = 16, 17, 18, 19, 20, 21
# visibility: visible|hidden|collapse|initial|inherit;
const IsHidden     = 22
const AlignBase, AlignMiddle = 23, 24
const TextCenter, TextRight, TextJustify = 25, 26, 27
const TextItalic, TextOblique = 28, 29
const TextBold       = 30
const IsHScroll, IsVScroll = 31, 32

const HasImage      = 33
const TextPath      = 34
const LinearGrad, RadialGrad = 35, 36

const LineBreakBefore, LineBreakAfter = 37, 38

const Percent      = 39
const Marked       = 39

const HasPath      = 40

#.==============================================================================


text:
"""
```
filterProperNodes( distNodes::Dict{T, Float64} ) where T <: Union{Int, Float64}
```
This function filters the bad nodes from the list of distribution nodes in `distNodes`. A bad node is a node with a negative node value or a negative weight.

The function returns a `Dict{T, Float64}`, the filtered list of distribution nodes.
"""
function filterProperNodes( distNodes::Dict{T,Float64} ) where T <: Union{Int, Float64}

    tmpNodes = deepcopy( distNodes )

    for node in keys( distNodes )
        if ( node < 0 ) || ( distNodes[node] < 0 )
            delete!( tmpNodes, node )
        end  # if ( node < 0 ) || ...
    end  # for node in keys( distNodes )

    return tmpNodes

end  # filterProperNodes( distNodes )

function setRecDist!( recruitment::Recruitment, distType::Symbol,
    distNodes::Dict{T,Float64} ) where T <: Union{Int, Float64}

    # Check if the distribution type is known.
    if !haskey( recruitmentDists, distType )
        @warn "Unknown distribution type, not making any changes."
        return false
    end  # if !haskey( recruitmentDists, distType )

    tmpDistNodes = filterProperNodes( distNodes )
    tmpNodes = sort( collect( keys( tmpDistNodes ) ) )

    # Check if there are sufficient nodes.
    if length( tmpNodes ) < recruitmentDists[distType][2]
        @warn "Not enough proper distribution nodes, not making any changes."
        return false
    end  # if length( tmpNodes ) < recruitmentDists[distType][2]

    # Check if the total probability mass is non-zero.
    pMass = sum( map( node -> tmpDistNodes[node], tmpNodes ) )

    if ( pMass == 0.0 ) || ( ( distType === :pUnif ) &&
        ( distNodes[tmpNodes[end]] == pMass ) )
        @warn "Proposed distribution has 0 probability mass, not making any changes."
        return false
    end  # if ( pMass == 0.0 ) || ...

    if T === Int
        recDist = recruitmentDists[distType][3]
        recruitment.isAdaptive = false
        recruitment.recruitmentDistType = distType
        recruitment.recruitmentDistNodes = tmpDistNodes
        recruitment.recruitmentDist = recDist( recruitment, tmpDistNodes,
            tmpNodes )
    else
        ageDist = recruitmentDists[distType][4]
        recruitment.ageDistType = distType
        recruitment.ageDistNodes = tmpDistNodes
        recruitment.ageDist = ageDist( recruitment, tmpDistNodes, tmpNodes )
    end  # if T === Int

    return true

end  # setRecDist!( recruitment, distType, distNodes )

function recDiscDist( recruitment::Recruitment, distNodes::Dict{Int,Float64},
    nodes::Vector{Int} )

    # No need to involve a distribution if there's only one node.
    if length( nodes ) == 1
        return function() return nodes[1] end
    end  # if length( nodes ) == 1

    # Get the point probabilities of the nodes.
    pNodes = map( node -> distNodes[node], nodes )
    pNodes /= sum( pNodes )

    return function()
        return nodes[rand( recruitment.recRNG, Categorical( pNodes ) )]
    end  # anonymous function()

end  # recDiscDist( recruitment, distNodes, nodes )

function recPUnifDist( recruitment::Recruitment, distNodes::Dict{Int,Float64},
    nodes::Vector{Int} )

    if length( nodes ) == 2
        intLength = nodes[2] - nodes[1]

        return function()
            return intLength == 1 ? nodes[1] :
                rand( recruitment.recRNG, nodes[1]:( nodes[2] - 1 ) )
        end  # anonymous function()
    end  # if length( nodes ) == 2

    # Get the point probabilities of the intervals.
    pInts = map( node -> distNodes[node], nodes[1:(end-1)] )
    pInts /= sum( pInts )

    return function()
        intInd = rand( recruitment.recRNG, Categorical( pInts ) )
        intLength = nodes[intInd + 1] - nodes[intInd]
        return intLength == 1 ? nodes[intInd] :
            rand( recruitment.recRNG, nodes[intInd]:( nodes[intInd + 1] - 1 ) )
    end  # anonymous function()

end  # recPUnifDist( recruitment, distNodes, nodes )

function recPLinDist( recruitment::Recruitment, distNodes::Dict{Int,Float64},
    nodes::Vector{Int} )

    # Get the point probabilities of the intervals.
    pointWeights = map( node -> distNodes[node], nodes )
    nodeDiffs = nodes[2:end] .- nodes[1:(end-1)] .- 1
    bracketWeights = map( eachindex( nodeDiffs ) ) do ii
        return pointWeights[ii] + nodeDiffs[ii] * ( pointWeights[ii] + pointWeights[ii + 1] ) / 2
    end  # map( eachindex( nodeDiffs ) ) do ii
    push!( bracketWeights, pointWeights[end] )
    bracketProbs = bracketWeights / sum( bracketWeights )

    return function()
        intInd = rand( recruitment.recRNG, Categorical( bracketProbs ) )

        if intInd == length( bracketProbs )
            return nodes[end]
        end  # if intInd == length( bracketProbs )

        posProbs = ( pointWeights[intInd + 1] - pointWeights[intInd] ) /
            ( nodeDiffs[intInd] + 1 )
        posProbs *= collect( 0:nodeDiffs[intInd] )
        posProbs = posProbs .+ pointWeights[intInd]
        posProbs /= bracketWeights[intInd]
        return nodes[intInd] + rand( recruitment.recRNG,
            Categorical( posProbs ) ) - 1
    end  # anonymous function()

end  # recPLinDist( recruitment, distNodes, nodes )

function ageDiscDist( recruitment::Recruitment,
    distNodes::Dict{Float64,Float64}, nodes::Vector{Float64} )

    # No need to involve a distribution if there's only one node.
    if length( nodes ) == 1
        return function( n::Integer ) return fill( nodes[1], n ) end
    end  # if length( nodes ) == 1

    # Get the point probabilities of the nodes.
    pNodes = map( node -> distNodes[node], nodes )
    pNodes /= sum( pNodes )

    return function( n::Integer )
        return nodes[rand( recruitment.ageRNG, Categorical( pNodes ), n )]
    end  # anonymous function( n )

end  # ageDiscDist( recruitment, distNodes, nodes )

function agePUnifDist( recruitment::Recruitment,
    distNodes::Dict{Float64,Float64}, nodes::Vector{Float64} )

    if length( nodes ) == 2
        intLength = nodes[2] - nodes[1]

        return function( n::Integer )
            return rand( recruitment.ageRNG, n ) * intLength .+ nodes[1]
        end  # anonymous function( n )
    end  # if length( nodes ) == 2

    # Get the point probabilities of the intervals.
    pInts = map( node -> distNodes[node], nodes[1:(end-1)] )
    pInts /= sum( pInts )

    return function( n::Integer )
        intInds = rand( recruitment.ageRNG, Categorical( pInts ), n )
        intLengths = nodes[intInds .+ 1] - nodes[intInds]
        return rand( recruitment.ageRNG, n ) .* intLengths + nodes[intInds]
    end  # anonymous function( n )

end  # agePUnifDist( recruitment, distNodes, nodes )

function agePLinDist( recruitment::Recruitment,
    distNodes::Dict{Float64,Float64}, nodes::Vector{Float64} )

    # Get the probabilities of each interval.
    pointWeights = map( node -> distNodes[node], nodes )
    bracketWeights = ( nodes[2:end] - nodes[1:(end - 1)] ) .*
        ( pointWeights[1:(end - 1)] + pointWeights[2:end] ) / 2
    bracketProbs = bracketWeights / sum( bracketWeights )

    return function( n::Integer )
        intInds = rand( recruitment.ageRNG, Categorical( bracketProbs ), n )
        isDiffWeights = pointWeights[intInds .+ 1] .!= pointWeights[intInds]
        diffInds = intInds[isDiffWeights]
        diffIndsP = diffInds .+ 1
        sameInds = intInds[.!isDiffWeights]
        result = rand( recruitment.ageRNG, n )

        # First, add the entries with equal weights in both endpoints of the
        #   containing interval.
        result[.!isDiffWeights] = nodes[sameInds] +
            result[.!isDiffWeights] .*
            ( nodes[sameInds .+ 1] - nodes[sameInds] )

        # Then, do the same for the other entries.
        a = ( pointWeights[diffIndsP] - pointWeights[diffInds] ) / 2
        b = pointWeights[diffInds] .* nodes[diffIndsP] -
            pointWeights[diffIndsP] .* nodes[diffInds]
        c = ( pointWeights[diffInds] + pointWeights[diffIndsP] ) .*
            nodes[diffInds].^2 / 2 - pointWeights[diffInds] .*
            nodes[diffInds] .* nodes[diffIndsP] -
            bracketWeights[diffInds] .* result[isDiffWeights] .*
            ( nodes[diffIndsP] - nodes[diffInds] )
        d = b.^2 - 4 * a .* c
        result[isDiffWeights] = ( sqrt.( d ) - b ) ./ ( 2 * a )

        return result
    end  # anonymous function( n )

end  # agePLinDist( recruitment, distNodes, nodes )

const recruitmentDists = Dict(
    :disc => ("discrete", 1, recDiscDist, ageDiscDist),
    :pUnif => ("piecewise uniform", 2, recPUnifDist, agePUnifDist),
    :pLin => ("piecewise linear", 2, recPLinDist, agePLinDist) )


text:
using OutlierDetectionPython
using OutlierDetectionTest

# Test the metadata of all exported detectors
test_meta.(eval.(OutlierDetectionPython.MODELS))

data = TestData()
run_test(detector) = test_detector(detector, data)

# ABOD
run_test(ABODDetector())

# CBLOF
run_test(CBLOFDetector(random_state=0))

# COF
run_test(COFDetector())

# COPOD
run_test(COPODDetector())

# HBOS
run_test(HBOSDetector())

# IForest
run_test(IForestDetector(random_state=0))

# KNN
run_test(KNNDetector())

# LMDD
run_test(LMDDDetector(random_state=0))

# LOD
run_test(LODADetector())

# LOF
run_test(LOFDetector())

# LOCI
run_test(LOCIDetector())

# MCD
run_test(MCDDetector(random_state=0))

# SVM
run_test(OCSVMDetector())

# PCA
run_test(PCADetector(random_state=0))

# ROD
run_test(RODDetector())

# SOD
run_test(SODDetector())

# SOS
run_test(SOSDetector())


text:
# ChainRulesTestUtils requires some method override to support sparse array
# BEGIN ChainRulesTestUtils

# we use this to generate Tangent and pass to test_rrule
# for example: A ⊢ sprand_tangent(A)
using SparseADRules:Normal_QR

function sprand_tangent(A::AT) where AT<:SparseMatrixCSC
    return ChainRulesTestUtils.rand_tangent(A)
end

function sprand_tangent(A::Adjoint{T, <:AbstractSparseMatrix}) where T
    # when ProjectTo is supported for Adjoint, copy is not necessary
    # https://github.com/JuliaDiff/ChainRulesCore.jl/issues/447
    return ChainRulesTestUtils.rand_tangent(copy(A))
end

function sprand_tangent(x::SparseVector)
    return ChainRulesTestUtils.rand_tangent(x)
end

function Base.:(+)(a::P, b::Tangent{P}) where P<:SparseMatrixCSC
    b[1] .+ a
    return b
end

function FiniteDifferences.to_vec(A::SparseMatrixCSC{Tv, Ti}) where {Tv, Ti}
    v = reinterpret(real(Tv), A.nzval)
    function SparseMatrixCSC_from_vec(v)
        m = A.m
        n = A.n
        colptr = A.colptr
        rowval = A.rowval
        cv = Vector(reinterpret(Tv, v))
        SparseMatrixCSC(m, n, colptr, rowval, cv)
    end
    return v, SparseMatrixCSC_from_vec
end

function FiniteDifferences.to_vec(x::SparseVector{Tv, Ti}) where {Tv, Ti}
    v= reinterpret(real(Tv), x.nzval)
    function SparseVector_from_vec(v)
        n = x.n
        nzind = x.nzind
        cv = Vector(reinterpret(Tv, v))
        SparseVector(n, nzind, cv)
    end
    return v, SparseVector_from_vec
end

function FiniteDifferences.to_vec(A::Adjoint{T, <:AbstractSparseMatrix}) where T
    Ā = copy(A)
    return FiniteDifferences.to_vec(Ā)
end

function FiniteDifferences.to_vec(x::S) where {S <: Normal_QR }
    q_vec, q_back = FiniteDifferences.to_vec(x.Q)
    r_vec, r_back = FiniteDifferences.to_vec(x.R)
    function Normal_QR_from_vec(v)
        Q_new = q_back(v.Q)
        R_new = r_back(v.R)
        return S(Q_new, R_new)
    end
    return [q_vec; r_vec], Normal_QR_from_vec
end

function FiniteDifferences._j′vp(fdm, f, ȳ::Vector{<:Number}, x::Vector{<:Complex})
    isempty(x) && return eltype(ȳ)[] # if x is empty, then so is the jacobian and x̄
    r_jcb = transpose(first(FiniteDifferences.jacobian(fdm, f, real(x)))) * real(ȳ)
    i_jcb = transpose(first(FiniteDifferences.jacobian(fdm, f, imag(x)))) * imag(ȳ)
    return Complex(r_jcb, i_jcb)
end

# END ChainRulesTestUtils


text:
using GreatCircle
using Base.Test

# One test decimal degree is 111000m
latitude = 40.0
longitude = -76.0

"""
    Great Circle Tests
"""
azimuth = 90.
new_gc = great_circle(111000., azimuth, latitude, longitude)
# We should have gone to the right
@test new_gc["longitude"] > longitude + 0.9

azimuth = 270.
new_gc = great_circle(111000., azimuth, latitude, longitude)
# We should have gone to the left
@test new_gc["longitude"] < longitude - 0.9

azimuth = 180.
new_gc = great_circle(111000., azimuth, latitude, longitude)
# We should have gone down
@test new_gc["latitude"] < latitude - 0.9

azimuth = 0.
new_gc = great_circle(111000., azimuth, latitude, longitude)
# We should have gone up
@test new_gc["latitude"] > latitude + 0.9

azimuth = 315.
new_gc = great_circle(111000., azimuth, latitude, longitude)
# We should have gone up and to the left
@test new_gc["latitude"] > latitude + 0.45
@test new_gc["longitude"] < longitude - 0.45

# One decimal degree is 111000m
latitude = [40.0:5.:60.0]
longitude = [-40.0:-5.:-60.0]

azimuth = ones(longitude) * 90.
new_gc = great_circle(111000., azimuth, latitude, longitude)
# We should have gone to the right
@test all(new_gc["longitude"] .> longitude + 0.9)

azimuth = 270.
distance = 111000. * ones(longitude)
new_gc = great_circle(distance, azimuth, latitude, longitude)
# We should have gone to the left
@test all(new_gc["longitude"] .< longitude - 0.9)

azimuth = 180.
new_gc = great_circle(111000., azimuth, latitude, longitude)
# We should have gone down
@test all(new_gc["latitude"] .< latitude - 0.9)

azimuth = 0.
new_gc = great_circle(111000., azimuth, latitude, longitude)
# We should have gone up
@test all(new_gc["latitude"] .> latitude + 0.9)

azimuth = 315.
new_gc = great_circle(111000., azimuth, latitude, longitude)
# We should have gone up and to the left
@test all(new_gc["latitude"] .> latitude + 0.45)
@test all(new_gc["longitude"] .< longitude - 0.45)

"""
    Great Distance Tests
"""
# One decimal degree at the equator is about 111.32km
latitude_start  = 0.
latitude_end    = 0.
longitude_start = 50.
longitude_end   = 52.
gd = great_distance(latitude_start, longitude_start, latitude_end, longitude_end)
@test round(gd["distance"] ./ 1000, 2) .== 111.32 .* 2

# One decimal degree is 111000m
latitude_start  = 0.
latitude_end    = 0.
longitude_start = [49., 75.]
longitude_end   = [50., 76.]
gd = great_distance(latitude_start, longitude_start, latitude_end, longitude_end)
@test all( round(gd["distance"] ./ 1000, 2) .== 111.32 )


text:
module MyExample

greet()=print("Hello World")
include("extra_file.jl")
end


text:
using NDTools # this is needed for the select_region! function below.

@testset "chirp z-transformation" begin
    @testset "czt" begin
        x = randn(ComplexF32, (5,6,7))
        y = randn(ComplexF32, (5,6))
        zoom = (1.0,1.0,1.0)
        @test ≈(czt(x, zoom), ft(x),rtol=1e-4)
        @test ≈(czt(y, zoom), ft(y),rtol=1e-5)
        @test ≈(iczt(czt(y,zoom),zoom),  y, rtol=1e-5)
        zoom = (2.0,2.0)
        @test ≈(czt(y,zoom),  select_region(upsample2(ft(y), fix_center=true),new_size=size(y)), rtol=1e-5)
        # @vt czt(y,zoom)  select_region(upsample2(ft(y), fix_center=true), new_size=size(y))
    end
end


text:
const shift = 1
const reducel = 2
const reducer = 3

type State
    tokens::Vector
    step::Int
    top::Int
    left
    right::Int
    lcs::Vector{Int}
    rcs::Vector{Int}
    prev
    act::Int
end
State(tokens) = State(tokens, 1, 1, nothing, 2, Int[], Int[], nothing, 0)

isfinal(s::State) = s.step == ((length(s.tokens)-1) * 2 + 1)

function expand(s::State)
    isfinal(s) && return State[]
    acts = Int[]
    s.left == nothing || push!(acts, reducel, reducer)
    s.right <= length(s.tokens) && push!(acts, shift)
    a = acts[rand(1:length(acts))]
    State(s, a)
end

function State(s::State, act::Int)
    top, lcs, rcs, left, right = begin
        if act == shift
            s.right, Int[], Int[], s, s.right+1
        elseif act == reducel
            _lcs = length(s.lcs) == 0 ? [s.left.top] : [s.left.top,s.lcs[1]]
            s.top, _lcs, s.rcs, s.left.left, s.right
        elseif act == reducer
            _rcs = length(s.left.rcs) == 0 ? [s.top] : [s.top,s.left.rcs[1]]
            s.left.top, s.left.lcs, _rcs, s.left.left, s.right
        else
            throw("Invalid action: $(act)")
        end
    end
    State(s.tokens, s.step+1, top, left, right, lcs, rcs, s, act)
end

function depparse(tokens::Vector)
    state = State(tokens)
    while !isfinal(state)
        state = expand(state)
    end
    state
end

function toheads(s::State)
    heads = fill(0, length(s.tokens))
    while s != nothing
        length(s.lcs) > 0 && (heads[s.lcs[1]] = s.top)
        length(s.rcs) > 0 && (heads[s.rcs[1]] = s.top)
        s = s.prev
    end
    heads
end

max_stars_repo_name: hshindo/react-nlp

max_issues_repo_name: hshindo/react-nlp

max_forks_repo_name: hshindo/react-nlp


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""DeploymentStrategy describes how to replace existing pods with new ones.

    IoK8sApiExtensionsV1beta1DeploymentStrategy(;
        rollingUpdate=nothing,
        type=nothing,
    )

    - rollingUpdate::IoK8sApiExtensionsV1beta1RollingUpdateDeployment : Rolling update config params. Present only if DeploymentStrategyType = RollingUpdate.
    - type::String : Type of deployment. Can be \"Recreate\" or \"RollingUpdate\". Default is RollingUpdate.
"""
mutable struct IoK8sApiExtensionsV1beta1DeploymentStrategy <: SwaggerModel
    rollingUpdate::Any # spec type: Union{ Nothing, IoK8sApiExtensionsV1beta1RollingUpdateDeployment } # spec name: rollingUpdate
    type::Any # spec type: Union{ Nothing, String } # spec name: type

    function IoK8sApiExtensionsV1beta1DeploymentStrategy(;rollingUpdate=nothing, type=nothing)
        o = new()
        validate_property(IoK8sApiExtensionsV1beta1DeploymentStrategy, Symbol("rollingUpdate"), rollingUpdate)
        setfield!(o, Symbol("rollingUpdate"), rollingUpdate)
        validate_property(IoK8sApiExtensionsV1beta1DeploymentStrategy, Symbol("type"), type)
        setfield!(o, Symbol("type"), type)
        o
    end
end # type IoK8sApiExtensionsV1beta1DeploymentStrategy

const _property_map_IoK8sApiExtensionsV1beta1DeploymentStrategy = Dict{Symbol,Symbol}(Symbol("rollingUpdate")=>Symbol("rollingUpdate"), Symbol("type")=>Symbol("type"))
const _property_types_IoK8sApiExtensionsV1beta1DeploymentStrategy = Dict{Symbol,String}(Symbol("rollingUpdate")=>"IoK8sApiExtensionsV1beta1RollingUpdateDeployment", Symbol("type")=>"String")
Base.propertynames(::Type{ IoK8sApiExtensionsV1beta1DeploymentStrategy }) = collect(keys(_property_map_IoK8sApiExtensionsV1beta1DeploymentStrategy))
Swagger.property_type(::Type{ IoK8sApiExtensionsV1beta1DeploymentStrategy }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiExtensionsV1beta1DeploymentStrategy[name]))}
Swagger.field_name(::Type{ IoK8sApiExtensionsV1beta1DeploymentStrategy }, property_name::Symbol) =  _property_map_IoK8sApiExtensionsV1beta1DeploymentStrategy[property_name]

function check_required(o::IoK8sApiExtensionsV1beta1DeploymentStrategy)
    true
end

function validate_property(::Type{ IoK8sApiExtensionsV1beta1DeploymentStrategy }, name::Symbol, val)
end


text:
mutable struct Parallel <: Geometry
    label::String
    shape::Vector{Line}
    nsample::Int
end

function Parallel( ; label, shape, nsample)
    return Parallel(label, shape, nsample)
end


text:
function betaBinomial(prior::P, bernoulliTrials::Vector{Int}) where {P <: Beta}
  n = length(bernoulliTrials)
  successes = sum(bernoulliTrials)
  failures = n - successes
  a = prior.α + successes
  b = prior.β + failures
  return Beta(a, b)
end
#betaBinomial(Beta(2,3), [1,0,1,1])


text:
# Autogenerated wrapper script for Minuit2_jll for x86_64-w64-mingw32-cxx11
export libMinuit2, libMinuit2Math

JLLWrappers.@generate_wrapper_header("Minuit2")
JLLWrappers.@declare_library_product(libMinuit2, "libMinuit2.dll")
JLLWrappers.@declare_library_product(libMinuit2Math, "libMinuit2Math.dll")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libMinuit2,
        "bin\\libMinuit2.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libMinuit2Math,
        "bin\\libMinuit2Math.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using Random # Random stdlib
import ForwardDiff: derivative
import StatsBase: sample, weights # To include the basic sample from array function
# export Line, Envelop, RejectionSampler # Structures/classes
# export run_sampler!, sample_envelop, eval_envelop, add_segment! # Methods

"""
    Line(slope::Float64, intercept::Float64)
Basic tangent ensamble-unit for an envelop.
"""
struct Line
    slope::Float64
    intercept::Float64
end
eval_line(l::Line, x::Float64) = l.slope * x + l.intercept

"""
    intersection(l1::Line, l2::Line)
Finds the horizontal coordinate of the intersection between lines
"""
function intersection(l1::Line, l2::Line)
    @assert l1.slope ≥ l2.slope "slopes should be weakly ordered (l1: $(l1.slope), l2: $(l2.slope))"
    @assert l1.slope != l2.slope
    - (l2.intercept - l1.intercept) / (l2.slope - l1.slope)
end

"""
    exp_integral(l::Line, x1::Float64, x2::Float64)
Computes the integral
    ``LaTeX \\int_{x_1} ^ {x_2} \\exp\\{ax + b\\} dx. ``
The resulting value is the weight assigned to the segment [x1, x2] in the envelop
"""
function exp_integral(l::Line, x1::Float64, x2::Float64)::Tuple{Float64, Float64}
    (x1 == x2) && (return (0.0, 0.0))
    @assert x1 ≤ x2 "invalid x1=$x1, x2=$x2, $l"
    a, b = l.slope, l.intercept
    # for numerical stability: multiply by 0.5 * a(x1 - x2)
    if a != 0.0
        z1, z2 = eval_line(l, x1), eval_line(l, x2)
        scale = max(z1, z2)
        out = (exp(z2 - scale) - exp(z1 - scale)) / a
    else
        scale = 0.0
        out = (x2 - x1) * exp(b)
    end
    @assert out ≥ 0 "negative weight (out: $out, x1: $x1, x2: $x2, l: $l, scale: $scale"
    # debug info
    #  println("  @exp_integral:")
    #  println("  |-- line: $l")
    #  println("  |-- x1: $x1")
    #  println("  |-- x2: $x2")
    #  println("  |-- exp(b + a * x2): $(exp(b + a * x2))")
    #  println("  |-- exp(b + a * x1): $(exp(b + a * x1))")
    #  println("  |-- exp(a * x2): $(exp(a * x2))")
    #  println("  |-- exp(a * x1): $(exp(a * x1))")
    #  println("  |-- out: $out")
    return out, scale
end

"""
    Envelop(lines::Vector{Line}, support::Tuple{Float64, Float64})
A piecewise linear function with k segments defined by the lines `L_1, ..., L_k` and cutpoints
`c_1, ..., c_k+1` with `c1 = support[1]` and `c2 = support[2]`. A line L_k is active in the segment
[c_k, c_k+1], and it's assigned a weight w_k based on [exp_integral](@exp_integral). The weighted integral
over c_1 to c_k+1 is one, so that the envelop is interpreted as a density.
"""
struct Envelop
    lines::Vector{Line}
    cutpoints::Vector{Float64}
    weights::Vector{Tuple{Float64, Float64}}
    secant_knots::Vector{Float64}
    secant_values::Vector{Float64}

    Envelop(
        lines::Vector{Line},
        support::Tuple{Float64, Float64},
        secant_knots::Vector{Float64} = Float64[],
        secant_values::Vector{Float64} = Float64[]
    ) = begin
        @assert length(secant_knots) == length(secant_values)
        #  println("lines: $lines")
        @assert issorted([l.slope for l in lines], rev = true) "line slopes must be decreasing"
        intersections = [intersection(lines[i], lines[i + 1]) for i in 1:(length(lines) - 1)]
        #  println("intersections: $intersections")
        cutpoints = [support[1]; intersections; support[2]]
        #  println("cutpoints: $cutpoints")
        weights = [exp_integral(l, cutpoints[i], cutpoints[i + 1]) for (i, l) in enumerate(lines)]
        # println("weights: $weights")
        new(lines, cutpoints, weights, secant_knots, secant_values)
    end
end
Base.length(e::Envelop) = length(e.lines)

function eval_secants(e::Envelop, x::Float64; log::Bool=false)
    pos = searchsortedfirst(e.secant_knots, x)
    if pos == 1 || pos == length(e) + 1
        out =- Inf
    else
        x1 = e.secant_knots[pos - 1]
        y1 = e.secant_values[pos - 1]
        x2 = e.secant_knots[pos]
        y2 = e.secant_values[pos]
        out = (y2 - y1) / (x2 - x1) * (x - x1) + y1
    end
    !log && (out = exp(out))
    return out
end

# struct SecantBuffer
#     secants::Vector{Buffer}
# end
# len(s::SecantBuffer) = length(s.secants)

# function eval_secant_buffer(s::SecatBuffer, x::Float64)
#     for sec in s.secants
#         (sec.x1 ≥ x) && (return eval_secant(sec, x))
#     end
#     error("point $x outside secant range $(s.secants[1].x1, s.secants[end].x2)")
# end

"""
    add_segment!(e::Envelop, l::Line)
Adds a new line segment to an envelop based on the value of its slope (slopes must be decreasing
always in the envelop). The cutpoints are automatically determined by intersecting the line with
the adjacent lines.
"""
function add_segment!(e::Envelop, l::Line)
    #  println("@add_segment")
    #  println("|-- adding segment $l")
    # printstyled("|-- envelop: $e\n", color=:blue)
    # Find the position in sorted array with binary search
    pos = searchsortedfirst([-line.slope for line in e.lines], -l.slope)
    # printstyled("|-- adding $l to position $pos\n", color=:green)
    # remember cutpoints has the support at the beginning and end!
    # Find the new cutpoints
    if pos == 1
        new_cut = intersection(l, e.lines[pos])  # smaller than all values
        # printstyled("|-- adding cut $(new_cut)\n", color=:yellow)
        # Insert in second position, first one is the support bound
        # 1. measure weight from support to new first cut
        # 2. measure weight from new cut to previous first cut after support
        ws1 = exp_integral(l, e.cutpoints[1], new_cut)  # support to new
        ws2 = exp_integral(e.lines[1], new_cut, e.cutpoints[2])  # new to next
        insert!(e.cutpoints, 2, new_cut)
        splice!(e.weights, 1, [ws1, ws2])
        # printstyled("|-- new weights: ($ws1, $ws2)\n", color=:green)
        # printstyled("|-- new envelop: $e\n", color=:blue)
        # error("early exit")
        # must update weights for first two segments only
    elseif pos == length(e) + 1  # greater than all values
        new_cut = intersection(e.lines[end], l)
        # printstyled("|-- adding cut $(new_cut)\n", color=:yellow)
        # 1. measure weight from new last cut to support
        # 2. measure weight from new cut to previous first cut after support
        ws1 = exp_integral(e.lines[end], e.cutpoints[end - 1], new_cut)
        ws2 = exp_integral(l, new_cut, e.cutpoints[end])
        insert!(e.cutpoints, length(e) + 1, new_cut)  # just before support end
        splice!(e.weights, length(e), [ws1, ws2])
        # printstyled("|-- new weights: ($ws1, $ws2)\n", color=:green)
        # printstyled("|-- new envelop: $e\n", color=:blue)
        # error("early exit")
    else
        (l.slope == e.lines[pos].slope) && return -1 # can't add segment
        new_cut1 = intersection(e.lines[pos - 1], l)
        new_cut2 = intersection(l, e.lines[pos])
        supp = e.cutpoints[1], e.cutpoints[end]
        @assert supp[1] ≤ new_cut1 ≤ new_cut2 "new_cut1 ≤ supp[2]: $new_cut1, new_cut2: $new_cut2, supp: $supp"
        # println("Current cuts: $(e.cutpoints)")
        # printstyled("|-- adding cuts $(new_cut1) and $(new_cut2) from tangent $(l.tangent_point)\n", color=:yellow)
        # must update three weights!
        ws1 = exp_integral(e.lines[pos - 1], e.cutpoints[pos - 1], new_cut1)
        ws2 = exp_integral(l, new_cut1, new_cut2)
        ws3 = exp_integral(e.lines[pos], new_cut2, e.cutpoints[pos + 1])
        # printstyled("|-- new weights: ($ws1, $ws2, $ws3)\n", color=:green)
        splice!(e.cutpoints, pos, [new_cut1, new_cut2])
        splice!(e.weights, (pos - 1):pos, [ws1, ws2, ws3])
        # printstyled("|-- new envelop: $e\n", color=:blue)
        # error("early exit")
        # @assert issorted(e.cutpoints)  "incompatible line: resulting intersection points aren't sorted"
    end
    # Insert the new line
    insert!(e.lines, pos, l)
    return pos
    # e.weights = [exp_integral(line, e.cutpoints[i], e.cutpoints[i + 1]) for (i, line) in enumerate(e.lines)]
end

"""
    sample_envelop(rng::AbstractRNG, e::Envelop)
    sample_envelop(e::Envelop)
Samples an element from the density defined by the envelop `e` with it's exponential weights.
See [`Envelop`](@Envelop) for details.
"""
function sample_envelop(rng::AbstractRNG, e::Envelop)
    # Randomly select a segment from the weights
    # scale each weight accordingly for numerical stability
    max_scale = maximum(s for (_, s) in e.weights)
    rel_wts = weights([w * exp(s - max_scale) for (w, s) in e.weights])
    i = sample(rng, 1:length(e), rel_wts)

    # use the inverse CDF method to sample within segment
    #  println("Sampling $i from current weights $(e.weights)")
    # println("Sampling with weights $rel_wts, max_scale: $max_scale, x1: $x1")
    l = e.lines[i]
    denom, scale = e.weights[i]
    a, b = e.lines[i].slope, e.lines[i].intercept
    x1, x2 = e.cutpoints[i:(i + 1)]
    u = rand(rng)
    if a == 0.0
        out = x1 + u * (x2 - x1)
    else
        # out = (log(a * U * denom + exp(eval_line(e.lines[i], x1) - scale)) - b + scale) / a
        # out = (
        #     log(exp(-b + scale) * u * denom * a + exp(eval_line(l, x1) - scale))
        #     - b + scale
        # ) / a
        out = (
            log(u * e.weights[i][1] * a + exp(a * x1 + b - scale))
            -b + scale
        ) / a

    end

    #  println("Sampled $out from segment [$(e.cutpoints[i]), $(e.cutpoints[i+1])]")
    return out
end

function sample_envelop(e::Envelop)
    sample_envelop(Random.GLOBAL_RNG, e)
end

"""
    eval_envelop(e::Envelop, x::Float64)
Eval point a point `x` in the piecewise linear function defined by `e`. Necessary for evaluating
the density assigned to the point `x`.
"""
function eval_envelop(e::Envelop, x::Float64; log::Bool = false)
    # searchsortedfirst is the proper method for and ordered list
    pos = searchsortedfirst(e.cutpoints, x)
    if pos == 1 || pos == length(e.cutpoints) + 1
        return log ? - Inf : 0.0
    else
        pos = searchsortedfirst(e.cutpoints, x)
        logv = eval_line(e.lines[pos - 1], x)
        return log ? logv : exp(logv)
    end
end

# --------------------------------

"""
    RejectionSampler(f::Function, support::Tuple{Float64, Float64}[ ,δ::Float64])
    RejectionSampler(f::Function, support::Tuple{Float64, Float64}, init::Tuple{Float64, Float64})
An adaptive rejection sampler to obtain iid samples from a logconcave function `f`, supported in the
domain `support` = (support[1], support[2]). To create the object, two initial points `init = init[1], init[2]`
such that `logf'(init[1]) > 0` and `logf'(init[2]) < 0` are necessary. If they are not provided, the constructor
will perform a greedy search based on `δ`.

The argument `support` must be of the form `(-Inf, Inf), (-Inf, a), (b, Inf), (a,b)`, and it represent the
interval in which f has positive value, and zero elsewhere.

## Keyword arguments
- `max_segments::Int = 10` : max size of envelop, the rejection-rate is usually slow with a small number of segments
- `max_failed_factor::Float64 = 0.001`: level at which throw an error if one single sample has a rejection rate
    exceeding this value
"""
struct RejectionSampler
    objective::Function
    envelop::Envelop
    max_segments::Int
    max_failed_rate::Float64
    support::Tuple{Float64, Float64}
    use_secants::Bool
    # Constructor when initial points are provided
    RejectionSampler(
        f::Function,
        support::Tuple{Float64, Float64},
        init::Tuple{Union{Float64, Nothing}, Union{Float64, Nothing}};
        apply_log::Bool = true,
        max_segments::Int = 25,
        max_failed_rate::Float64 = 0.001,
        autograd::Bool = true,
        use_secants::Bool = false
    ) = begin
        # check support
        logf = apply_log ? (x -> log(f(x))) : f
        objective = autograd ? (x -> (logf(x), derivative(logf, x))) : logf
        @assert support[1] < support[2] "invalid support, not an interval"

        # check initial cutpoints
        x1, x2 = init
        isnothing(x1) && (x1 = support[1])
        isnothing(x2) && (x2 = support[2])

        @assert !isinf(x1) "initial left point can't be infinite if support unbounded below"
        @assert !isinf(x2) "initial right point can't be infinite if support unbounded above"
        @assert x1 < x2 "Initial points ($x1, $x2) must be ordered"

        # validate slopes
        lines = Line[]
        secant_knots = Float64[]
        secant_values = Float64[]

        if x1 != support[1]
            # println(objective(x1))
            y1, a1 = objective(x1)
            @assert !isinf(support[1]) || a1 > 0.0 "slope must be positive in left init when unbounded below"
            b1 = y1 - a1 * x1
            line1 = Line(a1, b1)
            push!(lines, line1)
            if use_secants
                push!(secant_knots, x1)
                push!(secant_values, y1)
            end
        end
        if x2 != support[2]
            y2, a2 = objective(x2)
            @assert !isinf(support[2]) || a2 < 0.0 "slope must be positive in right init when unbounded above"
            b2 = y2 - a2 * x2
            line2 = Line(a2, b2)
            push!(lines, line2)
            if use_secants
                push!(secant_knots, x2)
                push!(secant_values, y2)
            end
        end

        if isempty(lines)
             # special case when starting with support only bounds
            @assert !isinf(support[1]) && !isinf(support[2]) "when not providing initial points, support must be bounded"
            xmid = 0.5(support[1] + support[2])
            ymid, amid = objective(xmid)
            bmid = ymid - amid * xmid
            line12 = Line(amid, bmid)
            push!(lines, line12)
            if use_secants
                push!(secant_knots, xmid)
                push!(secant_values, ymid)
            end
        end
        envelop = Envelop(lines, support, secant_knots, secant_values)
        # println("initial envelop: $envelop")
        new(
            objective,
            envelop,
            max_segments,
            max_failed_rate,
            support,
            use_secants
        )
    end

    # Constructor for greedy search of starting points
    RejectionSampler(
            f::Function,
            support::Tuple{Float64, Float64},
            δ::Float64;
            search_range::Tuple{Float64, Float64} = (-10.0,10.0),
            apply_log::Bool = true,
            grad::Union{Function, Nothing} = nothing,
            autograd::Bool = true,
            kwargs...
    ) = begin
        logf = apply_log ? (x -> log(f(x))) : f
        objective = autograd ? (x -> (logf(x), derivative(logf, x))) : logf
        grid_lims = max(search_range[1], support[1]), min(search_range[2], support[2])
        grid = grid_lims[1]:δ:grid_lims[2]
        i1, i2 = findfirst(x -> objective(x)[2] > 0.0, grid), findfirst(objective -> f(x)[2] < 0., grid)
        @assert !isnothing(i1) &&  !isnothing(i2) "couldn't find initial points, please provide them or change `search_range`"
        x1, x2 = grid[i1], grid[i2]
        # pass autograd false since we already checked it
        RejectionSampler(objective, support, (x1, x2); apply_log=false, autograd=false, kwargs...)
    end
end

"""
    run_sampler!(rng::AbstractRNG, sampler::RejectionSampler, n::Int)
    run_sampler!(sampler::RejectionSampler, n::Int)
It draws `n` iid samples of the objective function of `sampler`, and at each iteration it adapts the envelop
of `sampler` by adding new segments to its envelop.
"""
function run_sampler!(rng::AbstractRNG, s::RejectionSampler, n::Int)
    i = 0
    failed, max_failed = 0, trunc(Int, n / s.max_failed_rate)
    out = zeros(n)
    env = s.envelop
    while i < n
        c = sample_envelop(rng, env)
        u = rand(rng)
        majorant = eval_envelop(env, c, log=true)
        if s.use_secants # squeeze test
            if u < exp(eval_secants(env, c, log=true) - majorant)
                i += 1
                out[i] = c
                continue
            end
        end
        failed += 1
        @assert failed < max_failed "max_failed_factor reached"
        y, a = s.objective(c)
        if u < exp(y - majorant) # evalf
            i += 1
            out[i] = c
            (i == n) && break
            # println(c)
            # error("early stop")
        end
        # add to segments
        if length(env) <= s.max_segments
            b = y - a * c
            pos = add_segment!(env, Line(a, b))
            if pos > 0 && s.use_secants
                insert!(env.secant_knots, pos, c)
                insert!(env.secant_values, pos, y)
            end
        end
    end
    out
end

function run_sampler!(s::RejectionSampler, n::Int)
    run_sampler!(Random.GLOBAL_RNG, s, n)
end

function eval_envelop(s::RejectionSampler, x::Float64, kwargs...)
    eval_envelop(s.envelop, x, kwargs...)
end
function eval_envelop(s::RejectionSampler, x::AbstractVector{Float64}, kwargs...)
    [eval_envelop(s.envelop, xᵢ, kwargs...) for xᵢ in x]
end
function eval_secants(s::RejectionSampler, x::Float64, kwargs...)
    eval_secants(s.envelop, x, kwargs...)
end
function eval_secants(s::RejectionSampler, x::AbstractVector{Float64}, kwargs...)
    [eval_secants(s.envelop, xᵢ, kwargs...) for xᵢ in x]
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

include("model_CloudError.jl")
include("model_ErrorAdditionalInfo.jl")
include("model_ErrorResponse.jl")
include("model_ParameterDefinitions.jl")
include("model_ParameterDefinitionsValue.jl")
include("model_PolicyDefinition.jl")
include("model_PolicyDefinitionListResult.jl")
include("model_PolicyDefinitionProperties.jl")