text:
using StatProfilerHTML
using ThreadPools
using ProgressMeter

"""
    compute_batched_pt2(ci_vector::ClusteredState{T,N,R}, cluster_ops, clustered_ham::ClusteredOperator;
        nbody=4,
        H0="Hcmf",
        E0=nothing, #pass in <0|H0|0>, or compute it
        thresh_foi=1e-8,
        prescreen=true,
        verbose=1) where {T,N,R}
"""
function compute_batched_pt2(ci_vector_in::ClusteredState{T,N,R}, cluster_ops, clustered_ham::ClusteredOperator;
        nbody=4,
        H0="Hcmf",
        E0=nothing, #pass in <0|H0|0>, or compute it
        thresh_foi=1e-9,
        prescreen=true,
        verbose=1) where {T,N,R}
    #={{{=#

    println()
    println(" |........................do batched PT2............................")
    println(" thresh_foi    :", thresh_foi   )
    println(" prescreen     :", prescreen   )
    println(" H0            :", H0   )
    println(" nbody         :", nbody   )

    e2 = zeros(T,R)

    ci_vector = deepcopy(ci_vector_in)
    clusters = ci_vector.clusters
    norms = norm(ci_vector);
    @printf(" Norms of input states:\n")
    [@printf(" %12.8f\n",i) for i in norms]
    orthonormalize!(ci_vector)

    clustered_ham_0 = extract_1body_operator(clustered_ham, op_string = H0)
    if E0 == nothing
        @printf(" %-50s", "Compute <0|H0|0>:")
        @time E0 = compute_expectation_value_parallel(ci_vector, cluster_ops, clustered_ham_0)
        #E0 = diag(E0)
        flush(stdout)
    end
    @printf(" %-50s", "Compute <0|H|0>:")
    @time Evar = compute_expectation_value_parallel(ci_vector, cluster_ops, clustered_ham)
    #Evar = diag(Evar)
    flush(stdout)

    #
    # define batches (FockConfigs present in resolvant)
    jobs = Dict{FockConfig{N},Vector{Tuple}}()
    for (fock_ket, configs_ket) in ci_vector.data
        for (ftrans, terms) in clustered_ham
            fock_x = ftrans + fock_ket

            #
            # check to make sure this fock config doesn't have negative or too many electrons in any cluster
            all(f[1] >= 0 for f in fock_x) || continue
            all(f[2] >= 0 for f in fock_x) || continue
            all(f[1] <= length(clusters[fi]) for (fi,f) in enumerate(fock_x)) || continue
            all(f[2] <= length(clusters[fi]) for (fi,f) in enumerate(fock_x)) || continue

            job_input = (terms, fock_ket, configs_ket)
            if haskey(jobs, fock_x)
                push!(jobs[fock_x], job_input)
            else
                jobs[fock_x] = [job_input]
            end

        end
    end

    #
    # prepare scratch arrays to help cut down on allocation in the threads
    jobs_vec = []
    for (fock_x, job) in jobs
        push!(jobs_vec, (fock_x, job))
    end

    scr_f = Vector{Vector{Vector{Float64}} }()
    scr_i = Vector{Vector{Vector{Int16}} }()
    scr_m = Vector{Vector{MVector{N,Int16}} }()
    nscr = 20

    scr1 = Vector{Vector{Float64}}()
    scr2 = Vector{Vector{Float64}}()
    scr3 = Vector{Vector{Float64}}()
    scr4 = Vector{Vector{Float64}}()
    tmp1 = Vector{MVector{N,Int16}}()
    tmp2 = Vector{MVector{N,Int16}}()

    e2_thread = Vector{Vector{Float64}}()
    for tid in 1:Threads.nthreads()
        push!(e2_thread, zeros(R))
        push!(scr1, zeros(1000))
        push!(scr2, zeros(1000))
        push!(scr3, zeros(1000))
        push!(scr4, zeros(1000))
        push!(tmp1, zeros(Int16,N))
        push!(tmp2, zeros(Int16,N))

        tmp = Vector{Vector{Float64}}()
        [push!(tmp, zeros(Float64,10000)) for i in 1:nscr]
        push!(scr_f, tmp)

        tmp = Vector{Vector{Int16}}()
        [push!(tmp, zeros(Int16,10000)) for i in 1:nscr]
        push!(scr_i, tmp)

        tmp = Vector{MVector{N,Int16}}()
        [push!(tmp, zeros(Int16,N)) for i in 1:nscr]
        push!(scr_m, tmp)
    end

    println(" Number of jobs:    ", length(jobs_vec))
    println(" Number of threads: ", Threads.nthreads())
    BLAS.set_num_threads(1)
    flush(stdout)

    tmp = Int(round(length(jobs_vec)/100))
    verbose < 1 || println("   |----------------------------------------------------------------------------------------------------|")
    verbose < 1 || println("   |0%                                                                                100%|")
    verbose < 1 || print("   |")
    #@profilehtml @Threads.threads for job in jobs_vec
    t = @elapsed begin
        #@qthreads for job in jobs_vec
        #@time for job in jobs_vec

        @Threads.threads for (jobi,job) in collect(enumerate(jobs_vec))
        #for (jobi,job) in collect(enumerate(jobs_vec))
            fock_bra = job[1]
            tid = Threads.threadid()
            e2_thread[tid] .+= _pt2_job(job[2], fock_bra, cluster_ops, nbody, thresh_foi,
                                        scr_f[tid], scr_i[tid], scr_m[tid],  prescreen, verbose,
                                        ci_vector, clustered_ham_0, E0)
            if verbose > 0
                if  jobi%tmp == 0
                    print("-")
                    flush(stdout)
                end
            end
        end
    end
    verbose < 1 || println("|")
    flush(stdout)

    @printf(" Time spent computing E2 %12.1f (s)\n",t)
    e2 = sum(e2_thread)

    #BLAS.set_num_threads(Threads.nthreads())

    @printf(" %5s %12s %12s\n", "Root", "E(0)", "E(2)")
    for r in 1:R
        @printf(" %5s %12.8f %12.8f\n",r, Evar[r], Evar[r] + e2[r])
    end
    println(" ..................................................................|")

    return e2
end
#=}}}=#

function _pt2_job(job, fock_x, cluster_ops, nbody, thresh,
                  scr_f, scr_i, scr_m, prescreen, verbose,
                  ci_vector::ClusteredState{T,N,R}, clustered_ham_0, E0) where {T,N,R}
    #={{{=#

    sig = ClusteredState(ci_vector.clusters, T=T, R=R)
    add_fockconfig!(sig, fock_x)

    for jobi in job

        terms, fock_ket, configs_ket = jobi

        for term in terms

            length(term.clusters) <= nbody || continue

            for (config_ket, coeff_ket) in configs_ket

                #term isa ClusteredTerm2B || continue
                #if (length(sig[fock_x]) > 0) && (term isa ClusteredTerm4B)
                if (term isa ClusteredTerm4B)
                    #println("1: ", length(sig[fock_x]))
                    #@btime contract_matvec_thread($term, $cluster_ops, $fock_x, $fock_ket, $config_ket, $coeff_ket,
                    #                       $sig[$fock_x],
                    #                       $scr_f, $scr_i, $scr_m,
                    #                       thresh=$thresh, prescreen=$prescreen)
                    #contract_matvec_thread(term, cluster_ops, fock_x, fock_ket, config_ket, coeff_ket,
                    #                   sig[fock_x],
                    #                   scr_f, scr_i, scr_m,
                    #                   thresh=thresh, prescreen=prescreen)
                    #println("2: ", length(sig[fock_x]))
                    #@profilehtml contract_matvec_thread(term, cluster_ops, fock_x, fock_ket, config_ket, coeff_ket,
                    #                   sig[fock_x],
                    #                   scr_f, scr_i, scr_m,
                    #                   thresh=thresh, prescreen=prescreen)
                    #error("we good?")
                end

                contract_matvec_thread(term, cluster_ops, fock_x, fock_ket, config_ket, coeff_ket,
                                       sig[fock_x],
                                       scr_f, scr_i, scr_m,
                                       thresh=thresh, prescreen=prescreen)
                #if term isa ClusteredTerm3B
                    #@code_warntype contract_matvec_thread(term, cluster_ops, fock_x, fock_ket, config_ket,
                    #                                coeff_ket,  sig[fock_x],scr1, scr2, thresh=thresh)
                    #@btime contract_matvec_thread($term, $cluster_ops, $fock_x, $fock_ket, $config_ket,
                    #                        $coeff_ket, $sig[$fock_x], $scr_f, $scr_i, $scr_m, thresh=$thresh)
                #end
            end
        end
    end

    #@btime project_out!($sig, $ci_vector)
    #@time project_out!(sig, ci_vector)

    project_out!(sig, ci_vector)
    #verbose > 0 || println(" Fock(X): ", fock_x)
    #verbose > 1 || println(" Length of FOIS vector: ", length(sig))
    #verbose > 1 || println(" Compute diagonal")

    nx = length(sig)

    Hd = scr_f[9]
    resize!(Hd, nx)
    Hd = reshape2(Hd, (nx, 1))

    #Hd = compute_diagonal(sig, cluster_ops, clustered_ham_0)
    fill!(Hd,0.0)
    compute_diagonal!(Hd, sig, cluster_ops, clustered_ham_0)
    #@btime compute_diagonal!($Hd, $sig, $cluster_ops, $clustered_ham_0)

    sig_v = scr_f[10]
    resize!(sig_v, nx*R)
    sig_v = reshape2(sig_v, size(sig))
    fill!(sig_v,0.0)

    get_vectors!(sig_v, sig)

    #sig_v = get_vectors(sig)

    e2 = zeros(T,R)

    _sum_pt2(sig_v, e2, Hd, E0, R)

    return e2
end
#=}}}=#

function _sum_pt2(sig_v, e2, Hd, E0, R)
    # put into a function to let compiler specialize
    nx = length(Hd)
    sig_vx = 0.0
    Hdx = 0.0
    @inbounds for x in 1:nx
        Hdx = Hd[x]
        @simd for r in 1:R
            sig_vx = sig_v[x,r]
            e2[r] += sig_vx*sig_vx / (E0[r] - Hdx)
            #verbose > 0 || @printf(" %5s %12.8f\n",r, e2[r])
        end
    end
end

"""
    open_matvec_thread2(ci_vector::ClusteredState{T,N,R}, cluster_ops, clustered_ham; thresh=1e-9, nbody=4) where {T,N,R}

Compute the action of the Hamiltonian on a tpsci state vector. Open here, means that we access the full FOIS
(restricted only by thresh), instead of the action of H on v within a subspace of configurations.
This is essentially used for computing a PT correction outside of the subspace, or used for searching in TPSCI.

This parallellizes over FockConfigs in the output state, so it's not the most fine-grained, but it avoids data races in
filling the final vector
"""
function open_matvec_thread2(ci_vector::ClusteredState{T,N,R}, cluster_ops, clustered_ham;
                             thresh=1e-9,
                             prescreen=true,
                             nbody=4) where {T,N,R}
#={{{=#
    println(" In open_matvec_thread2")
    sig = deepcopy(ci_vector)
    zero!(sig)
    clusters = ci_vector.clusters
    jobs = Dict{FockConfig{N},Vector{Tuple}}()
    #sig = ClusteredState(clusters)
    #sig = OrderedDict{FockConfig{N}, OrderedDict{NTuple{N,Int16}, MVector{T} }}()

    @printf(" %-50s", "Setup threaded jobs: ")
    @time for (fock_ket, configs_ket) in ci_vector.data
        for (ftrans, terms) in clustered_ham
            fock_bra = ftrans + fock_ket

            #
            # check to make sure this fock config doesn't have negative or too many electrons in any cluster
            all(f[1] >= 0 for f in fock_bra) || continue
            all(f[2] >= 0 for f in fock_bra) || continue
            all(f[1] <= length(clusters[fi]) for (fi,f) in enumerate(fock_bra)) || continue
            all(f[2] <= length(clusters[fi]) for (fi,f) in enumerate(fock_bra)) || continue

            job_input = (terms, fock_ket, configs_ket)
            if haskey(jobs, fock_bra)
                push!(jobs[fock_bra], job_input)
            else
                jobs[fock_bra] = [job_input]
            end

        end
    end

    jobs_vec = []
    for (fock_bra, job) in jobs
        push!(jobs_vec, (fock_bra, job))
    end

    scr_f = Vector{Vector{Vector{Float64}} }()
    scr_i = Vector{Vector{Vector{Int16}} }()
    scr_m = Vector{Vector{MVector{N,Int16}} }()
    nscr = 20

    scr1 = Vector{Vector{Float64}}()
    scr2 = Vector{Vector{Float64}}()
    scr3 = Vector{Vector{Float64}}()
    scr4 = Vector{Vector{Float64}}()
    tmp1 = Vector{MVector{N,Int16}}()
    tmp2 = Vector{MVector{N,Int16}}()

    jobs_out = Vector{ClusteredState{T,N,R}}()
    for tid in 1:Threads.nthreads()
        push!(jobs_out, ClusteredState(clusters, T=T, R=R))
        push!(scr1, zeros(1000))
        push!(scr2, zeros(1000))
        push!(scr3, zeros(1000))
        push!(scr4, zeros(1000))
        push!(tmp1, zeros(Int16,N))
        push!(tmp2, zeros(Int16,N))

        tmp = Vector{Vector{Float64}}()
        [push!(tmp, zeros(Float64,10000)) for i in 1:nscr]
        push!(scr_f, tmp)

        tmp = Vector{Vector{Int16}}()
        [push!(tmp, zeros(Int16,10000)) for i in 1:nscr]
        push!(scr_i, tmp)

        tmp = Vector{MVector{N,Int16}}()
        [push!(tmp, zeros(Int16,N)) for i in 1:nscr]
        push!(scr_m, tmp)
    end

    println(" Number of jobs:    ", length(jobs_vec))
    println(" Number of threads: ", Threads.nthreads())
    BLAS.set_num_threads(1)
    #Threads.@threads for job in jobs_vec

    flush(stdout)

    #@time for job in jobs_vec
    #@qthreads for job in jobs_vec
    @printf(" %-50s", "Compute matrix-vector: ")
    @time @Threads.threads for job in jobs_vec
        fock_bra = job[1]
        tid = Threads.threadid()
        _open_matvec_thread2_job(job[2], fock_bra, cluster_ops, nbody, thresh,
                                 jobs_out[tid], scr_f[tid], scr_i[tid], scr_m[tid], prescreen)
    end
    flush(stdout)

    @printf(" %-50s", "Now collect thread results: ")
    flush(stdout)
    @time for threadid in 1:Threads.nthreads()
        add!(sig, jobs_out[threadid])
    end

    #BLAS.set_num_threads(Threads.nthreads())
    return sig
end
#=}}}=#

function open_matvec_serial2(ci_vector::ClusteredState{T,N,R}, cluster_ops, clustered_ham;
                             thresh=1e-9,
                             prescreen=true,
                             nbody=4) where {T,N,R}
#={{{=#
    println(" In open_matvec_serial2")
    sig = deepcopy(ci_vector)
    zero!(sig)
    clusters = ci_vector.clusters
    jobs = Dict{FockConfig{N},Vector{Tuple}}()

    @printf(" %-50s", "Setup threaded jobs: ")
    @time for (fock_ket, configs_ket) in ci_vector.data
        for (ftrans, terms) in clustered_ham
            fock_bra = ftrans + fock_ket

            #
            # check to make sure this fock config doesn't have negative or too many electrons in any cluster
            all(f[1] >= 0 for f in fock_bra) || continue
            all(f[2] >= 0 for f in fock_bra) || continue
            all(f[1] <= length(clusters[fi]) for (fi,f) in enumerate(fock_bra)) || continue
            all(f[2] <= length(clusters[fi]) for (fi,f) in enumerate(fock_bra)) || continue

            job_input = (terms, fock_ket, configs_ket)
            if haskey(jobs, fock_bra)
                push!(jobs[fock_bra], job_input)
            else
                jobs[fock_bra] = [job_input]
            end

        end
    end

    jobs_vec = []
    for (fock_bra, job) in jobs
        push!(jobs_vec, (fock_bra, job))
    end

    scr_f = Vector{Vector{Vector{Float64}} }()
    scr_i = Vector{Vector{Vector{Int16}} }()
    scr_m = Vector{Vector{MVector{N,Int16}} }()
    nscr = 20

    scr1 = Vector{Vector{Float64}}()
    scr2 = Vector{Vector{Float64}}()
    scr3 = Vector{Vector{Float64}}()
    scr4 = Vector{Vector{Float64}}()
    tmp1 = Vector{MVector{N,Int16}}()
    tmp2 = Vector{MVector{N,Int16}}()

    jobs_out = Vector{ClusteredState{T,N,R}}()
    for tid in 1:1
        push!(jobs_out, ClusteredState(clusters, T=T, R=R))
        push!(scr1, zeros(1000))
        push!(scr2, zeros(1000))
        push!(scr3, zeros(1000))
        push!(scr4, zeros(1000))
        push!(tmp1, zeros(Int16,N))
        push!(tmp2, zeros(Int16,N))

        tmp = Vector{Vector{Float64}}()
        [push!(tmp, zeros(Float64,10000)) for i in 1:nscr]
        push!(scr_f, tmp)

        tmp = Vector{Vector{Int16}}()
        [push!(tmp, zeros(Int16,10000)) for i in 1:nscr]
        push!(scr_i, tmp)

        tmp = Vector{MVector{N,Int16}}()
        [push!(tmp, zeros(Int16,N)) for i in 1:nscr]
        push!(scr_m, tmp)
    end

    println(" Number of jobs:    ", length(jobs_vec))
    flush(stdout)
    @printf(" %-50s", "Compute matrix-vector: ")
    @time for job in jobs_vec
        fock_bra = job[1]
        tid = 1
        _open_matvec_thread2_job(job[2], fock_bra, cluster_ops, nbody, thresh,
                                 jobs_out[tid], scr_f[tid], scr_i[tid], scr_m[tid],prescreen)
    end
    flush(stdout)

    flush(stdout)
    @printf(" %-50s", "Now collect thread results: ")
    @time for threadid in 1:1
        add!(sig, jobs_out[threadid])
    end

    return sig
end
#=}}}=#

function _open_matvec_thread2_job(job, fock_bra, cluster_ops, nbody, thresh, sig, scr_f, scr_i, scr_m, prescreen)
#={{{=#

    haskey(sig, fock_bra) || add_fockconfig!(sig, fock_bra)
    #sizehint!(sig[fock_bra],1000)

    for jobi in job

        terms, fock_ket, configs_ket = jobi

        for term in terms

            length(term.clusters) <= nbody || continue

            for (config_ket, coeff_ket) in configs_ket

                #term isa ClusteredTerm2B || continue

                contract_matvec_thread(term, cluster_ops, fock_bra, fock_ket, config_ket, coeff_ket, sig[fock_bra],
                                       scr_f, scr_i, scr_m, thresh=thresh, prescreen=prescreen)
                if term isa ClusteredTerm3B
                    #@code_warntype contract_matvec_thread(term, cluster_ops, fock_bra, fock_ket, config_ket, coeff_ket,
                    #                              sig[fock_bra],scr1, scr2, thresh=thresh)
                    #@btime contract_matvec_thread($term, $cluster_ops, $fock_bra, $fock_ket, $config_ket, $coeff_ket,
                    #                              $sig[$fock_bra], $scr_f, $scr_i, $scr_m, thresh=$thresh)
                end
            end
        end
    end
    return
end
#=}}}=#

reshape2(a, dims) = invoke(Base._reshape, Tuple{AbstractArray,typeof(dims)}, a, dims)

"""
    contract_matvec_thread(   term::ClusteredTerm1B,
                                    cluster_ops::Vector{ClusterOps},
                                    fock_bra::FockConfig{N},
                                    fock_ket::FockConfig{N}, conf_ket::ClusterConfig{N}, coef_ket::MVector{R,T},
                                    sig,
                                    scr_f::Vector{Vector{Float64}},
                                    scr_i::Vector{Vector{Int16}},
                                    scr_m::Vector{MVector{N,Int16}};
                                    thresh=1e-9) where {T,R,N}
"""
function contract_matvec_thread(   term::ClusteredTerm1B,
                                    cluster_ops::Vector{ClusterOps},
                                    fock_bra::FockConfig{N},
                                    fock_ket::FockConfig{N}, conf_ket::ClusterConfig{N}, coef_ket::MVector{R,T},
                                    sig,
                                    scr_f::Vector{Vector{Float64}},
                                    scr_i::Vector{Vector{Int16}},
                                    scr_m::Vector{MVector{N,Int16}};
                                    thresh=1e-9, prescreen=true) where {T,R,N}
#={{{=#
    c1 = term.clusters[1]

    #
    # determine sign from rearranging clusters if odd number of operators
    state_sign = compute_terms_state_sign(term, fock_ket)

    #
    # <:|p'|J> h(pq) <:|q|L>
    #
    # <:|p'|J> h(pq) <:|q|L>
    g1a = cluster_ops[c1.idx][term.ops[1]]

    # the following type declarations are needed to prevent allocations. This ultimately is due to the fact
    # that the ClusterOps type isn't formed in an ideal concrete manner. We just currently use Array.
    # This should be cleaned up eventually, preferably with a distinction between contracted, and uncontracted
    # operators.
    g1::Array{Float64,2} = g1a[(fock_bra[c1.idx],fock_ket[c1.idx])]

    @views gamma1 = g1[:,conf_ket[c1.idx]]

    scr1 = scr_f[1]
    resize!(scr1, size(gamma1,1))

    scr1 .= gamma1 .* state_sign
    newI = 1:size(scr1,1)
    #new_coeffs = scr1
    #newI = 1:size(new_coeffs,1)

    _collect_significant_thread!(sig, conf_ket, scr1, coef_ket, c1.idx,  newI,  thresh, scr_m[1])
    #_collect_significant_thread!(sig, conf_ket, new_coeffs, coef_ket, c1.idx,  newI,  thresh)

    return
end
#=}}}=#

"""
    contract_matvec_thread(   term::ClusteredTerm2B,
                                    cluster_ops::Vector{ClusterOps},
                                    fock_bra::FockConfig{N},
                                    fock_ket::FockConfig{N}, conf_ket::ClusterConfig{N}, coef_ket::MVector{R,T},
                                    sig,
                                    scr_f::Vector{Vector{Float64}},
                                    scr_i::Vector{Vector{Int16}},
                                    scr_m::Vector{MVector{N,Int16}};
                                    thresh=1e-9) where {T,R,N}
"""
function contract_matvec_thread(   term::ClusteredTerm2B,
                                    cluster_ops::Vector{ClusterOps},
                                    fock_bra::FockConfig{N},
                                    fock_ket::FockConfig{N}, conf_ket::ClusterConfig{N}, coef_ket::MVector{R,T},
                                    sig,
                                    scr_f::Vector{Vector{Float64}},
                                    scr_i::Vector{Vector{Int16}},
                                    scr_m::Vector{MVector{N,Int16}};
                                    thresh=1e-9, prescreen=true) where {T,R,N}
#={{{=#

    c1 = term.clusters[1]
    c2 = term.clusters[2]

    #
    # determine sign from rearranging clusters if odd number of operators
    state_sign = compute_terms_state_sign(term, fock_ket)

    #
    # <:|p'|J> h(pq) <:|q|L>
    g1a = cluster_ops[c1.idx][term.ops[1]]
    g2a = cluster_ops[c2.idx][term.ops[2]]

    # the following type declarations are needed to prevent allocations. This ultimately is due to the fact
    # that the ClusterOps type isn't formed in an ideal concrete manner. We just currently use Array.
    # This should be cleaned up eventually, preferably with a distinction between contracted, and uncontracted
    # operators.

    haskey(g1a, (fock_bra[c1.idx],fock_ket[c1.idx])) || return
    haskey(g2a, (fock_bra[c2.idx],fock_ket[c2.idx])) || return

    g1::Array{Float64,3} = g1a[(fock_bra[c1.idx],fock_ket[c1.idx])]
    g2::Array{Float64,3} = g2a[(fock_bra[c2.idx],fock_ket[c2.idx])]

    @views gamma1 = g1[:,:,conf_ket[c1.idx]]
    @views gamma2 = g2[:,:,conf_ket[c2.idx]]

    scr1 = scr_f[1]
    scr2 = scr_f[2]

    resize!(scr1, size(term.ints,2) * size(gamma1,2))
    resize!(scr2, size(gamma1,2) * size(gamma2,2))

    scr1 = reshape2(scr1, (size(term.ints,2), size(gamma1,2)))
    scr2 = reshape2(scr2, (size(gamma1,2), size(gamma2,2)))

    mul!(scr1, term.ints', gamma1)
    mul!(scr2, scr1', gamma2)
    new_coeffs = scr2

    newI = 1:size(new_coeffs,1)
    newJ = 1:size(new_coeffs,2)

    _collect_significant_thread!(sig, conf_ket, new_coeffs, coef_ket, c1.idx, c2.idx, newI, newJ, thresh, scr_m[1], state_sign)
    #@btime _collect_significant_thread!($sig, $conf_ket, $new_coeffs, $coef_ket, $c1.idx, $c2.idx, $newI, $newJ, $thresh, $scr3)

    return
end
#=}}}=#

"""
    contract_matvec_thread(   term::ClusteredTerm3B,
                                    cluster_ops::Vector{ClusterOps},
                                    fock_bra::FockConfig{N},
                                    fock_ket::FockConfig{N}, conf_ket::ClusterConfig{N}, coef_ket::MVector{R,T},
                                    sig,
                                    scr_f::Vector{Vector{Float64}},
                                    scr_i::Vector{Vector{Int16}},
                                    scr_m::Vector{MVector{N,Int16}};
                                    thresh=1e-9, prescreen=true) where {T,R,N}

This version should only use M^2N^2 storage, and n^4 scaling n={MN}
"""
function contract_matvec_thread(   term::ClusteredTerm3B,
                                    cluster_ops::Vector{ClusterOps},
                                    fock_bra::FockConfig{N},
                                    fock_ket::FockConfig{N}, conf_ket::ClusterConfig{N}, coef_ket::MVector{R,T},
                                    sig,
                                    scr_f::Vector{Vector{Float64}},
                                    scr_i::Vector{Vector{Int16}},
                                    scr_m::Vector{MVector{N,Int16}};
                                    thresh=1e-9, prescreen=true) where {T,R,N}
#={{{=#
    c1 = term.clusters[1]
    c2 = term.clusters[2]
    c3 = term.clusters[3]

    #
    # determine sign from rearranging clusters if odd number of operators
    state_sign = compute_terms_state_sign(term, fock_ket)

    #
    # h(pqr) <I|p'|_> <J|q|_> <K|r|_>
    #
    # X(p,J,K) = h(pqr) <J|q|_> <K|r|_>
    #
    #
    haskey(cluster_ops[c1.idx][term.ops[1]],  (fock_bra[c1.idx],fock_ket[c1.idx])) || return
    haskey(cluster_ops[c2.idx][term.ops[2]],  (fock_bra[c2.idx],fock_ket[c2.idx])) || return
    haskey(cluster_ops[c3.idx][term.ops[3]],  (fock_bra[c3.idx],fock_ket[c3.idx])) || return

    g1::Array{Float64,3} = cluster_ops[c1.idx][term.ops[1]][(fock_bra[c1.idx],fock_ket[c1.idx])]
    g2::Array{Float64,3} = cluster_ops[c2.idx][term.ops[2]][(fock_bra[c2.idx],fock_ket[c2.idx])]
    g3::Array{Float64,3} = cluster_ops[c3.idx][term.ops[3]][(fock_bra[c3.idx],fock_ket[c3.idx])]
    @views gamma1 = g1[:,:,conf_ket[c1.idx]]
    @views gamma2 = g2[:,:,conf_ket[c2.idx]]
    @views gamma3 = g3[:,:,conf_ket[c3.idx]]

    scr1 = scr_f[1]
    scr2 = scr_f[2]
    scr3 = scr_f[3]
    tmp1 = scr_m[1]
    tmp2 = scr_m[2]

    newI = UnitRange{Int16}(1,size(g1,2))
    newJ = UnitRange{Int16}(1,size(g2,2))
    newK = UnitRange{Int16}(1,size(g3,2))

    if prescreen
        up_bound = upper_bound_thread(term.ints, gamma1, gamma2, gamma3,
                              scr1, scr2, scr3,
                              c=maximum(abs.(coef_ket)))
        if up_bound < thresh
            return
        end
        #@btime newI, newJ, newK = upper_bound2($term.ints, $gamma1, $gamma2, $gamma3, $thresh, c=maximum(abs.($coef_ket)))
        newI, newJ, newK = upper_bound2_thread(term.ints, gamma1, gamma2, gamma3,
                                                     scr_f, scr_i, thresh, c=maximum(abs, coef_ket))
        length(newI) > 0 || return
        length(newJ) > 0 || return
        length(newK) > 0 || return

        if true
            gamma1 = scr_f[4]
            gamma2 = scr_f[5]
            gamma3 = scr_f[6]

            resize!(gamma1, size(g1,1)*length(newI))
            resize!(gamma2, size(g2,1)*length(newJ))
            resize!(gamma3, size(g3,1)*length(newK))

            gamma1 = reshape2(gamma1, (size(g1,1), length(newI)))
            gamma2 = reshape2(gamma2, (size(g2,1), length(newJ)))
            gamma3 = reshape2(gamma3, (size(g3,1), length(newK)))

            _fill1!(gamma1, newI, g1, conf_ket[c1.idx])
            _fill1!(gamma2, newJ, g2, conf_ket[c2.idx])

            #gamma1 .= g1[:,newI,conf_ket[c1.idx]]
            #gamma2 .= g2[:,newJ,conf_ket[c2.idx]]
            #gamma3 .= g3[:,newK,conf_ket[c3.idx]]
        else
            @views gamma1 = g1[:,newI,conf_ket[c1.idx]]
            @views gamma2 = g2[:,newJ,conf_ket[c2.idx]]
            @views gamma3 = g3[:,newK,conf_ket[c3.idx]]
        end
    end

    #
    # for K in G3
    #   scr1(p,q) = h(p,q,r) * G3(r;K)          N^3M^1
    #   scr2(p,J) = scr1(p,q) * G2(q,J)
    #   scr3(I,J) = G1(p,I)' * scr2(p,J)
    #
    #   collect(scr3(I,J))

    np = size(term.ints,1)
    nq = size(term.ints,2)
    nr = size(term.ints,3)
    nI = length(newI)
    nJ = length(newJ)
    nK = length(newK)

    resize!(scr1, np*nq)
    resize!(scr2, np*nJ)
    resize!(scr3, nI*nJ)

    scr2 = reshape2(scr2, (np,nJ))
    scr3 = reshape2(scr3, (nI,nJ))

    v = reshape2(term.ints, (np*nq,nr))

    #a = @allocated begin
    #end; if a > 0 println(a); error("here we are"); end

    tmp1 .= conf_ket.config

    for K::Int16 in newK

        @views gamma3 = g3[:,K,conf_ket[c3.idx]]

        tmp1[c3.idx] = K

        tmp_conf = ClusterConfig(SVector(tmp1))

        mul!(scr1, v, gamma3)
        mul!(scr2, reshape2(scr1, (np,nq)), gamma2)
        #@allocated mul!(scr2, reshape2(scr1, (np,nq)), gamma2) == 0 || error(" we allocated")

        if prescreen
            #println("3b")
            #@btime upper_bound_thread($gamma1, $scr2, c=maximum(abs.($coef_ket)))
            #println(size(gamma1))
            #println(size(scr2))
            #println(typeof(gamma1))
            #println(typeof(scr2))
            upper_bound_thread(gamma1, scr2, c=maximum(abs.(coef_ket))) > thresh || continue
        end

        mul!(scr3, gamma1', scr2)
        #@btime mul!($scr3, $gamma1', $scr2)
        #a = @allocated begin
        #    mul!(scr3, gamma1', scr2)
        #end
        #a == 0 || error(" we allocated: ", a)

        #@btime _collect_significant_thread!($sig, $tmp_conf, $scr3, $coef_ket, $c1.idx, $c2.idx, $newI, $newJ, $thresh, $tmp2, $state_sign)
        _collect_significant_thread!(sig, tmp_conf, scr3, coef_ket, c1.idx, c2.idx, newI, newJ, thresh, tmp2, state_sign)
    end

    return
end
#=}}}=#

"""
    contract_matvec_thread(   term::ClusteredTerm4B,
                                    cluster_ops::Vector{ClusterOps},
                                    fock_bra::FockConfig{N},
                                    fock_ket::FockConfig{N}, conf_ket::ClusterConfig{N}, coef_ket::MVector{R,T},
                                    sig,
                                    scr_f::Vector{Vector{Float64}},
                                    scr_i::Vector{Vector{Int16}},
                                    scr_m::Vector{MVector{N,Int16}};
                                    thresh=1e-9, prescreen=true) where {T,R,N}

This version should only use M^2N^2 storage, and n^5 scaling n={MN}
"""
function contract_matvec_thread(   term::ClusteredTerm4B,
                                    cluster_ops::Vector{ClusterOps},
                                    fock_bra::FockConfig{N},
                                    fock_ket::FockConfig{N}, conf_ket::ClusterConfig{N}, coef_ket::MVector{R,T},
                                    sig,
                                    scr_f::Vector{Vector{Float64}},
                                    scr_i::Vector{Vector{Int16}},
                                    scr_m::Vector{MVector{N,Int16}};
                                    thresh=1e-9, prescreen=true) where {T,R,N}
#={{{=#
    c1 = term.clusters[1]
    c2 = term.clusters[2]
    c3 = term.clusters[3]
    c4 = term.clusters[4]

    #
    # determine sign from rearranging clusters if odd number of operators
    state_sign = compute_terms_state_sign(term, fock_ket)

    #
    # h(pqr) <I|p'|_> <J|q|_> <K|r|_>
    #
    # X(p,J,K) = h(pqr) <J|q|_> <K|r|_>
    #
    #
    haskey(cluster_ops[c1.idx][term.ops[1]],  (fock_bra[c1.idx],fock_ket[c1.idx])) || return
    haskey(cluster_ops[c2.idx][term.ops[2]],  (fock_bra[c2.idx],fock_ket[c2.idx])) || return
    haskey(cluster_ops[c3.idx][term.ops[3]],  (fock_bra[c3.idx],fock_ket[c3.idx])) || return
    haskey(cluster_ops[c4.idx][term.ops[4]],  (fock_bra[c4.idx],fock_ket[c4.idx])) || return

    g1::Array{Float64,3} = cluster_ops[c1.idx][term.ops[1]][(fock_bra[c1.idx],fock_ket[c1.idx])]
    g2::Array{Float64,3} = cluster_ops[c2.idx][term.ops[2]][(fock_bra[c2.idx],fock_ket[c2.idx])]
    g3::Array{Float64,3} = cluster_ops[c3.idx][term.ops[3]][(fock_bra[c3.idx],fock_ket[c3.idx])]
    g4::Array{Float64,3} = cluster_ops[c4.idx][term.ops[4]][(fock_bra[c4.idx],fock_ket[c4.idx])]
    @views gamma1 = g1[:,:,conf_ket[c1.idx]]
    @views gamma2 = g2[:,:,conf_ket[c2.idx]]
    @views gamma3 = g3[:,:,conf_ket[c3.idx]]
    @views gamma4 = g4[:,:,conf_ket[c4.idx]]

    np = size(term.ints,1)
    nq = size(term.ints,2)
    nr = size(term.ints,3)
    ns = size(term.ints,4)

    scr1 = scr_f[1]
    scr2 = scr_f[2]
    scr3 = scr_f[3]
    scr4 = scr_f[4]
    tmp1 = scr_m[1]
    tmp2 = scr_m[2]

    newI = UnitRange{Int16}(1,size(g1,2))
    newJ = UnitRange{Int16}(1,size(g2,2))
    newK = UnitRange{Int16}(1,size(g3,2))
    newL = UnitRange{Int16}(1,size(g4,2))
    if prescreen
        up_bound = upper_bound_thread(term.ints, gamma1, gamma2, gamma3, gamma4,
                              scr1, scr2, scr3, scr4,
                              c=maximum(abs.(coef_ket)))
        if up_bound < thresh
            return
        end
        #println("a")
        #
        # screen phase 2: ignore indices for each cluster which will produce discarded terms

        #@btime newI, newJ, newK, newL = upper_bound2_thread($term.ints, $gamma1, $gamma2, $gamma3, $gamma4,
        #                                             $scr_f, $scr_i, $thresh, c=maximum(abs, $coef_ket))

        newI, newJ, newK, newL = upper_bound2_thread(term.ints, gamma1, gamma2, gamma3, gamma4,
                                                     scr_f, scr_i, thresh, c=maximum(abs, coef_ket))
        length(newI) > 0 || return
        length(newJ) > 0 || return
        length(newK) > 0 || return
        length(newL) > 0 || return
        #println("b")

        #
        # While we could simply use a view, this creates a non-contiguous array and BLAS then does allocations

        if true
            g1t = typeof(gamma1)
            g2t = typeof(gamma1)
            g3t = typeof(gamma1)
            g4t = typeof(gamma1)
            gamma1 = scr_f[5]
            gamma2 = scr_f[6]
            gamma3 = scr_f[7]
            gamma4 = scr_f[8]

            #g1t == typeof(gamma1) || println(" here: ", g1t, " vs ", typeof(gamma1))
            #g2t == typeof(gamma2) || println(" here: ", g2t, " vs ", typeof(gamma2))
            #g3t == typeof(gamma3) || println(" here: ", g3t, " vs ", typeof(gamma3))
            #g4t == typeof(gamma4) || println(" here: ", g4t, " vs ", typeof(gamma4))

            resize!(gamma1, size(g1,1)*length(newI))
            resize!(gamma2, size(g2,1)*length(newJ))
            resize!(gamma3, size(g3,1)*length(newK))
            resize!(gamma4, size(g4,1)*length(newL))

            gamma1 = reshape2(gamma1, (size(g1,1), length(newI)))
            gamma2 = reshape2(gamma2, (size(g2,1), length(newJ)))
            gamma3 = reshape2(gamma3, (size(g3,1), length(newK)))
            gamma4 = reshape2(gamma4, (size(g4,1), length(newL)))

            _fill1!(gamma1, newI, g1, conf_ket[c1.idx])
            _fill1!(gamma2, newJ, g2, conf_ket[c2.idx])
            #_fill1!(gamma3, newK, g3, conf_ket[c3.idx])
            #_fill1!(gamma4, newL, g4, conf_ket[c4.idx])

            #gamma1 .= g1[:,newI,conf_ket[c1.idx]]
            #gamma2 .= g2[:,newJ,conf_ket[c2.idx]]
            #gamma3 .= g3[:,newK,conf_ket[c3.idx]]
            #gamma4 .= g4[:,newL,conf_ket[c4.idx]]
        else
            @views gamma1 = g1[:,newI,conf_ket[c1.idx]]
            @views gamma2 = g2[:,newJ,conf_ket[c2.idx]]
            @views gamma3 = g3[:,newK,conf_ket[c3.idx]]
            @views gamma4 = g4[:,newL,conf_ket[c4.idx]]
        end
    end

    #
    #   for L in G4
    #       scr1(p,q,r) = h(p,q,r,s) * G4(s;L)          N^4M^1
    #
    #       for K in G3
    #           scr2(p,q) = scr1(p,q,r) * G3(r;K)       N^3M^2
    #           scr3(p,J) = scr2(p,q) * G2(q,J)         N^2M^3
    #           scr4(I,J) = G1'(p,I) * scr3(p,J)        N^1M^4
    #
    #           collect(scr4(I,J))

    nI = length(newI)
    nJ = length(newJ)
    nK = length(newK)
    nL = length(newL)

    resize!(scr1, np*nq*nr)
    resize!(scr2, np*nq)
    resize!(scr3, np*nJ)
    resize!(scr4, nI*nJ)

    scr3 = reshape2(scr3, (np,nJ))
    scr4 = reshape2(scr4, (nI,nJ))

    v = reshape2(term.ints, (np*nq*nr,ns))
    #@btime v = reshape2($term.ints, ($np*$nq*$nr,$ns))

    tmp1 .= conf_ket.config

    a = 0.0
    #@btime begin
    #    for L::Int16 in $newL
    #        $a += 1.0
    #    end
    #end

    for L::Int16 in newL

        @views gamma4 = g4[:,L,conf_ket[c4.idx]]
        #@btime @views gamma4 = $g4[:,$L,$conf_ket[$c4.idx]]

        tmp1[c4.idx] = L
        #@btime $tmp1[$c4.idx] = $L

        mul!(scr1, v, gamma4)
        #@btime mul!($scr1, $v, $gamma4)

        for K::Int16 in newK

            @views gamma3 = g3[:,K,conf_ket[c3.idx]]
            #@btime @views $gamma3 = $g3[:,$K,$conf_ket[$c3.idx]]

            tmp1[c3.idx] = K
            #@btime $tmp1[$c3.idx] = $K

            tmp_conf = ClusterConfig(SVector(tmp1))
            #@btime $tmp_conf = ClusterConfig(SVector($tmp1))

            mul!(scr2, reshape2(scr1, (np*nq,nr)), gamma3)
            #@btime mul!($scr2, reshape2($scr1, ($np*$nq,$nr)), $gamma3)

            #@views gamma2 = g2[:,:,conf_ket[c2.idx]]

            #
            # this is allocating :(  .... but only when threaded? actually @allocated doesn't work when threaded
            mul!(scr3, reshape2(scr2, (np,nq)), gamma2)
            #a = @allocated begin
            #    mul!(scr3, reshape2(scr2, (np,nq)), gamma2)
            #end
            #@allocated mul!(scr3, reshape2(scr2, (np,nq)), gamma2) == 0 || error("why allocate?")
            #@btime mul!($scr3, reshape2($scr2, ($np,$nq)), $gamma2)

            #display(typeof(gamma3))
            #display(typeof(gamma2))
            #println()
            if prescreen
                #println("4b")
                #@btime upper_bound_thread($gamma1, $scr3, c=maximum(abs.($coef_ket)))
                #println(size(gamma1))
                #println(size(scr3))
                #println(typeof(gamma1))
                #println(typeof(scr3))
                upper_bound_thread(gamma1, scr3, c=maximum(abs.(coef_ket))) > thresh || continue
            end
            #println("c")

            mul!(scr4, gamma1', scr3)
            #a = @allocated begin
            #    mul!(scr4, gamma1', scr3)
            #end
            #a == 0 || error(" we allocated: ", a)
            #@btime mul!($scr4, $gamma1', $scr3)

            _collect_significant_thread!(sig, tmp_conf, scr4, coef_ket, c1.idx, c2.idx, newI, newJ, thresh, tmp2, state_sign)
        end
    end

    return
end
#=}}}=#

function _fill1!(gamma, idx, g, k)
    @inbounds for i in 1:size(g,1)
        for (II,I) in enumerate(idx)
            gamma[i,II] = g[i,I,k]
            #gamma[i,II] = deepcopy(g[i,I,k])
        end
    end
end

function _collect_significant_thread!(out, conf_ket, new_coeffs, coeff, c1idx, newI, thresh, tmp1)
    #={{{=#
    tmp1 .= conf_ket.config
    thresh_curr = thresh / maximum(abs.(coeff))
    for i::Int16 in newI
        if (new_coeffs[i] > thresh_curr) || (new_coeffs[i] < -thresh_curr)
            tmp1[c1idx] = i

            tmp = ClusterConfig(SVector(tmp1))
            if haskey(out, tmp)
                out[tmp] .+= new_coeffs[i].*coeff
            else
                out[tmp] = new_coeffs[i].*coeff
            end
        end
    end
end
#=}}}=#

function _collect_significant_thread!(out, conf_ket, new_coeffs, coeff, c1idx, c2idx, newI, newJ, thresh, tmp1, sign)
#={{{=#
    tmp1 .= conf_ket.config
    thresh_curr = thresh / maximum(abs.(coeff))

    ii::Int = 1
    nI = length(newI)
    nJ = length(newJ)
    @inbounds for j::Int16 in 1:nJ
        tmp1[c2idx] = newJ[j]
        for i::Int16 in 1:nI
            if abs(new_coeffs[i,j]) > thresh_curr
                tmp1[c1idx] = newI[i]
                tmp = ClusterConfig(SVector(tmp1))
                #@btime haskey($out, $tmp)
                if haskey(out, tmp)
                    ii += 1
                    if sign == 1
                        out[tmp] .+= new_coeffs[i,j] .* coeff
                    elseif sign == -1
                        out[tmp] .-= new_coeffs[i,j] .* coeff
                    else
                        error(" only 1 or -1")
                    end
                else
                    if sign == 1
                        out[tmp] = new_coeffs[i,j].*coeff
                    elseif sign == -1
                        out[tmp] = -new_coeffs[i,j].*coeff
                    else
                        error(" only 1 or -1")
                    end
                end
            end
        end
    end
end
#=}}}=#

"""
    upper_bound_thread(g1, g2; c::Float64=1.0)

Return upper bound on the size of matrix elements resulting from matrix multiply

    V[I,J] =  g1[i,I] * g2[i,J] * c

    max(|V|) <= sum_i max|g1[i,:]| * max|g2[i,:]| * |c|
"""
function upper_bound_thread(g1, g2; c::Float64=1.0)
#={{{=#
    bound::Float64 = 0.0
    n1 = size(g1,1)
    n2 = size(g2,1)
    n1 == n2 || throw(DimensionMismatch)

    absc = abs(c)
    @inbounds @simd for p in 1:n1
        @views pmax = maximum(abs, g1[p,:])
        @views qmax = maximum(abs, g2[p,:])
        bound += pmax * qmax
    end
    bound *= absc
    return bound
end
#=}}}=#

"""
    upper_bound_thread(v::AbstractArray{Float64,3}, g1, g2, g3, scr1, scr2, scr3; c::Float64=1.0)

Return upper bound on the size of tensor elements resulting from the following contraction

    V[I,J,K] = v[i,j,k] * g1[i,I] * g2[j,J] * g3[k,K]

    max(|V|) <= sum_ijk |v[ijk]| * |g1[i,:]|_8 * |g2[j,:]|_8 * |g3[k,:]|_8
"""
function upper_bound_thread(v::AbstractArray{Float64,3}, g1, g2, g3, scr1, scr2, scr3; c::Float64=1.0)
    #={{{=#
        bound = 0
        n1 = size(g1,1)
        n2 = size(g2,1)
        n3 = size(g3,1)

        resize!(scr1, n1)
        resize!(scr2, n2)
        resize!(scr3, n3)
        pmax = scr1
        qmax = scr2
        rmax = scr3

        fill!(pmax, 0.0)
        fill!(qmax, 0.0)
        fill!(rmax, 0.0)

        for p in 1:n1
            @views pmax[p] = maximum(abs(i) for i in g1[p,:])
            #@views pmax[p] = maximum(abs.(g1[p,:]))
        end
        for p in 1:n2
            @views qmax[p] = maximum(abs(i) for i in g2[p,:])
        end
        for p in 1:n3
            @views rmax[p] = maximum(abs(i) for i in g3[p,:])
        end

        tmp = 0.0
        @inbounds for r in 1:n3
            for q in 1:n2
                tmp = abs(c) * qmax[q] * rmax[r]
                @simd for p in 1:n1
                    bound += tmp * abs(v[p,q,r]) * pmax[p]
                end
            end
        end
    return bound
end
#=}}}=#

"""
    upper_bound_thread(v::Array{Float64,4}, g1, g2, g3, g4, scr1, scr2, scr3, scr4; c::Float64=1.0)

Return upper bound on the size of tensor elements resulting from the following contraction

    V[I,J,K,L] = v[i,j,k,l] * g1[i,I] * g2[j,J] * g3[k,K] * g4[l,L]

    max(|V|) <= sum_ijkl |v[ijkl]| * |g1[i,:]|_8 * |g2[j,:]|_8 * |g3[k,:]|_8 * |g4[l,:]|_8
"""
function upper_bound_thread(v::Array{Float64,4}, g1, g2, g3, g4, scr1, scr2, scr3, scr4; c::Float64=1.0)
    #={{{=#
        bound = 0
        n1 = size(g1,1)
        n2 = size(g2,1)
        n3 = size(g3,1)
        n4 = size(g4,1)

        resize!(scr1, n1)
        resize!(scr2, n2)
        resize!(scr3, n3)
        resize!(scr4, n4)
        pmax = scr1
        qmax = scr2
        rmax = scr3
        smax = scr4

        fill!(pmax, 0.0)
        fill!(qmax, 0.0)
        fill!(rmax, 0.0)
        fill!(smax, 0.0)

        for p in 1:n1
            @views pmax[p] = maximum(abs(i) for i in g1[p,:])
        end
        for p in 1:n2
            @views qmax[p] = maximum(abs(i) for i in g2[p,:])
        end
        for p in 1:n3
            @views rmax[p] = maximum(abs(i) for i in g3[p,:])
        end
        for p in 1:n4
            @views smax[p] = maximum(abs(i) for i in g4[p,:])
        end

        tmp = 0.0
        @inbounds for s in 1:n4
            for r in 1:n3
                for q in 1:n2
                    tmp = abs(c) * qmax[q] * rmax[r] * smax[s]
                    @simd for p in 1:n1
                        bound += tmp * abs(v[p,q,r,s]) * pmax[p]
                    end
                end
            end
        end
    return bound
end
#=}}}=#

"""
    upper_bound2_thread(v::Array{Float64,4}, g1, g2, g3,
        scr_f::Vector{Vector{Float64}}, scr_i::Vector{Vector{Int16}}, thresh;
        c::Float64=1.0)

    max(H_IJ(K)|_K <= sum_s (sum_pqr vpqrs max(g1[p,:]) * max(g2[q,:]) * |c| ) * |g3(r,K)|
"""
function upper_bound2_thread(v::Array{Float64,3}, g1, g2, g3,
        scr_f::Vector{Vector{Float64}}, scr_i::Vector{Vector{Int16}}, thresh;
        c::Float64=1.0)
    #={{{=#

        newI = scr_i[1]
        newJ = scr_i[2]
        newK = scr_i[3]
        resize!(newI,0)
        resize!(newJ,0)
        resize!(newK,0)

        n1 = size(v,1)
        n2 = size(v,2)
        n3 = size(v,3)

        n1 == size(g1,1) || throw(DimensionMismatch)
        n2 == size(g2,1) || throw(DimensionMismatch)
        n3 == size(g3,1) || throw(DimensionMismatch)

        pmax = scr_f[5]
        qmax = scr_f[6]
        rmax = scr_f[7]

        resize!(pmax,n1)
        resize!(qmax,n2)
        resize!(rmax,n3)

        fill!(pmax, 0.0)
        fill!(qmax, 0.0)
        fill!(rmax, 0.0)

        for p in 1:n1
            @views pmax[p] = maximum(abs, g1[p,:])
        end
        for p in 1:n2
            @views qmax[p] = maximum(abs, g2[p,:])
        end
        for p in 1:n3
            @views rmax[p] = maximum(abs, g3[p,:])
        end

        tmp = 0.0

        mI = scr_f[9]
        resize!(mI,size(g1,2))
        fill!(mI, 0.0)
        @inbounds for r in 1:n3
            for q in 1:n2
                tmp = qmax[q] * rmax[r] * abs(c)
                for p in 1:n1
                    @views @. mI += abs(v[p,q,r]) * abs.(g1[p,:]) * tmp
                end
            end
        end

        mJ = scr_f[10]
        resize!(mJ,size(g2,2))
        fill!(mJ, 0.0)
        @inbounds for r in 1:n3
            for p in 1:n1
                tmp = pmax[p] * rmax[r] * abs(c)
                for q in 1:n2
                    @views @. mJ += abs(v[p,q,r]) * abs.(g2[q,:])  * tmp
                end
            end
        end

        mK = scr_f[11]
        resize!(mK,size(g3,2))
        fill!(mK, 0.0)
        @inbounds for q in 1:n2
            for p in 1:n1
                tmp = pmax[p] * qmax[q] * abs(c)
                for r in 1:n3
                    @views @. mK += abs(v[p,q,r]) * abs.(g3[r,:]) * tmp
                end
            end
        end

        for I in 1:size(g1,2)
            if abs(mI[I]) > thresh
                push!(newI,I)
            end
        end

        for J in 1:size(g2,2)
            if abs(mJ[J]) > thresh
                push!(newJ,J)
            end
        end

        for K in 1:size(g3,2)
            if abs(mK[K]) > thresh
                push!(newK,K)
            end
        end

    return newI, newJ, newK
end
#=}}}=#

"""
    upper_bound2_thread(v::Array{Float64,4}, g1, g2, g3, g4,
        scr_f::Vector{Vector{Float64}}, scr_i::Vector{Vector{Int16}}, thresh;
        c::Float64=1.0)

    max(H_IJK(L)|_L <= sum_s (sum_pqr vpqrs max(g1[p,:]) * max(g2[q,:]) * max(g3[r,:]) * |c| ) * |g4(s,L)|
"""
function upper_bound2_thread(v::Array{Float64,4}, g1, g2, g3, g4,
        scr_f::Vector{Vector{Float64}}, scr_i::Vector{Vector{Int16}}, thresh;
        c::Float64=1.0)
    #={{{=#

        newI = scr_i[1]
        newJ = scr_i[2]
        newK = scr_i[3]
        newL = scr_i[4]
        resize!(newI,0)
        resize!(newJ,0)
        resize!(newK,0)
        resize!(newL,0)

        n1 = size(v,1)
        n2 = size(v,2)
        n3 = size(v,3)
        n4 = size(v,4)

        n1 == size(g1,1) || throw(DimensionMismatch)
        n2 == size(g2,1) || throw(DimensionMismatch)
        n3 == size(g3,1) || throw(DimensionMismatch)
        n4 == size(g4,1) || throw(DimensionMismatch)

        pmax = scr_f[5]
        qmax = scr_f[6]
        rmax = scr_f[7]
        smax = scr_f[8]

        resize!(pmax,n1)
        resize!(qmax,n2)
        resize!(rmax,n3)
        resize!(smax,n4)

        fill!(pmax, 0.0)
        fill!(qmax, 0.0)
        fill!(rmax, 0.0)
        fill!(smax, 0.0)

        for p in 1:n1
            @views pmax[p] = maximum(abs, g1[p,:])
        end
        for p in 1:n2
            @views qmax[p] = maximum(abs, g2[p,:])
        end
        for p in 1:n3
            @views rmax[p] = maximum(abs, g3[p,:])
        end
        for p in 1:n4
            @views smax[p] = maximum(abs, g4[p,:])
        end

        tmp = 0.0

        mI = scr_f[9]
        resize!(mI,size(g1,2))
        fill!(mI, 0.0)
        @inbounds for s in 1:n4
            for r in 1:n3
                for q in 1:n2
                    tmp = qmax[q] * rmax[r] * smax[s] * abs(c)
                    for p in 1:n1
                        @views @. mI += abs(v[p,q,r,s]) * abs.(g1[p,:]) * tmp
                    end
                end
            end
        end

        mJ = scr_f[10]
        resize!(mJ,size(g2,2))
        fill!(mJ, 0.0)
        @inbounds for s in 1:n4
            for r in 1:n3
                for p in 1:n1
                    tmp = pmax[p] * rmax[r] * smax[s] * abs(c)
                    for q in 1:n2
                        @views @. mJ += abs(v[p,q,r,s]) * abs.(g2[q,:])  * tmp
                    end
                end
            end
        end

        mK = scr_f[11]
        resize!(mK,size(g3,2))
        fill!(mK, 0.0)
        @inbounds for s in 1:n4
            for q in 1:n2
                for p in 1:n1
                    tmp = pmax[p] * qmax[q] * smax[s] * abs(c)
                    for r in 1:n3
                        @views @. mK += abs(v[p,q,r,s]) * abs.(g3[r,:]) * tmp
                    end
                end
            end
        end

        mL = scr_f[12]
        resize!(mL,size(g4,2))
        fill!(mL, 0.0)
        @inbounds for r in 1:n3
            for q in 1:n2
                for p in 1:n1
                    tmp =  pmax[p] * qmax[q] * rmax[r] * abs(c)
                    for s in 1:n4
                        @views @. mL += abs(v[p,q,r,s]) * abs.(g4[s,:]) * tmp
                    end
                end
            end
        end

        for I in 1:size(g1,2)
            if abs(mI[I]) > thresh
                push!(newI,I)
            end
        end

        for J in 1:size(g2,2)
            if abs(mJ[J]) > thresh
                push!(newJ,J)
            end
        end

        for K in 1:size(g3,2)
            if abs(mK[K]) > thresh
                push!(newK,K)
            end
        end

        for L in 1:size(g4,2)
            if abs(mL[L]) > thresh
                push!(newL,L)
            end
        end

    return newI, newJ, newK, newL
end
#=}}}=#

max_stars_repo_name: nmayhall-vt/FermiCG

max_issues_repo_name: nmayhall-vt/FermiCG

max_forks_repo_name: nmayhall-vt/FermiCG


text:
export DaxConstrained

mutable struct DaxConstrained{matT,T,Tsparse,U} <: AbstractLinearSolver
  S::matT
  u::Vector{T}
  λ::Float64
  B::Tsparse
  Bnorm²::Vector{Float64}
  denom::Vector{U}
  rowindex::Vector{Int64}
  zk::Vector{T}
  bk::Vector{T}
  bc::Vector{T}
  xl::Vector{T}
  yl::Vector{T}
  yc::Vector{T}
  δc::Vector{T}
  εw::Vector{T}
  τl::T
  αl::T
  weights::Vector{U}
  iterations::Int64
  iterationsInner::Int64
end

"""This function solves a constrained linear least squares problem using an algorithm proposed in [1].
Returns an approximate solution to Sᵀx = u s.t. Bx>=0 (each component >=0).

[1] Dax, A. On Row Relaxation Methods for Large Constrained Least Squares Problems. SIAM J. Sci. Comput. 14, 570–584 (1993).

### Input Arguments
* `S::AbstractMatrix{T}`: Problem Matrix S.
* `u::Vector{T}`: Righthandside of the linear equation.
* `B::AbstractMatrix{T}`: Transpose of Basistransformation if solving in dual space.

### Keyword/Optional Arguments

* `iterations::Int`: Number of Iterations of outer dax scheme.
* `iterationsInner::Int`: Number of Iterations of inner dax scheme.
* `λ::Float64`: The regularization parameter ɛ>0 influences the speed of convergence but not the solution.
* `weights::Bool`: Use weights in vector to weight equations. The larger the weight the more one 'trusts' a sqecific equation.
* `B`: Basistransformation if solving in dual space.
"""
function DaxConstrained(S, b=nothing; λ::Real=0.0
              , weights::Vector{R}=ones(Float64,size(S,1))
              , sparseTrafo=nothing
              , enforceReal::Bool=false
              , enforcePositive::Bool=false
              , iterations::Int=3
              , iterationsInner::Int=2
              , kargs...) where R <: Real

  T = typeof(real(S[1]))
  M,N = size(S)

  # setup denom and rowindex
  denom, rowindex = initkaczmarzconstraineddax(S,λ,weights)

  # set basis transformation
  sparseTrafo==nothing ? B=Matrix{T}(I, size(S,2), size(S,2)) : B=sparseTrafo
  Bnorm² = [rownorm²(B,i) for i=1:size(B,2)]

  if b != nothing
    u = b
  else
    u = zeros(eltype(S),M)
  end

  zk = zeros(eltype(S),N)
  bk = zeros(eltype(S),M)
  bc = zeros(T,size(B,2))
  xl = zeros(eltype(S),N)
  yl = zeros(eltype(S),M)
  yc = zeros(eltype(S),N)
  δc = zeros(eltype(S),N)
  εw = zeros(eltype(S),length(rowindex))
  τl = zero(eltype(S))
  αl = zero(eltype(S))

  return DaxConstrained(S,u,Float64(λ),B,Bnorm²,denom,rowindex,zk,bk,bc,xl,yl,yc,δc,εw,τl,αl
                  ,T.(weights),iterations,iterationsInner)
end

function init!(solver::DaxConstrained
              ; S::matT=solver.S
              , λ::Real=solver.λ
              , u::Vector{T}=eltype(S)[]
              , zk::Vector{T}=eltype(S)[]
              , weights::Vector{Float64}=solver.weights) where {matT,T}

  if S != solver.S
    denom, rowindex = initkaczmarzconstraineddax(S,λ,weights)
  end
  solver.λ = Float64(λ)

  solver.u[:] .= u
  solver.weights=weights

  # start vector
  if isempty(zk)
    solver.zk[:] .= zeros(T,size(S,2))
  else
    solver.zk[:] .= x
  end

  solver.bk[:] .= zero(T)
  solver.bc[:] .= zero(T)
  solver.xl[:] .= zero(T)
  solver.yl[:] .= zero(T)
  solver.yc[:] .= zero(T)
  solver.δc[:] .= zero(T)
  solver.αl = zero(T)        #temporary scalar
  solver.τl = zero(T)        #temporary scalar

  for i=1:length(solver.rowindex)
    j = solver.rowindex[i]
    solver.ɛw[i] = sqrt(solver.λ)/weights[j]
  end
end

function solve(solver::DaxConstrained, u::Vector{T}; λ::Real=solver.λ
                , S::matT=solver.S, startVector::Vector{T}=eltype(S)[]
                , weights::Vector=solver.weights
                , solverInfo=nothing, kargs...) where {T,matT}

  # initialize solver parameters
  init!(solver; S=S, λ=λ, u=u, zk=startVector, weights=weights)

  # log solver information
  solverInfo != nothing && storeInfo(solverInfo,solver.zk,norm(solver.bk))

  # perform CGNR iterations
  for (iteration, item) = enumerate(solver)
    solverInfo != nothing && storeInfo(solverInfo,solver.zk,norm(solver.bk))
  end

  return solver.zk
end

function iterate(solver::DaxConstrained, iteration::Int=0)
  if done(solver,iteration) return nothing end

  # bk = u-S'*zk
  copyto!(solver.bk,solver.u)
  gemv!('N',-1.0,solver.S,solver.zk,1.0,solver.bk)

  # solve min ɛ|x|²+|W*A*x-W*bk|² with weightingmatrix W=diag(wᵢ), i=1,...,M.
  for l=1:solver.iterationsInner
    for i=1:length(solver.rowindex) # perform kaczmarz for all rows, which receive an update.
      j = solver.rowindex[i]
      solver.τl = dot_with_matrix_row(solver.S,solver.xl,j)
      solver.αl = solver.denom[i]*(solver.bk[j]-solver.τl-solver.ɛw[i]*solver.yl[j])
      kaczmarz_update!(solver.S,solver.xl,j,solver.αl)
      solver.yl[j] += solver.αl*solver.ɛw[i]
    end

    #Lent-Censor scheme for ensuring B(xl+zk) >= 0
    # copyto!(solver.δc,solver.xl)
    # BLAS.axpy!(1.0,solver.zk,solver.δc)
    # lmul!(solver.B, solver.δc)
    # lentcensormin!(solver.δc,solver.yc)
    #
    # BLAS.axpy!(1.0,solver.δc,solver.yc) # yc += δc
    #
    # δc = gemv('T',solver.B,solver.δc)
    # BLAS.axpy!(1.0,solver.solver.δc,solver.xl) # xl += Bᵀ*δc

    #Lent-Censor scheme for solving Bx >= 0
    # bc = xl + zk
    copyto!(solver.bc,solver.zk)
    BLAS.axpy!(1.0,solver.xl,solver.bc)

    for i=1:size(solver.B,2)
      δ = dot_with_matrix_row(solver.B,solver.bc,i)/solver.Bnorm²[i]
      δ = δ<solver.yc[i] ? -δ : -solver.yc[i]
      solver.yc[i] += δ
      kaczmarz_update!(solver.B,solver.xl,i,δ)  # update xl
      kaczmarz_update!(solver.B,solver.bc,i,δ)  # update bc
    end
  end
  BLAS.axpy!(1.0,solver.xl,solver.zk)  # zk += xl

  # reset xl and yl for next Kaczmarz run
  rmul!(solver.xl,0.0)
  rmul!(solver.yl,0.0)
  rmul!(solver.yc,0.0)

  return solver.bk, iteration+1

end

@inline done(solver::DaxConstrained,iteration::Int) = iteration>=solver.iterations

function lentcensormin!(x::Vector{T},y::Vector{T}) where {T<:Real}
  α = -one(T)
  @simd for i=1:length(x)
    @inbounds x[i]<=y[i] ? x[i] = α*x[i] : x[i] = α*y[i]
  end
end

"""This funtion saves the denominators to compute αl in denom and the rowindices,
  which lead to an update of cl in rowindex."""
function initkaczmarzconstraineddax(S::AbstractMatrix,ɛ::Number,weights::Vector)
  length(weights)==size(S,1) ? nothing : error("number of weights must equal number of equations")
  denom = Float64[]
  rowindex = Int64[]

  for i=1:size(S,1)
    s² = rownorm²(S,i)*weights[i]^2
    if s²>0
      push!(denom,weights[i]^2/(s²+ɛ))
      push!(rowindex,i)
    end
  end
  return denom, rowindex
end

"""
This funtion saves the denominators to compute αl in denom and the rowindices,
which lead to an update of cl in rowindex.
"""
function initkaczmarzconstraineddaxfft(S::AbstractMatrix,ɛ::Number,weights::Vector)
  length(weights)==size(S,1) ? nothing : error("number of weights must equal number of equations")
  denom = Float64[]
  rowindex = Int64[]

  for i=1:2:size(S,1)
    s²a = (rownorm²(S,i)+rownorm²(S,i+1))*weights[i]^2
    s²b = (rownorm²(S,i)+rownorm²(S,i+1))*weights[i+1]^2
    if s²a>0 && s²b>0
      push!(denom,weights[i]^2/(s²a+ɛ))
      push!(denom,weights[i+1]^2/(s²b+ɛ))
      push!(rowindex,i)
      push!(rowindex,i+1)
    end
  end
  return denom, rowindex
end

ext: jl


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

mutable struct DiskInstanceView <: SwaggerModel
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    encryptionSettings::Any # spec type: Union{ Nothing, Vector{DiskEncryptionSettings} } # spec name: encryptionSettings
    statuses::Any # spec type: Union{ Nothing, Vector{InstanceViewStatus} } # spec name: statuses

    function DiskInstanceView(;name=nothing, encryptionSettings=nothing, statuses=nothing)
        o = new()
        validate_property(DiskInstanceView, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(DiskInstanceView, Symbol("encryptionSettings"), encryptionSettings)
        setfield!(o, Symbol("encryptionSettings"), encryptionSettings)
        validate_property(DiskInstanceView, Symbol("statuses"), statuses)
        setfield!(o, Symbol("statuses"), statuses)
        o
    end
end # type DiskInstanceView

const _property_map_DiskInstanceView = Dict{Symbol,Symbol}(Symbol("name")=>Symbol("name"), Symbol("encryptionSettings")=>Symbol("encryptionSettings"), Symbol("statuses")=>Symbol("statuses"))
const _property_types_DiskInstanceView = Dict{Symbol,String}(Symbol("name")=>"String", Symbol("encryptionSettings")=>"Vector{DiskEncryptionSettings}", Symbol("statuses")=>"Vector{InstanceViewStatus}")
Base.propertynames(::Type{ DiskInstanceView }) = collect(keys(_property_map_DiskInstanceView))
Swagger.property_type(::Type{ DiskInstanceView }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_DiskInstanceView[name]))}
Swagger.field_name(::Type{ DiskInstanceView }, property_name::Symbol) =  _property_map_DiskInstanceView[property_name]

function check_required(o::DiskInstanceView)
    true
end

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


text:
## Abstract  types for  orthogonal  polynomials

export Basis

## Has An(P), Bn(P), Cn(P)
abstract type AbstractOrthogonalPolynomial{T,X} <: AbstractSpecialPolynomial{T,X} end
abstract type AbstractContinuousOrthogonalPolynomial{T,X} <:
              AbstractOrthogonalPolynomial{T,X} end
abstract type AbstractDiscreteOrthogonalPolynomial{T,X} <: AbstractOrthogonalPolynomial{T,X} end

"""
    AbstractOrthogonalPolynomial{T,X}

Type to represent systems of orthogonal polynomials. These polynomials have  several properties, including an accompanying inner product satsifying  `⟨yᵢ, yⱼ⟩ = cᵢδᵢⱼ`.

In addition to methods inherited from the underlying `AbstractPolynomial`  type, orthogonal polynomial  types may have methods   `weight_function`, `generating_function`, `leading_term`, `norm2`, `jacobi_matrix`, and `gauss_nodes_weights`,  though none are  exported.

Subtypes of `AbstractCOP <: AbstractOrthogonalPolynomial` utilize the fact that the basis  polynomials  satisfy

`(ax² + bx + c)yᵢ''(x) + (dx+e)yᵢ'(x) + λᵢyᵢ(x) = 0` (or a discrete analogue)

where the structural relations are functions of `a,b,c,d,e`. These allow default definitions for polynomial evaluation,   addition, multiplication, differentiation, integration, and  conversion to and from  the `Polynomial` type (the `FallingFactorial` type in the discrete  c case),

A key structural relation is the three-term recursion,  `yᵢ₊₁ =  (Aᵢx +  Bᵢ)yᵢ -  Cᵢyᵢ₋₁`. For systems  specfied by  a  weight function, the  values of `Aᵢ`, `Bᵢ`, and `Cᵢ` can  be  generated, yielding formulas for polynomial evaluation, addition, and conversion to the `Polynomial`  type throughe evaluation.
"""
AbstractOrthogonalPolynomial

##
## --------------------------------------------------
##
function Polynomials.fromroots(
    P::Type{<:AbstractOrthogonalPolynomial},
    roots::AbstractVector;
    var::Polynomials.SymbolLike=:x,
)
    x = variable(P, var)
    p = prod(x - r for r in roots)
end

Base.convert(P::Type{<:AbstractOrthogonalPolynomial}, c::Number) = c * one(P)
Base.one(P::Type{<:AbstractOrthogonalPolynomial}, var::Polynomials.SymbolLike=:x) =
    basis(P, 0, var) / k0(P)

Polynomials.variable(
    P::Type{<:AbstractOrthogonalPolynomial},
    var::Polynomials.SymbolLike=:x,
) = (basis(P, 1, var) / k0(P) - Bn(P, 0)) / An(P, 0)

## Evaluation

# from type, cs, x
function clenshaw_eval(P::Type{<:AbstractOrthogonalPolynomial{T}}, cs, x::S) where {T,S}
    N = length(cs)
    p₀ = k0(P)
    R = promote_type(promote_type(T, S), typeof(An(P, 0)))
    N == 0 && return zero(R)
    N == 1 && return (cs[1] * p₀) * one(R)

    Δ0::R = cs[end - 1]
    Δ1::R = cs[end]
    @inbounds for i in (N - 1):-1:2
        Δ0, Δ1 =
            cs[i - 1] - Δ1 * Cn(P, i - 1), Δ0 + Δ1 * muladd(x, An(P, i - 1), Bn(P, i - 1))
    end
    p₁ = muladd(x, An(P, 0), Bn(P, 0)) * p₀
    return Δ0 * p₀ + Δ1 * p₁
end

## Connection/Linearization

#  Structs for connection, linearization (connection.jl)
struct Connection{P,Q}
    n::Int
    k::Int
end

struct Linearization{P,V}
    l::Int
    n::Int
    m::Int
end

## Collect facts  about  the orthogonal polynomial system

"""
    weight_function(p)
    weight_function(::Type{P})

For an orthogonal polynomial type, a function `w` with `∫ B_n(t) B_m(t) w(t) dt = 0` when `n` and `m` are not equal.

"""
weight_function(::Type{P}) where {P<:AbstractOrthogonalPolynomial} =
    throw(MethodError("Not implemented"))
weight_function(::P) where {P<:AbstractOrthogonalPolynomial} = weight_function(P)

"""
    generating_function(p)
    generating_function(::Type{P})

The generating function is a function defined by: `(t,x) -> sum(t^n Pn(x) for n in 0:oo)`.
"""
generating_function(::Type{P}) where {P<:AbstractOrthogonalPolynomial} =
    throw(MethodError("Not implemented"))
generating_function(::P) where {P<:AbstractOrthogonalPolynomial} = generating_function(P)

"""
    leading_term(::Type{P},n)

Return leading term of `basis(P,n)` in the  standard basis. By default this is generated through the three-point recursion.
"""
function leading_term(::Type{P}, n::Int) where {P<:AbstractOrthogonalPolynomial}
    n < 0 && throw(ArgumentError("n must be a non-negative integer"))
    n == 0 && return one(eltype(P))
    prod(An(P, i) for i in (n - 1):-1:0)
end

# is P a monic polynomial system?
ismonic(::Type{P}) where {P<:AbstractOrthogonalPolynomial} = false
isorthonormal(::Type{P}) where {P<:AbstractOrthogonalPolynomial} = false

# cf. https://en.wikipedia.org/wiki/Orthogonal_polynomials#Recurrence_relation
# Orthogonal polynomials have a three-point recursion formula
# parameterized here through:
# P_{n+1} = (An*x + Bn) * P_n + Cn P_{n-1}
# also commonly  written as
# An⋅x⋅P_n = -P_{n+1} + Bn⋅P_n + C_n P_{n-1}
"""
    An(::Type{P},n)

Orthogonal polynomials defined by a weight function satisfy a three point recursion formula of the form:

`P_{n+1} = (A_n x + B_n) P_{n} - C_n P_{n-1}`

If the polynomials are monic, this is usually parameterized as:

`π_{n+1} = (x - α̃_n) π_n - β̃_n π_{n-1}`

These functions are used through recursion when evaluating the polynomials, converting to `Polynomial` format, for constructing the Vandermonde matrix, for construction the Jacobi matrix, and elsewhere.
"""
An(::Type{P}, n) where {P<:AbstractOrthogonalPolynomial} =
    throw(ArgumentError("No default method"))

"""
    Bn(::Type{P},n)
    Bn(p::P, n)

cf. [`An()`](@ref)
"""
Bn(::Type{P}, n) where {P<:AbstractOrthogonalPolynomial} =
    throw(ArgumentError("No default method"))

"""
    Cn(::Type{P},n)
    Cn(p::P, n)

cf. [`An()`](@ref)
"""
Cn(::Type{P}, n) where {P<:AbstractOrthogonalPolynomial} =
    throw(ArgumentError("No default method"))

## For monic polynomials, we have
## π_{n+1} = (x - α̃(n)) π_{n} - β̃(n)π_{n-1}
"""
    π̃αn(::Type{P}, n)

cf. [`An()`](@ref)
"""
π̃αn(P::Type{<:AbstractOrthogonalPolynomial}, n) = -Bn(P, n) / An(P, n)

"""
    β̃n(::Type{P}, n)

cf. [`An()`](@ref)
"""
function π̃βn(P::Type{<:AbstractOrthogonalPolynomial}, n)
    iszero(n) && return innerproduct(P, one, one)
    Cn(P, n) / An(P, n) / An(P, n - 1)
end

"""
    monic(p::AbstractOrthogonalPolynomial)

Return `p` as a monic polynomial *when* represented in the standard basis. Returns the zero polynomial if the degree of `p` is `-1`.
"""
function monic(p::P) where {P<:AbstractOrthogonalPolynomial}
    n = degree(p)
    n == -1 && return ⟒(P)(0 / one(eltype(p)))
    p / (p[end] * leading_term(P, n))
end

##
## Work with compact  basis
##
"""
     Basis(P,n)
     Basis{P}(n)

The  command `basis(P,n, [var])` realizes the polynomial. `Basis(P,n)` does  not.  This  can  be  useful for  evaluation.
"""
struct Basis{Π}
    n::Int
end
Basis(P, n) = Basis{P}(n)
Base.show(io::IO, mimetype::MIME"text/plain", b::Basis{P}) where {P} =
    print(io, "$(P)($(b.n))")

function innerproduct(::P, f::Basis{P}, g::Basis{P}) where {P}
    n, m = f.n, g.n
    if n == m
        return norm2(P, n)
    else
        return zero(P)
    end
end

##
## Vandermonde matrix can be generated through the 3-point recursion formula
##
function Polynomials.vander(
    p::Type{P},
    x::AbstractVector{T},
    n::Integer,
) where {P<:AbstractOrthogonalPolynomial,T<:Number}
    A = Matrix{T}(undef, length(x), n + 1)

    # for i in 0:n
    #     A[:, i+1] = basis(P, i).(x)
    # end
    # return A

    A[:, 1] .= basis(P, 0)(one(T)) #P0(P, one(T))

    if n > 0
        A[:, 2] .= basis(P, 1).(x) #P1(P, x)
        @inbounds for i in 1:(n - 1)
            # `P_{n+1} = (A_n x + B_n) P_{n} - C_n P_{n-1}`
            n′ = i + 1
            A[:, n′ + 1] =
                (An(p, n′) * x .+ Bn(p, n′)) .* A[:, n′] .- (Cn(p, n′) * A[:, n′ - 1])
        end
    end

    return A
end

# Jacobi matrix
# roots(basis(P,n)) = eigvals(jacobi_matrix(P,n)), but this is more stable
# https://en.wikipedia.org/wiki/Gaussian_quadrature#The_Golub-Welsch_algorithm
"""
    jacobi_matrix(::Type{P}, n)

The Jacobi Matrix is a symmetric tri-diagonal matrix. The diagonal entries are the `alphaᵢ` values, the off diagonal entries,
the square root of the  `betaᵢ` values. This matrix has the properties that

* the eigenvalues are the roots of the corresponding basis vector. As these roots are important in quadrature, and finding eigenvalues of a symmetric tri-diagonal matrix yields less error than finding the eigenvalues of the companion matrix, this can be used for higher degree basis polynomials.
* the normalized eigenvectors have initial term proportional to the weights in a quadrature formula

"""
function jacobi_matrix(::Type{P}, n) where {P<:AbstractOrthogonalPolynomial}
    LinearAlgebra.SymTridiagonal(
        [π̃αn(P, i) for i in 0:(n - 1)],
        [sqrt(π̃βn(P, i)) for i in 1:(n - 1)],
    )
end
jacobi_matrix(p::P, n) where {P<:AbstractOrthogonalPolynomial} = jacobi_matrix(P, n)

##  Compute weights and nodes for quadrature
"""
    gauss_nodes_weights(::Type{P}, n)

Returns a tuple of nodes and weights for Gauss quadrature for the given orthogonal type.

When available, the values are computed through  the `FastGaussQuadratures` package.

For some types, a method from  A. Glaser, X. Liu, and V. Rokhlin. "A fast algorithm for the calculation of the roots of special functions." SIAM J. Sci. Comput., 29 (2007), 1420-1438. is used.

For others the Jacobi matrix, J_n, for which the Golub-Welsch] algorithm The nodes  are computed from the eigenvalues of J_n, the weights a scaling of the first component of the normalized eigen vectors (β_0 * [v[1] for v in vs])

"""
function gauss_nodes_weights(p::Type{P}, n) where {P<:AbstractOrthogonalPolynomial}
    J = jacobi_matrix(P, n)
    eig = eigen(J, extrema(P)...)
    # Is this necessary?
    nm  = 1 #diag(eig.vectors * eig.vectors')
    wts = π̃βn(P, 0) * (eig.vectors[1, :] ./ nm) .^ 2
    eig.values, wts
end

gauss_nodes_weights(B::Basis{P}) where {P} = gauss_nodes_weights(B.P, B.n)

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

"""
      innerproduct(::Type{P}, f, g)

Compute  `<f,g> = ∫ f⋅g⋅w dx` where  `w` is the weight function of the  type  `P`  and the integral is  taken  over  the domain of the type `P`.
"""
function innerproduct(
    P::Type{<:Union{AbstractOrthogonalPolynomial}},
    f,
    g;
    atol=sqrt(eps(float(eltype(P)))),
)
    dom = domain(P)
    a, b = first(dom), last(dom)
    if first(bounds_types(dom)) == Open
        a += eps(float(one(a)))
    end
    if last(bounds_types(dom)) == Open
        b -= eps(float(one(b)))
    end
    fn = x -> f(x) * g(x) * weight_function(P)(x)

    return quadgk(fn, a, b; atol=atol)[1]
end

## Compute <p_i, p_i> = \| p \|^2; allows export, but work is in norm2
Base.abs2(::Type{P}, n) where {P<:AbstractOrthogonalPolynomial} = norm2(P, n)

## Compute <p_i, p_i> = \| p \|^2
## Slow default; generally should  be directly expressed for each type
function norm2(::Type{P}, n) where {P<:AbstractOrthogonalPolynomial}
    p = basis(P, n)
    innerproduct(P, p, p)
end

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

"""
    fit([Val(S)], P::Type{<:AbstractOrthogonalPolynomial}, f, n::Int; var=:x)

Find an approximating polynomial of degree `n` or less for a function `f`, that is returns `p(x) = ∑ᵢⁿ cᵢ Pᵢ(x)` for some coefficients `cᵢ`.

Defaults to an interpolating polynomial. To specify others, use one of `Val(:interpolating)`, `Val(:lsq)` (least squares), or `Val(:series)` (trunated series expansion) as the first argument. See [`SpecialPolynomials.cks`](@ref) for some more detail.

"""
Polynomials.fit(P::Type{<:AbstractOrthogonalPolynomial}, f, n::Int; var=:x) =
    fit(Val(:interpolating), P, f, n, var=var)

"""
    fit(val::Val{:interpolating}, P::Type{<:AbstractOrthogonalPolynomial}, f, deg::Int; var=:x)

Fit `f` with an interpolating polynomial of degree `n` or less using nodes
`x0,x1, ..., xn`, the  zeros of `P_{n+1} = basis(P, n+1)`. and `p(xᵢ)=f(xᵢ)`.
"""
Polynomials.fit(
    val::Val{:interpolating},
    P::Type{<:AbstractOrthogonalPolynomial},
    f,
    n::Int;
    var=:x,
) = P(cks(val, P, f, n), var)

"""
    fit(Val(:lsq), P::Type{<:AbstractOrthogonalPolynomial}, f, n::Int)

Fit `f` with `p(x)=∑ d_i P_i(x)` where `p` had degree `n` or less using least squares.
"""
Polynomials.fit(
    val::Val{:lsq},
    P::Type{<:AbstractOrthogonalPolynomial},
    f,
    n::Int;
    var=:x,
) = P(cks(val, P, f, n), var)

"""
    fit(Val(:series), P::Type{<:AbstractOrthogonalPolynomial}, f, n::Int)

If `f(x)` is written as an infinite sum `∑ c_kP_k(x)`, this returns a truncated sum `∑_0^n c̃_k P_k(x)` where `n` is chosen algorithmically and `c̃_k` is chosen using an efficient manner, not necessarily through the orthogonality condition, `<f,p_k>/<p_k,p_k>`.

"""
Polynomials.fit(
    val::Val{:series},
    P::Type{<:AbstractOrthogonalPolynomial},
    f;
    var=:x,
    kwargs...,
) = P(cks(val, P, f; kwargs...), var)

"""
    cks(::Val{:interpolating}, ::Type{P}, f, n::Int)

Fit `f` with the interpolating polynomial using roots of `P_{n+1}`.

Let `xs`, `ws` be the gauss nodes and weights of `P` (`xs` are the zeros of `P_{n+1}`).

Then if we interpolate `f` at the `xs` to get `p`, then
`p(xᵢ) = f(xᵢ)` and

`∑ p_n(xᵢ) p_m(xᵢ) wᵢ = K_m δ_{nm}` and

`p(x) = ∑ᵢ cᵢ Pᵢ(x)`.

Using this:
`∑ᵢ f(xᵢ)P_k(xᵢ) wᵢ =`
`∑ᵢ (∑_j c_j P_j(xᵢ)) P_k(xᵢ) wᵢ =`
`∑_j c_j (K_k δ_{ik}) = c_k K_k`, So
`c_k = (1/K_k) ∑ᵢ f(xᵢ)P_k(xᵢ) wᵢ`
"""
function cks(
    ::Val{:interpolating},
    ::Type{P},
    f,
    n::Int,
) where {P<:AbstractOrthogonalPolynomial}
    xs, ws = gauss_nodes_weights(P, n)
    return [
        sum(f(xⱼ) * basis(P, k)(xⱼ) * wⱼ for (xⱼ, wⱼ) in zip(xs, ws)) / norm2(P, k) for
        k in 0:n
    ]
end

"""
    cks(::Val{:lsq}, ::Type{P}, f, n::Int)

Fit `f` with a polynomial `∑ᵢⁿ cᵢ Pᵢ` chosen so `<f-p,f-p>_w` is as small as possible.
 Using the normal equations, the coefficients are found to be `c_k = <f,P_k>_w / <P_k,P_k>_w.`
For some types an approximation to the inner product, `<f,P_k>_w` may be used.

ref: [http://www.math.niu.edu/~dattab/MATH435.2013/APPROXIMATION](http://www.math.niu.edu/~dattab/MATH435.2013/APPROXIMATION)
"""
function cks(::Val{:lsq}, ::Type{P}, f, n::Int) where {P<:AbstractOrthogonalPolynomial}
    ## return ck =  <f,P_k>/<P_k,P_k>, k =  0...n
    [innerproduct(P, f, basis(P, k)) / norm2(P, k) for k in 0:n]
end

"""
    cks(::Val{:series}, ::Type{P}, f, n::Int)

If `f(x)` is written as an infinite sum `∑ c_kP_k(x)`, then this
tries to identify an `n` for which the series expansion is a good approximation and returns the coefficients.

"""
function cks(::Val{:series}, ::Type{P}, f, n::Int) where {P<:AbstractOrthogonalPolynomial}
    throw(ArgumentError("No default method"))
end


text:
using BinDeps

@BinDeps.setup

nlopt = library_dependency("libnlopt")
bayesopt = library_dependency("libbayesopt")
version = v"0.8.2"
commit = "6f2b74381f97"
srcdir = joinpath(BinDeps.depsdir(bayesopt), "src", "rmcantin-bayesopt-$(commit)")
prefix = joinpath(BinDeps.depsdir(bayesopt), "usr")

provides(Sources,
         URI("https://bitbucket.org/rmcantin/bayesopt/downloads/bayesopt-v$version.zip"),
         [nlopt, bayesopt],
         unpacked_dir="rmcantin-bayesopt-$(commit)")

provides(SimpleBuild,
         (@build_steps begin
              GetSources(bayesopt)
              @build_steps begin
                  ChangeDirectory(srcdir)
                  `cmake -DCMAKE_INSTALL_PREFIX="$prefix" -DCMAKE_INSTALL_RPATH="$prefix/lib" -DCMAKE_INSTALL_RPATH_USE_LINK_PATH=ON -DCMAKE_MACOSX_RPATH=ON -DBAYESOPT_BUILD_SHARED=ON -DNLOPT_BUILD_SHARED=ON`
                  `make`
                  `make install`
              end
          end), [nlopt, bayesopt], os=:Unix)

@BinDeps.install Dict(:libnlopt => :libnlopt, :libbayesopt => :libbayesopt)


text:
## parsing
td = :(@traitdef Cr20{X} begin
    length(X)
end)
a,b,c = Traits.parsebody(td.args[end])
# a is not hard to test because of the random gensym
@test a.head==:call
@test a.args[1]==:(Traits.FDict)
@test a.args[2].head==:(=>)
@test a.args[2].args[1] == :length
@test a.args[2].args[2].args[1].args[2] == :(::Any)
@test a.args[2].args[2].args[1].args[3] == :(::X)
@test a.args[2].args[2].args[2] == :(nothing)
@test b==:(Bool[])
@test c.args[1]==:(assoctyps = Any[])

td0 = :(@traitdef Cr20{X} begin
    length(X)

    @constraints begin
        string(X.name)[1]=='I'
    end
end)
a,b = Traits.parsebody(td0.args[end])
@test b==:(Bool[(string(X.name))[1] == 'I'])

td1 = :(@traitdef Cr20{X} begin
    length(X) -> Int

    @constraints begin
        string(X.name)[1]=='I'
    end
end)
a,b = Traits.parsebody(td1.args[end])
@test b==:(Bool[(string(X.name))[1] == 'I'])

td2 = :(@traitdef Cr20{X,Y} begin
    X + Y -> Int,Float64
    -(X,Y) -> Int
    (/)(X,Y) -> Int

    @constraints begin
        string(X.name)[1]=='I'
    end
end)
a,b,c = Traits.parsebody(td2.args[end])
@test b==:(Bool[(string(X.name))[1] == 'I'])
@test c.head==:block

td3 = :(@traitdef Cr20{X,Y} begin
    fn(X) -> Type{X}
end)
a,b,c = Traits.parsebody(td3.args[end])

# td4 = :(@traitdef Cr20{X} begin
#     fn{Y<:II}(X,Y) -> Type{X}
#     fn76{K<:AbstractFloat, I<:Integer}(X, Vector{I}, Vector{K}) -> I
# end)
# a,b,c = Traits.parsebody(td4.args[end])
# v = :(TypeVar(symbol("Y"),II))
# t = :(TypeVar(symbol("I"),Integer))
# k = :(TypeVar(symbol("K"),AbstractFloat))

# @test a==Expr(:dict, :(fn=>((X,$v),Type{X})),
#                      :(fn76=>((X,Vector{$t},Vector{$k}),$t))
#               )

## test making traits

@traitdef MyIter{X}  begin
    start(X)
end

## Testing trait definitions in commontraits.jl
@test istrait(Cmp{Int,Int})
@test istrait(Cmp{Int,Float64})
@test !istrait(Cmp{Int,AbstractString})

#coll = [Vector, Vector{Int}, Dict{Int}, Dict{Int,Int}, Set{Int}]
coll = [Vector{Int}, Dict{Int,Int}, Set{Int}]
iter = [Traits.GenerateTypeVars{:upcase},  Int] #todo: add AbstractString,
dicts = [Dict{Int,Int}] # Dict and Dict{Int} does not work, ObjectIdDict does not fulfill the trait
index = [Array{Int,2}, StepRange{Int,Int}]
c=1
for c in coll
#    @show IsCollection{c}() # heisenbug protection
    @test istrait(IsCollection{c}, verbose=verbose)
    @test istrait(IsIterable{c}, verbose=verbose)
    @test istrait(IsIterColl{c}, verbose=verbose)
end
@test !istrait(IsIndexable{Set})

for c in iter
    @test istrait(IsIterable{c}, verbose=verbose)
end

for c in dicts
    @test istrait(IsAssociative{c}, verbose=verbose)
end
#
for c in index
    @test istrait(IsIndexable{c}, verbose=verbose)
end

@test istrait(IsIterable{Array}, verbose=verbose)
@test istrait(IsIterable{ASCIIString}, verbose=verbose)
@test istrait(IsIterable{Int}, verbose=verbose)
@test !istrait(IsIterable{Void})

arith = [Int, Float64, Rational{Int}]
a1,a2 = 1,1
for a1 in arith
    for a2 in arith
        @test istrait(Arith{a1,a2}, verbose=verbose)
    end
end

## test trait definition
@traitdef FF{X} begin
    f948576()
end
@test !istrait(FF{Int})
f948576() = 1
@test istrait(FF{Int})

@traitdef Tr20{X} begin
    length(X) -> Bool
end
@traitdef Tr21{X} <: Tr20{X} begin
    size(X) -> Bool
end
@traitdef Tr211{X} <: Tr21{X} begin
    size(X) -> Bool
end
@traitdef Tr2111{X} <: Tr211{X} begin
    size(X) -> Bool
end

@traitdef Tr10{X,Y}  begin
    isless(X,Y) -> Bool
end
@traitdef Tr11{X,Y}  <: Tr10{X,Y} begin
   ==(X,Y) -> Bool
end

@traitdef Tr13{X,Y}  <: Tr11{X,Y}, Tr20{X}, Tr21{Y} begin
   ==(X,Y) -> Bool
end

@test traitgetsuper(Tr20)==Tuple{}
@test traitgetsuper(Tr21)==Tuple{Tr20}
@test traitgetsuper(Tr13)==Tuple{Tr11, Tr20, Tr21}

@test issubtrait(Tr21, Tr20)
@test issubtrait(Tr211, Tr20)
@test issubtrait(Tr2111, Tr20)
@test issubtrait(Tr13, Tr11)
@test issubtrait(Tr13, Tr10)
@test issubtrait(Tr13, Tr21)
@test issubtrait(Tr13, Tr20)

@test issubtrait(Tuple{Tr21}, Tuple{Tr20})
@test  issubtrait(Tuple{Tr21,Tr11}, Tuple{Tr20,Tr10})
@test !issubtrait(Tuple{Tr21,Tr11}, Tuple{Tr10,Tr20}) # todo: this should be true, as order shouldn't matter
@test issubtrait(Tuple{Tr11,Tr21}, Tuple{Tr10,Tr20})

@test !issubtrait(Tr21{Int}, Tr20{Float64})
@test !issubtrait(Tuple{Tr21{Int}}, Tuple{Tr20{Float64}})

####
# Test functions parameterized on non-trait parameters.
###
@traitdef Pr0{X} begin
    fn75{Y <: Integer}(X, Y) -> Y
end
fn75{Y <: Integer}(x::UInt8, y::Y) = y+x
## julia> Base.return_types(fn75, (UInt8, TypeVar(:Y,Integer)))
## 1-element Array{Any,1}:
##  Any
# thus:
if return_types_bug1
    @test !istrait(Pr0{UInt8})
else
    @test istrait(Pr0{UInt8})
end
@test !istrait(Pr0{Int8})

@traitdef Pr1{X}  begin
    fn76{I<:Integer}(X, Vector{I}) -> I
end
fn76{I<:Integer}(x::UInt8, v::Vector{I}) = v[x]
@test istrait(Pr1{UInt8})

@traitdef Pr2{X} begin
    fn77{Y<:Number}(X,Y,Y) -> Number
#    fn77{Y}(X)
end
fn77(a::Array,b::Int, c::Float64) = a[1]
@test !istrait(Pr2{Array})
fn77{Y<:Real}(a::Array,b::Y, c::Y) = a[1]
@test !istrait(Pr2{Array})
fn77{Y<:Number}(a::Array,b::Y, c::Y) = a[1]
if return_types_bug1
    @test !istrait(Pr2{Array})
else
    @test istrait(Pr2{Array})
end

#####
# Trait functions parameterized on trait parameters
####
@traitdef Pr3{X} begin
    fn78{T<:X}(T,T)
end
fn78(b::Int, c::Int) = b
@test istrait(Pr3{Int})
fn78(b::Real, c::Real) = b
@test !istrait(Pr3{Real})
fn78{T}(b::T, c::T) = b
@test istrait(Pr3{Real})
@test istrait(Pr3{Any})

@traitdef Pr04{X} begin
    fnpr04{T<:X, S<:Integer}(T,T, S, S)
end
fnpr04(b::Int, c::Int, ::Int, ::Int) = b
@test !istrait(Pr04{Int})
fnpr04{I<:Integer}(b::Int, c::Int, ::I, ::I) = b
@test istrait(Pr04{Int})

@traitdef Pr05{X} begin
    fnpr05{T<:X, S<:Integer}(Dict{T,T}, Dict{S,T})
end
fnpr05{T<:AbstractFloat, S<:Integer}(::Dict{T,T}, ::Dict{S,T}) = 1
@test istrait(Pr05{Float64})

@traitdef Pr06{X} begin
    fnpr06{T<:X, S<:Integer}(Dict{T,S}, Dict{S,T})
end
fnpr06{T<:AbstractFloat, S<:Integer}(::Dict{T,S}, ::Dict{S,T}) = 1
@test istrait(Pr06{Float64})

@traitdef Pr07{X} begin
    fnpr07(X, X, Integer)
end
fnpr07{T<:Integer}(::T, ::T, ::Integer) = 1
@test !istrait(Pr07{Integer}) # not trait because fnpr07(Int8, UInt8, ...) is not callable
@test istrait(Pr07{Int})

# function parameters only one of the methods
@traitdef Pr08{X} begin
    fnpr08(X, Vector{X}, Integer)
end
fnpr08{T<:Integer}(::T, ::Vector{T}, ::Integer) = 1
@test !istrait(Pr08{Integer})
@test istrait(Pr08{Int})

@traitdef Pr10{X} begin
    fnpr10{T<:X}(T, Vector{T}, Integer)
end
fnpr10(::Int, ::Vector{Int}, ::Integer) = 1
@test !istrait(Pr10{Integer})
@test istrait(Pr10{Int})

@traitdef Pr11{X} begin
    fnpr11(Int, Vector{UInt}, X)
end
fnpr11{T<:Integer}(::T, ::Vector{T}, ::Integer) = 1
@test !istrait(Pr11{Integer})
@test !istrait(Pr11{Int})
@test !istrait(Pr11{UInt})

@traitdef Pr12{X} begin
    fnpr12(Int, Vector{UInt}, X)
end
fnpr12{T<:Integer}(::T, ::Vector{T}, ::Integer) = 1
@test !istrait(Pr12{Integer})
@test !istrait(Pr12{Int})

####
# Test constraints
####
@traitdef Cr20{X} begin
    length(X) -> Any

    @constraints begin
        string(X.name)[1]=='I'
    end
end

@test !istrait(Cr20{Float32})
@test istrait(Cr20{Int})

@traitdef Cr21{X,Y} begin
    ==(X,Y)
    @constraints begin
        string(X.name)[1]=='I'
    end
    # probably not recommended but there can be several constraint
    # blocks:
    @constraints begin
        X==Y
    end
end

@test !istrait(Cr21{Float32, Float32})
@test istrait(Cr21{Int, Int})

@traitdef Cr22{X,Y} begin
    @constraints begin
        X==Y
    end
end
@test istrait(Cr22{Float32, Float32})
@test istrait(Cr22{Int, Int})
@test !istrait(Cr22{Int, Float32})

######
# istrait
#####
check_return_types(false)
f12(x::Int) = 1
@traitdef UU{X} begin
    f12(X)
end

@test !istrait(UU{Any})     # ==false: this should be false
@test !istrait(UU{Integer}) # ==false: this should be false
@test !istrait(UU{Int8})    # ==false: this should be false

f13(x::Integer) = 1
@traitdef UU13{X} begin
    f13(X)
end

@test !istrait(UU13{Any})
@test istrait(UU13{Integer})
@test istrait(UU13{Int8})
check_return_types(true)
#####
# Associated types
####
@traitdef Iter2{X} begin
    # type-functions based on return_type:
    State = Base.return_types(start, (X,))[1]  # this is circular but that is ok, as trait needs to be implemented.
    Item =  Base.return_types(next, (X,State))[1][1]

    # interface functions
    start(X) -> State
    next(X, State) -> Item, State
    done(X, State) -> Bool
    # automatically provides:
    # zip, enumerated, in, map, reduce, ...
end
@test istrait(Iter2{Int})
@test istrait(Iter2{Array})

# isbits associated types, issue #9

@traitdef AssocIsBits{X} begin
    T = X.parameters[1]  # a type
    N = X.parameters[2]  # a isbits
    S = X.parameters[3]  # a symbol
    getindex(X, Int) -> T
end
type T3484675{T,N,S} end
Base.getindex(::T3484675, i::Int) = i

AssocIsBits{T3484675{Int,4.5,:a}}()
@test istrait(AssocIsBits{T3484675{Int,4.5,:a}}) # errors because it is assumed that all
                                                 # parameters are TypeVars
#####
# Varags
#####
@traitdef TT31{X} begin
    foo31(X, Int...)
end
foo31(::AbstractString, x::UInt...) = 1
@test !istrait(TT31{AbstractString})
foo31(::AbstractString) = 2 # to avoid ambiguity warnings
foo31(::AbstractString, x::Int...) = 2
@test istrait(TT31{AbstractString})

@traitdef TT32{X} begin
    foo32(X...)
end
foo32(::AbstractString) = 1
@test !istrait(TT32{AbstractString})
foo32(a::AbstractString...) = 2 # to avoid ambiguity warnings
@test istrait(TT32{AbstractString})

@traitdef TT33{X} begin
    foo33{Y<:X}(X, Y...)
end
foo33(::AbstractString) = 1
@test !istrait(TT33{AbstractString})
foo33{T<:AbstractString}(::AbstractString, a::T...) = 2
@test istrait(TT33{AbstractString})

####
# DataType constructors
####

@traitdef TT45{D} begin
    # This trait contains all datatypes which have a constructor with
    # no arguments.
    D() -> D
end
type A4758 end
type A4759
    a
end

@test istrait(TT45{A4758})
@test !istrait(TT45{A4759})
@test istrait(TT45{Dict{Int,Int}})
@test istrait(TT45{Set{Int}})
@test !istrait(TT45{Int})
@test istrait(TT45{Array{Int,1}})

@traitdef TT44{D} begin
    Array(Type{D},Integer) -> Array # the standard array constructor, should be working for all Types
end
@test istrait(TT44{A4758})
@test istrait(TT44{A4759})
@test istrait(TT44{Dict{Int,Int}})
@test istrait(TT44{Set{Int}})
@test istrait(TT44{Int})
@test istrait(TT44{Array{Int,1}})

# This is the trait for datatypes with Array like constructors:
@traitdef TT46{Ar} begin
    T = Type{eltype(Ar)}
    Arnp = deparameterize_type(Ar)  # Array stripped of type parameters

    Arnp(T, Int64) -> Ar
    Arnp(T, Int...) -> Ar
    @constraints begin
        length(Ar.parameters)>1 # need at least two parameters to be array-like, right?
    end
end
@test !istrait(TT46{A4758})
@test !istrait(TT46{Dict{Int,Int}})
# @test istrait(TT46{Set{Int}}, verbose=verbose) this actually works, but not as expected and gives a deprecation warning
@test !istrait(TT46{Int})
@test istrait(TT46{Array{Int,1}}, verbose=verbose)

# TODO: This does not pass currently because of:
## julia> f() = Array{Int}()
## f (generic function with 1 method)
#
## julia> Base.return_types(f, ())
## 1-element Array{Any,1}:
##  Array{Int64,0}
# @test istrait(TT46{Array{Int}}, verbose=verbose)
@test istrait(TT46{Array}, verbose=verbose)


text:
using NeutralLandscapes
using Test

# simply test it's running
@test DiscreteVoronoi(42) != π


text:
push!(LOAD_PATH,"../src/")

using Documenter, Oxygen

makedocs(
    sitename="Oxygen.jl",
    authors="Rushil Mallarapu",
    modules = [Oxygen],
    format = Documenter.HTML(prettyurls = false),
    pages = [
        "Home" => "index.md",
        "API" => [
            "api/atoms.md",
            "api/mols.md",
            "api/smiles.md",
            "api/molops.md",
        ],
        "TODO" => "todo.md"
    ]
)

deploydocs(
    repo = "github.com/sudo-rushil/Oxygen.jl.git",
)


text:
struct FLSynapseParameter
end

@with_kw mutable struct FLSynapse
    param::FLSynapseParameter = FLSynapseParameter()
    W::Matrix{SNNFloat}  # synaptic weight
    rI::Vector{SNNFloat} # postsynaptic rate
    rJ::Vector{SNNFloat} # presynaptic rate
    g::Vector{SNNFloat}  # postsynaptic conductance
    P::Matrix{SNNFloat}  # <rᵢrⱼ>⁻¹
    q::Vector{SNNFloat}  # P * r
    u::Vector{SNNFloat} # force weight
    w::Vector{SNNFloat} # output weight
    f::SNNFloat = 0 # postsynaptic traget
    z::SNNFloat = 0.5randn()  # output z ≈ f
    records::Dict = Dict()
end

function FLSynapse(pre, post; σ = 1.5, p = 0.0, α = 1)
    rI, rJ, g = post.r, pre.r, post.g
    W = σ * 1 / √pre.N * randn(post.N, pre.N) # normalized recurrent weight
    w = 1 / √post.N * (2rand(post.N) - 1) # initial output weight
    u = 2rand(post.N) - 1 # initial force weight
    P = α * eye(post.N) # initial inverse of C = <rr'>
    q = zeros(post.N)
    FLSynapse(;@symdict(W, rI, rJ, g, P, q, u, w)...)
end

function forward!(c::FLSynapse, param::FLSynapseParameter)
    z = dot(w, rI)
    BLAS.A_mul_B!(q, P, rJ)
    BLAS.A_mul_B!(g, W, rJ)
    BLAS.axpy!(z, u, g)
end

function plasticity!(c::FLSynapse, param::FLSynapseParameter, dt::SNNFloat, t::SNNFloat)
    C = 1 / (1 + dot(q, rI))
    BLAS.axpy!(C * (f - z), q, w)
    BLAS.ger!(-C, q, q, P)
end


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

mutable struct VirtualMachineScaleSetExtensionProfile <: SwaggerModel
    extensions::Any # spec type: Union{ Nothing, Vector{VirtualMachineScaleSetExtension} } # spec name: extensions

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

const _property_map_VirtualMachineScaleSetExtensionProfile = Dict{Symbol,Symbol}(Symbol("extensions")=>Symbol("extensions"))
const _property_types_VirtualMachineScaleSetExtensionProfile = Dict{Symbol,String}(Symbol("extensions")=>"Vector{VirtualMachineScaleSetExtension}")
Base.propertynames(::Type{ VirtualMachineScaleSetExtensionProfile }) = collect(keys(_property_map_VirtualMachineScaleSetExtensionProfile))
Swagger.property_type(::Type{ VirtualMachineScaleSetExtensionProfile }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_VirtualMachineScaleSetExtensionProfile[name]))}
Swagger.field_name(::Type{ VirtualMachineScaleSetExtensionProfile }, property_name::Symbol) =  _property_map_VirtualMachineScaleSetExtensionProfile[property_name]

function check_required(o::VirtualMachineScaleSetExtensionProfile)
    true
end

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


text:
""" TODO """
function readShaders(program::AbstractGraphicsShaderProgram, fileSource::FileSource)
	if program == nothing end

	global shaderAttributes = []
	shaderSource=""
	globalSource=""
	lokalSource=""
	shaderType=nothing
	pbegin=0
	comment=false
	count=0

	#o.logShader.Info("Shader %s\n", path)

	clear(program)

	shaderCatch = function()
		pend=length(shaderSource)

		if shaderType != nothing
			f=FileSourcePart(fileSource,range(pbegin+1,pend-pbegin))
			f.cache = globalSource * lokalSource #f.source.cache[f.range]
			shader = createShader(program,shaderType,f)
			shader.cache="" #DEFAULT_SHADER_CODE(shaderType) #default
			lokalSource=""
			shaderType=nothing
		end

		pbegin = pend
		count+=1
		#o.logShader.Info("(%v,%v) %s\n", sourceBegin, sourceEnd, shader.GetSource())
	end

	readLines = function(file)
		for line in eachline(file)
			#line=replace(line, "\r", "") #remove windows extra break sign
			line=rstrip(lstrip(line)) #trim
			#line=filter(x -> !isspace(x), line)

			len = length(line)

			if len < 1 continue end

			firstSign = line[1]
			inputFmt = firstSign

			if len > 1 inputFmt = line[1:2] end

			if comment
				if inputFmt == "*/" comment = false end
				#continue
			#elseif firstSign == '\n'
				#continue
			elseif inputFmt == "/*"
				comment = true
				#continue
			elseif inputFmt == "//" # Comment
				#continue

				if len < 3 continue	end

				if line[3] == '!'
					ln = line[3:len]
					len = length(ln)

					if len < 2 continue	end

					sp = split(ln[2:len], " ")
					len = length(sp)
					cmd = ""
					value = ""

					if len < 1 continue end

					cmd = uppercase(lstrip(sp[1]))

					if len > 1
						value = uppercase(rstrip(sp[2]))

						if cmd == ""
							cmd = value
						end
					end

					#o.logShader.Info("(%s):\n", value)

					if cmd == value
						value=Symbol(value)
						if haskey(LIST_SHADER,value)
							shaderCatch()
							shaderType=value
						end
						#for (k, v) in SHADER_FILE_TYPE
						#	if v == value
						#		shaderCatch()
						#		shaderType=v
						#		break
						#	end
						#end
					end
				end
				#continue

			else
				findShaderProperty(program, line)
				if count == 0
					globalSource = string(globalSource,line,"\n")
				else
					lokalSource = string(lokalSource,line,"\n")
				end
				shaderSource = string(shaderSource,line,"\n")
			end
		end
		shaderCatch()
		fileSource.cache=shaderSource
		shaderSource != ""
	end

	FileManager.waitForFileReady(fileSource.path,readLines)
end

include("GLSLRessources.jl")

max_stars_repo_name: Gilga/JGE

max_issues_repo_name: Gilga/JGE

max_forks_repo_name: Gilga/JGE


text:
"""
    correct_reduced_field_size(loc, arch, grid, dims, Tx, Ty, Tz)

Test that the ReducedField at `loc`ation on `arch`itecture and `grid`
and reduced along `dims` has size `(Tx, Ty, Tz)`.
"""
correct_reduced_field_size(loc, arch, grid, dims, Tx, Ty, Tz) =
    size(parent(ReducedField(loc, arch, grid; dims=dims))) == (Tx, Ty, Tz)

function correct_reduced_field_value_was_set(arch, grid, loc, dims, val::Number)
    f = ReducedField(loc, arch, grid; dims=dims)
    set!(f, val)
    return all(interior(f) .≈ val * arch_array(arch, ones(size(f))))
end

@testset "ReducedFields" begin
    @info "Testing ReducedFields..."

    N = (4, 6, 8)
    L = (2π, 3π, 5π)
    H = (1, 1, 1)

    @testset "ReducedField initialization" begin
        @info "  Testing ReducedField initialization..."
        for arch in archs, FT in float_types

            grid = RegularRectilinearGrid(FT, size=N, extent=L, halo=H, topology=(Bounded, Bounded, Bounded))

            @test correct_reduced_field_size((Center, Center, Center), arch, grid, 1,         1,               N[2] + 2 * H[2],     N[3] + 2 * H[3])
            @test correct_reduced_field_size((Face,   Center, Center), arch, grid, 1,         1,               N[2] + 2 * H[2],     N[3] + 2 * H[3])
            @test correct_reduced_field_size((Center, Face,   Center), arch, grid, 1,         1,               N[2] + 2 * H[2] + 1, N[3] + 2 * H[3])
            @test correct_reduced_field_size((Center, Face,   Face),   arch, grid, 1,         1,               N[2] + 2 * H[2] + 1, N[3] + 2 * H[3] + 1)
            @test correct_reduced_field_size((Center, Center, Center), arch, grid, 2,         N[1] + 2 * H[1], 1,                   N[3] + 2 * H[3])
            @test correct_reduced_field_size((Center, Center, Center), arch, grid, 2,         N[1] + 2 * H[1], 1,                   N[3] + 2 * H[3])
            @test correct_reduced_field_size((Center, Center, Center), arch, grid, 3,         N[1] + 2 * H[1], N[2] + 2 * H[2],     1)
            @test correct_reduced_field_size((Center, Center, Center), arch, grid, (1, 2),    1,               1,                   N[3] + 2 * H[3])
            @test correct_reduced_field_size((Center, Center, Center), arch, grid, (2, 3),    N[1] + 2 * H[1], 1,                   1)
            @test correct_reduced_field_size((Center, Center, Center), arch, grid, (1, 2, 3), 1,               1,                   1)
        end
    end

    reduced_dims = (1, 2, 3, (1, 2), (2, 3), (1, 3), (1, 2, 3))

    int_vals = Any[0, Int8(-1), Int16(2), Int32(-3), Int64(4), Int128(-5)]
    uint_vals = Any[6, UInt8(7), UInt16(8), UInt32(9), UInt64(10), UInt128(11)]
    float_vals = Any[0.0, -0.0, 6e-34, 1.0f10]
    rational_vals = Any[1//11, -23//7]
    other_vals = Any[π]
    vals = vcat(int_vals, uint_vals, float_vals, rational_vals, other_vals)

    @testset "Setting ReducedFields" begin
        @info "  Testing ReducedField setting..."
        for arch in archs, FT in float_types

            grid = RegularRectilinearGrid(FT, size=N, extent=L, halo=H, topology=(Periodic, Periodic, Bounded))

            for dims in reduced_dims, val in vals
                @test correct_reduced_field_value_was_set(arch, grid, (Center, Center, Center), dims, val)
            end

            for dims in reduced_dims
                reduced_field = ReducedField((Center, Center, Center), arch, grid, dims=dims)
                sz = size(reduced_field)
                A = rand(FT, sz...)
                set!(reduced_field, A)

                @test reduced_field[1, 1, 1] == A[1, 1, 1]

                fill_halo_regions!(reduced_field, arch)

                # No-flux boundary conditions at top and bottom
                @test reduced_field[1, 1, 0] == A[1, 1, 1]
                @test reduced_field[1, 1, grid.Nz+1] == A[1, 1, end]

                # Periodic boundary conditions in the horizontal directions
                @test reduced_field[1, 0, 1] == A[1, end, 1]
                @test reduced_field[1, grid.Ny+1, 1] == A[1, 1, 1]

                @test reduced_field[0, 1, 1] == A[end, 1, 1]
                @test reduced_field[grid.Nx+1, 1, 1] == A[1, 1, 1]
            end
        end
    end
end


text:
"Autumn Language"
module Autumn
using Reexport

include("parameters.jl")
@reexport using .Parameters

include("aexpr.jl")
@reexport using .AExpressions

include("util.jl")
@reexport using .Util

include("subexpr.jl")
@reexport using .SubExpressions

include("sexpr.jl")
@reexport using .SExpr

include("compileutils.jl")
@reexport using .CompileUtils

include("compile.jl")
@reexport using .Compile

include("abstractinterpretation.jl")
@reexport using .AbstractInterpretation

include("transform.jl")
@reexport using .Transform

include("best_effort_ranking.jl")
@reexport using .BestEffortRanking
end


text:
@inline sqr(x) = x*x
@inline cub(x) = x*x*x