text:
function distance_weigth(a::NTuple{N,T}, b::NTuple{N,T}; order=4) where {N,T}
    return 1 / distance(a, b)^order
end

## CPU 2D

@inbounds function _gathering!(upper, lower, Fpi, p, x, y, dxi, order)

    # check that the particle is inside the grid
    isinside(p, x, y)

    # indices of lowermost-left corner of
    # the cell containing the particlex
    idx_x, idx_y = parent_cell(p, dxi)

    ω = (
        distance_weigth((x[idx_x], y[idx_y]), p; order=order),
        distance_weigth((x[idx_x + 1], y[idx_y]), p; order=order),
        distance_weigth((x[idx_x], y[idx_y + 1]), p; order=order),
        distance_weigth((x[idx_x + 1], y[idx_y + 1]), p; order=order),
    )

    nt = Threads.threadid()

    upper[nt][idx_x, idx_y] += ω[1] * Fpi
    upper[nt][idx_x + 1, idx_y] += ω[2] * Fpi
    upper[nt][idx_x, idx_y + 1] += ω[3] * Fpi
    upper[nt][idx_x + 1, idx_y + 1] += ω[4] * Fpi
    lower[nt][idx_x, idx_y] += ω[1]
    lower[nt][idx_x + 1, idx_y] += ω[2]
    lower[nt][idx_x, idx_y + 1] += ω[3]

    return lower[nt][idx_x + 1, idx_y + 1] += ω[4]
end

function gathering!(F::Array{T,2}, Fp::Vector{T}, xi, particle_coords; order=4) where {T}

    # unpack tuples
    px, py = particle_coords
    x, y = xi
    dxi = (x[2] - x[1], y[2] - y[1])

    # number of particles
    np = length(Fp)

    # TODO think about pre-allocating these 2 buffers
    upper = [zeros(size(F)) for _ in 1:Threads.nthreads()]
    lower = [zeros(size(F)) for _ in 1:Threads.nthreads()]

    # compute ∑ωᵢFᵢ and ∑ωᵢ
    Threads.@threads for i in 1:np
        _gathering!(upper, lower, Fp[i], (px[i], py[i]), x, y, dxi, order)
    end

    # compute Fᵢ=∑ωᵢFpᵢ/∑ωᵢ
    Threads.@threads for i in eachindex(F)
        @inbounds F[i] =
            sum(upper[nt][i] for nt in 1:Threads.nthreads()) /
            sum(lower[nt][i] for nt in 1:Threads.nthreads())
    end
end

## CPU 3D

@inbounds function _gathering!(upper, lower, Fpi, p, x, y, z, dxi, order)
    # check that the particle is inside the grid
    isinside(p, x, y, z)

    # indices of lowermost-left corner of
    # the cell containing the particle
    idx_x, idx_y, idx_z = parent_cell(p, dxi)

    ω = (
        distance_weigth((x[idx_x], y[idx_y], z[idx_z]), p; order=order),
        distance_weigth((x[idx_x + 1], y[idx_y], z[idx_z]), p; order=order),
        distance_weigth((x[idx_x], y[idx_y + 1], z[idx_z]), p; order=order),
        distance_weigth((x[idx_x + 1], y[idx_y + 1], z[idx_z]), p; order=order),
        distance_weigth((x[idx_x], y[idx_y], z[idx_z + 1]), p; order=order),
        distance_weigth((x[idx_x + 1], y[idx_y], z[idx_z + 1]), p; order=order),
        distance_weigth((x[idx_x], y[idx_y + 1], z[idx_z + 1]), p; order=order),
        distance_weigth((x[idx_x + 1], y[idx_y + 1], z[idx_z + 1]), p; order=order),
    )

    nt = Threads.threadid()

    upper[nt][idx_x, idx_y, idx_z] += ω[1] * Fpi
    upper[nt][idx_x + 1, idx_y, idx_z] += ω[2] * Fpi
    upper[nt][idx_x, idx_y + 1, idx_z] += ω[3] * Fpi
    upper[nt][idx_x + 1, idx_y + 1, idx_z] += ω[4] * Fpi
    upper[nt][idx_x, idx_y, idx_z + 1] += ω[5] * Fpi
    upper[nt][idx_x + 1, idx_y, idx_z + 1] += ω[6] * Fpi
    upper[nt][idx_x, idx_y + 1, idx_z + 1] += ω[7] * Fpi
    upper[nt][idx_x + 1, idx_y + 1, idx_z + 1] += ω[8] * Fpi
    lower[nt][idx_x, idx_y, idx_z] += ω[1]
    lower[nt][idx_x + 1, idx_y, idx_z] += ω[2]
    lower[nt][idx_x, idx_y + 1, idx_z] += ω[3]
    lower[nt][idx_x + 1, idx_y + 1, idx_z] += ω[4]
    lower[nt][idx_x, idx_y, idx_z + 1] += ω[5]
    lower[nt][idx_x + 1, idx_y, idx_z + 1] += ω[6]
    lower[nt][idx_x, idx_y + 1, idx_z + 1] += ω[7]

    return lower[nt][idx_x + 1, idx_y + 1, idx_z + 1] += ω[8]
end

function gathering!(F::Array{T,3}, Fp::Vector{T}, xi, particle_coords; order=4) where {T}
    # unpack tuples
    px, py, pz = particle_coords
    x, y, z = xi
    dxi = (x[2] - x[1], y[2] - y[1], z[2] - z[1])

    # number of particles
    np = length(Fp)

    # TODO think about pre-allocating these 2 buffers
    upper = [zeros(size(F)) for _ in 1:Threads.nthreads()]
    lower = [zeros(size(F)) for _ in 1:Threads.nthreads()]

    # compute ∑ωᵢFᵢ and ∑ωᵢ
    Threads.@threads for i in 1:np
        _gathering!(upper, lower, Fp[i], (px[i], py[i], pz[i]), x, y, z, dxi, order)
    end

    # compute Fᵢ=∑ωᵢFpᵢ/∑ωᵢ
    Threads.@threads for i in eachindex(F)
        @inbounds F[i] =
            sum(upper[nt][i] for nt in 1:Threads.nthreads()) /
            sum(lower[nt][i] for nt in 1:Threads.nthreads())
    end
end

## CUDA 2D

function _gather1!(
    upper::CuDeviceMatrix{T,1},
    lower::CuDeviceMatrix{T,1},
    Fpd::CuDeviceVector{T,1},
    xi,
    dxi,
    p;
    order=4,
) where {T}
    idx = (blockIdx().x - 1) * blockDim().x + threadIdx().x

    # unpack tuples
    px, py = p
    x, y = xi

    @inbounds if idx ≤ length(px)
        # check that the particle is inside the grid
        isinside(px[idx], py[idx], x, y)

        p_idx = (px[idx], py[idx])

        # indices of lowermost-left corner of
        # the cell containing the particle
        idx_x, idx_y = parent_cell(p_idx, dxi)

        ω1::Float64 = distance_weigth((x[idx_x], y[idx_y]), p_idx; order=order)
        ω2::Float64 = distance_weigth((x[idx_x + 1], y[idx_y]), p_idx; order=order)
        ω3::Float64 = distance_weigth((x[idx_x], y[idx_y + 1]), p_idx; order=order)
        ω4::Float64 = distance_weigth((x[idx_x + 1], y[idx_y + 1]), p_idx; order=order)

        Fpi::Float64 = Fpd[idx]

        CUDA.@atomic upper[idx_x, idx_y] += ω1 * Fpi
        CUDA.@atomic upper[idx_x + 1, idx_y] += ω2 * Fpi
        CUDA.@atomic upper[idx_x, idx_y + 1] += ω3 * Fpi
        CUDA.@atomic upper[idx_x + 1, idx_y + 1] += ω4 * Fpi
        CUDA.@atomic lower[idx_x, idx_y] += ω1
        CUDA.@atomic lower[idx_x + 1, idx_y] += ω2
        CUDA.@atomic lower[idx_x, idx_y + 1] += ω3
        CUDA.@atomic lower[idx_x + 1, idx_y + 1] += ω4
    end

    return nothing
end

function _gather2!(Fd::CuDeviceArray{T,2}, upper, lower) where {T}
    idx = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    idy = (blockIdx().y - 1) * blockDim().y + threadIdx().y

    if (idx < size(Fd, 1)) && (idy < size(Fd, 2))
        @inbounds Fd[idx, idy] = upper[idx, idy] / lower[idx, idy]
    end

    return nothing
end

function gathering!(
    Fd::CuArray{T,2}, Fpd::CuArray{T,1}, xi, dxi, particle_coords; nt=512
) where {T}
    upper = CUDA.zeros(T, size(Fd))
    lower = CUDA.zeros(T, size(Fd))
    dxi = (x[2] - x[1], y[2] - y[1])

    # first kernel that computes ∑ωᵢFᵢ and ∑ωᵢ
    N = length(Fpd)
    numblocks = ceil(Int, N / nt)
    CUDA.@sync begin
        @cuda threads = nt blocks = numblocks _gather1!(
            upper, lower, Fpd, xi, dxi, particle_coords
        )
    end

    # seond and final kernel that computes Fᵢ=∑ωᵢFpᵢ/∑ωᵢ
    nx, ny = size(Fd)
    nblocksx = ceil(Int, nx / 32)
    nblocksy = ceil(Int, ny / 32)
    CUDA.@sync begin
        @cuda threads = (32, 32) blocks = (nblocksx, nblocksy, nblocksz) _gather2!(
            Fd, upper, lower
        )
    end
end

## CUDA 3D

function _gather1!(
    upper::CuDeviceArray{T,3},
    lower::CuDeviceArray{T,3},
    Fpd::CuDeviceVector{T,1},
    xi,
    dxi,
    p;
    order=4,
) where {T}
    idx = (blockIdx().x - 1) * blockDim().x + threadIdx().x

    # unpack tuples
    px, py, pz = p
    x, y, z = xi

    @inbounds if idx ≤ length(px)
        # check that the particle is inside the grid
        isinside(px[idx], py[idx], pz[idx], x, y, z)

        p_idx = (px[idx], py[idx], pz[idx])

        # indices of lowermost-left corner of
        # the cell containing the particle
        idx_x, idx_y, idx_z = parent_cell(p_idx, dxi)

        ω1::Float64 = distance_weigth((x[idx_x], y[idx_y], z[idx_z]), p_idx; order=order)
        ω2::Float64 = distance_weigth(
            (x[idx_x + 1], y[idx_y], z[idx_z]), p_idx; order=order
        )
        ω3::Float64 = distance_weigth(
            (x[idx_x], y[idx_y + 1], z[idx_z]), p_idx; order=order
        )
        ω4::Float64 = distance_weigth(
            (x[idx_x + 1], y[idx_y + 1], z[idx_z]), p_idx; order=order
        )
        ω5::Float64 = distance_weigth(
            (x[idx_x], y[idx_y], z[idx_z + 1]), p_idx; order=order
        )
        ω6::Float64 = distance_weigth(
            (x[idx_x + 1], y[idx_y], z[idx_z + 1]), p_idx; order=order
        )
        ω7::Float64 = distance_weigth(
            (x[idx_x], y[idx_y + 1], z[idx_z + 1]), p_idx; order=order
        )
        ω8::Float64 = distance_weigth(
            (x[idx_x + 1], y[idx_y + 1], z[idx_z + 1]), p_idx; order=order
        )

        Fpi::Float64 = Fpd[idx]

        CUDA.@atomic upper[idx_x, idx_y, idx_z] += ω1 * Fpi
        CUDA.@atomic upper[idx_x + 1, idx_y, idx_z] += ω2 * Fpi
        CUDA.@atomic upper[idx_x, idx_y + 1, idx_z] += ω3 * Fpi
        CUDA.@atomic upper[idx_x + 1, idx_y + 1, idx_z] += ω4 * Fpi
        CUDA.@atomic upper[idx_x, idx_y, idx_z + 1] += ω5 * Fpi
        CUDA.@atomic upper[idx_x + 1, idx_y, idx_z + 1] += ω6 * Fpi
        CUDA.@atomic upper[idx_x, idx_y + 1, idx_z + 1] += ω7 * Fpi
        CUDA.@atomic upper[idx_x + 1, idx_y + 1, idx_z + 1] += ω8 * Fpi
        CUDA.@atomic lower[idx_x, idx_y, idx_z] += ω1
        CUDA.@atomic lower[idx_x + 1, idx_y, idx_z] += ω2
        CUDA.@atomic lower[idx_x, idx_y + 1, idx_z] += ω3
        CUDA.@atomic lower[idx_x + 1, idx_y + 1, idx_z] += ω4
        CUDA.@atomic lower[idx_x, idx_y, idx_z + 1] += ω5
        CUDA.@atomic lower[idx_x + 1, idx_y, idx_z + 1] += ω6
        CUDA.@atomic lower[idx_x, idx_y + 1, idx_z + 1] += ω7
        CUDA.@atomic lower[idx_x + 1, idx_y + 1, idx_z + 1] += ω8
    end

    return nothing
end

function _gather2!(
    Fd::CuDeviceArray{T,3}, upper::CuDeviceArray{T,3}, lower::CuDeviceArray{T,3}
) where {T}
    idx = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    idy = (blockIdx().y - 1) * blockDim().y + threadIdx().y
    idz = (blockIdx().z - 1) * blockDim().z + threadIdx().z

    if (idx ≤ size(Fd, 1)) && (idy ≤ size(Fd, 2)) && (idz ≤ size(Fd, 3))
        @inbounds Fd[idx, idy, idz] = upper[idx, idy, idz] / lower[idx, idy, idz]
    end

    return nothing
end

function gathering!(
    Fd::CuArray{T,3}, Fpd::CuArray{T,1}, xi, particle_coords; nt=512
) where {T}
    x, y, z = xi
    dxi = (x[2] - x[1], y[2] - y[1], z[2] - z[1])
    upper = CUDA.zeros(T, size(Fd))
    lower = CUDA.zeros(T, size(Fd))

    # first kernel that computes ∑ωᵢFᵢ and ∑ωᵢ
    N = length(Fpd)
    numblocks = ceil(Int, N / nt)
    CUDA.@sync begin
        @cuda threads = nt blocks = numblocks _gather1!(
            upper, lower, Fpd, xi, dxi, particle_coords
        )
    end

    # second and final kernel that computes Fᵢ=∑ωᵢFpᵢ/∑ωᵢ
    nx, ny, nz = size(Fd)
    nblocksx = ceil(Int, nx / 8)
    nblocksy = ceil(Int, ny / 8)
    nblocksz = ceil(Int, nz / 8)
    CUDA.@sync begin
        @cuda threads = (8, 8, 8) blocks = (nblocksx, nblocksy, nblocksz) _gather2!(
            Fd, upper, lower
        )
    end
end


text:
# Correia, Sérgio P.C. and Dickinson, Anthony and Clayton, Nicola S., Western
# Scrub-Jays Anticipate Future Needs Independently of Their Current Motivational
# State, Current Biology 17 (10), 856-861, http://dx.doi.org/10.1016/j.cub.2007.03.063
#
function _summarize(::Experiment{:Correia07_exp1}, results)
	combine(groupby(results, [:prefed, :action, :foodtype]),
				 df -> DataFrame(μ = mean(df.count), sem = sem(df.count)))
end
function statistical_tests(exp::Experiment{:Correia07_exp1}, data)
    groups = groupby(data, [:prefed])
    tests = Test[]
    id = "interaction"
    push!(tests, Test(id,
                      splitplotanova(data,
                                     exp.tests[id].locals,
                                     betweenfactors = [],
                                     withinfactors = [:prefed, :foodtype],
                                     id = :id, data = :count),
                      exp.tests))
    for g in groups
        if g.prefed[1] in ("both", "none")
            push!(tests, Test("prefed $(g.prefed[1])",
                              isnot_significantly_different(g,
                                                            conditionkey = :foodtype,
                                                            condition1 = "kibble",
                                                            condition2 = "pinenut",
                                                            valuekey = :count),
                              exp.tests))
        else
            conditions = ["kibble", "pinenut"]
            if g.prefed[1] == "kibble"
                reverse!(conditions)
            end
            push!(tests, Test("prefed $(g.prefed[1])",
                              is_significantly_bigger(g, conditionkey = :foodtype,
                                                      condbigger = conditions[1],
                                                      condsmaller = conditions[2],
                                                      valuekey = :count),
                              exp.tests))
        end
    end
    TestCollection(tests)
end
function run!(exp::Experiment{:Correia07_exp1}, models; N = 11)
    results = DataFrame(prefed = String[], action = String[],
                        foodtype = String[], id = Int[],
						count = Int[])
	for id in 1:N
		m = models[id]
		for prefed in ("pinenut", "kibble", "both", "none")
			if prefed in ("pinenut", "both"); add!(m, Pinenut, 100); end
			if prefed in ("kibble", "both"); add!(m, Kibble, 100); end
			wait!(m, 3u"hr")
			remove!(m, Any)
			add!(m, Pinenut, 30)
			add!(m, Kibble, 30)
			wait!(m, 10u"minute")
			push!(results, [prefed, "eat", "pinenut", id, 30 - countfooditems(m, Pinenut)])
			push!(results, [prefed, "eat", "kibble", id, 30 - countfooditems(m, Kibble)])
			remove!(m, Any)
			add!(m, MaintenanceDiet)
			wait!(m, 50u"minute" + 20u"hr") # next day; not specified in paper
			remove!(m, MaintenanceDiet)
		end
	end
	results
end

function _summarize(::Experiment{:Correia07_exp2}, results)
	combine(groupby(results, [:trial, :group, :foodtype, :action]),
				 df -> DataFrame(μ = mean(df.count), sem = sem(df.count)))
end
function statistical_tests(exp::Experiment{:Correia07_exp2}, data)
    groups = groupby(data, :action)
    tests = Test[]
    for g in groups
        if g.action[1] == "eat"
            id = "eat"
            test = Test(id,
                        splitplotanova(g,
                                       exp.tests[id].locals,
                                       betweenfactors = :group,
                                       withinfactors = [:foodtype, :trial],
                                       id = :id, data = :count),
                        exp.tests)
            push!(tests, test)
        else
            id = "cache"
            test = Test(id,
                        splitplotanova(g,
                                       exp.tests[id].locals,
                                       betweenfactors = [:group, :prefed],
                                       withinfactors = [:foodtype, :trial],
                                       id = :id, data = :count),
                        exp.tests)
            push!(tests, test)
            ddata = combine(groupby(g, [:trial, :group, :action, :id]),
                       df ->  begin
#                            df.foodtype[1] != df.foodtype[2] &&
#                            length(df.foodtype) == 2 || error(df)
                           sum(df.count) > 0 ? df.count[df.foodtype[1] == "prefed" ? 1 : 2]/sum(df.count) : missing end)
            dropmissing!(ddata, disallowmissing = true)
            ddata.x1 = float.(ddata.x1)
            for i in 1:3 # trials
                x = @where ddata :trial .== i :group .== "same"
                y = @where ddata :trial .== i :group .== "different"
                if length(x.x1) > 0 && length(y.x1) > 0
                    test = Test("trial $i: proportion cache",
                                utest(x.x1, y.x1),
                                exp.tests)
                    push!(tests, test)
                elseif length(y.x1) > 0
                    push!(tests, Test("trial $i: proportion cache",
                                      NoTest("no caching in same condition"),
                                      exp.tests,
                                      [1.]))
                else
                    push!(tests, Test("trial $i: proportion cache",
                                      NoTest("no caching at all"),
                                      exp.tests,
                                      [1.]))
                end
            end
        end
    end
    TestCollection(tests)
end
function run!(exp::Experiment{:Correia07_exp2}, models; N = 11)
	results = DataFrame(trial = Int[], group = String[],
                        prefed = String[], foodtype = String[],
						action = String[], id = Int[], count = Int[])
	for id in 1:N
		m = models[id]
		for trial in 1:3
		    tray = Tray()
            if id % 2 == 0
                prefed = Pinenut
            else
                prefed = Kibble
            end
            add!(m, prefed, 100)
            wait!(m, 3u"hr")
            remove!(m, Any)
            add!(m, Pinenut, 30)
            add!(m, Kibble, 30)
            add!(m, tray)
            wait!(m, 10u"minute")
            group = id > N/2 ? "different" : "same"
            for foodtype in (Pinenut, Kibble)
                food = foodtype == prefed ? "prefed" : "nonprefed"
                push!(results, [trial, group, string(prefed), food,
                                "eat", id, 30 - countfooditems(m, foodtype)])
                push!(results, [trial, group, string(prefed), food,
                                "cache", id, countcache(tray, foodtype)])
            end
            remove!(m, Any)
            wait!(m, 30u"minute")
            if group == "same"
                add!(m, prefed, 100)
            else
                add!(m, setdiff([Pinenut, Kibble], [prefed])[1], 100)
            end
            wait!(m, 3u"hr")
            remove!(m, Any)
            add!(m, tray)
            wait!(m, 10u"minute")
            remove!(m, Any)
            add!(m, MaintenanceDiet)
            wait!(m, 10u"minute" + 16u"hr") # next day; not specified in paper
            remove!(m, MaintenanceDiet)
		end
	end
	results
end

ext: jl


text:
abstract type AbstractEvaluator end

mutable struct SameInstancesEvaluator <: AbstractEvaluator
    instances::Union{Array{CPModel}, Nothing}
    metrics::Union{Matrix{<:AbstractMetrics}, Nothing}
    evalFreq::Int64
    nbInstances::Int64
    nbHeuristics::Union{Int64, Nothing}
end

function SameInstancesEvaluator(valueSelectionArray::Array{H, 1}, generator::AbstractModelGenerator; seed=nothing, evalFreq::Int64 = 50, nbInstances::Int64 = 10, evalTimeOut::Union{Nothing,Int64} = nothing) where H<: ValueSelection
    instances = Array{CPModel}(undef, nbInstances)
    metrics = Matrix{AbstractMetrics}(undef,nbInstances, size(valueSelectionArray,1))

    for i in 1:nbInstances
        instances[i] = CPModel()
        fill_with_generator!(instances[i], generator; seed=seed)
        instances[i].limit.searchingTime = evalTimeOut
        for (j, value) in enumerate(valueSelectionArray)
            metrics[i,j]=BasicMetrics(instances[i],value;meanOver=1)   #no slidding mean on evaluation because the instance remains the same
        end
    end
    SameInstancesEvaluator(instances, metrics, evalFreq, nbInstances, size(valueSelectionArray,1))
end

function evaluate(eval::SameInstancesEvaluator, variableHeuristic::AbstractVariableSelection, strategy::S; verbose::Bool=true) where{S<:SearchStrategy}
    for j in 1:eval.nbHeuristics
        heuristic = eval.metrics[1,j].heuristic
        initsize = isa(heuristic, LearnedHeuristic) ? length(heuristic.agent.trajectory) : nothing

        testmode!(heuristic, true)
        for i in 1:eval.nbInstances
            model = eval.instances[i]
            reset_model!(model)

            dt = @elapsed search!(model, strategy, variableHeuristic, heuristic)
            eval.metrics[i,j](model,dt)

            verbose && println(typeof(heuristic), " evaluated with: ", model.statistics.numberOfNodes, " nodes, taken ", dt, "s, number of solutions found : ", model.statistics.numberOfSolutions)
        end
        testmode!(heuristic, false)
        if !isnothing(initsize)
            @assert length(heuristic.agent.trajectory) == initsize "You have leaks in your evaluation pipeline!"
        end
    end
end


text:
# Generate a random sparse matrix using SuiteSparseGraphBLAS Julia wrapper
# from https://github.com/abhinavmehndiratta/SuiteSparseGraphBLAS.jl

using GraphBLASInterface, SuiteSparseGraphBLAS

"""
    GraphBLAS_sparse_square_matrix!(A::GrB_Matrix{T}, dim::Int) where T

Build a random sparse square matrix of dimension `dim`x`dim` using GraphBLAS.
"""
function GraphBLAS_sparse_square_matrix!(
    A::GrB_Matrix{T},
    dim::Int) where T

    # Build A with at most sqrt(N) elements, N being the number of nonzero matrix elements
    ntuples = dim

    for k in 1:ntuples

        i = ZeroBasedIndex(rand(0:dim-1))
        j = ZeroBasedIndex(rand(0:dim-1))
        x = rand(T)

        GrB_Matrix_setElement(A, x, i, j)
    end

    # Force pending operations to complete and return
    GrB_wait()

    return (GrB_SUCCESS)
end


text:
module BioinformaticsBISC195

export normalizeDNA,
       composition,
       gc_content,
       complement,
       reverse_complement,
       parse_fasta

# # uncomment the following line if you intend to use BioSequences types
 using BioSequences
 import BioSequences: composition, gc_content, complement, reverse_complement

"""
    normalizeDNA(::AbstractString)

Ensures that a sequence only contains valid bases
(or `'N'` for unknown bases).
Returns a String.
"""
function normalizeDNA(seq)
    seq = uppercase(string(seq))
    for base in seq
        # note: `N` indicates an unknown base
        occursin(base, "AGCTN") || error("invalid base $base")
    end
    return seq # change to `return LongDNASeq(seq)` if you want to try to use BioSequences types
end

# Your code here.
# Don't forget to export your functions!

function composition(sequence::AbstractString)
    BioSequences.composition(normalizeDNA(sequence))
end

function gc_content(sequence::AbstractString)
    BioSequences.gc_content(normalizeDNA(sequence))
end

function complement(sequence::AbstractString)
    BioSequences.complement(normalizeDNA(sequence))
end

function reverse_complement(sequence::AbstractString)
    BioSequences.reverse_complement(normalizeDNA(sequence))
end

function parse_fasta(path)
    headers = []
    sequences = []
    str = ""
    for line in eachline(path)
        if '>' ∈ line
            push!(headers, string(lstrip(line, '>')))
            if str != ""
                push!(sequences, LongDNASeq(str))
            end
            str = ""
        elseif line != ""
            for base in line
                occursin(base, "AGCTNagctn") || error("invalid base $base")
            end
            str = str * line
        end
    end
    push!(sequences, LongDNASeq(str))
    parsedfile = (headers, sequences)
    return parsedfile
end

end # module Assignment07


text:
abstract AbstractMeasure{T}

abstract SymmetricMeasure{T} <: AbstractMeasure{T}

# Entropy
# =======

"""
Shannon entropy (H)
"""
immutable Entropy{T} <: SymmetricMeasure{T}
  base::T
end

call{T}(::Type{Entropy{T}}) = Entropy(T(Base.e))

## Estimate Entropy using ResidueProbability

"""
`estimate(Entropy(base), p)`

`p` should be a `ResidueProbability` table. The result type is determined by `base`.
"""
function estimate{B, T, N, UseGap}(measure::Entropy{B}, p::ResidueProbability{T, N, UseGap})
  H = zero(B)
  for i in 1:length(p)
    @inbounds pi = B(p[i])
    if pi != 0.0
      H += pi * log(pi)
    end
  end
  -H/log(measure.base)
end

"""
`estimate_on_marginal(Entropy{T}(base), p, marginal)`

This function estimate the entropy H(X) if marginal is 1, H(Y) for 2, etc.
The result type is determined by `base`.
"""
function estimate_on_marginal{B, T, N, UseGap}(measure::Entropy{B}, p::ResidueProbability{T, N, UseGap}, marginal::Int)
  H = zero(B)
  for i in 1:nresidues(p)
    @inbounds pi = B(p.marginals[i, marginal])
    if pi != 0.0
      H += pi * log(pi)
    end
  end
  -H/log(measure.base)
end

## Estimate Entropy using ResidueCount

"""
`estimate(Entropy{T}(base), n::ResidueCount)`

It's the fastest option (you don't spend time on probability calculations).
The result type is determined by the `base`.
"""
function estimate{T}(measure::Entropy{T}, n::ResidueCount)
  H = zero(T)
  total = T(n.total)
  for i in 1:length(n)
    @inbounds ni = T(n[i])
    if ni != 0.0
      H += ni * log(ni/total)
    end
  end
  (-H/total)/log(measure.base)
end

function estimate_on_marginal{T}(measure::Entropy{T}, n::ResidueCount, marginal::Int)
  H = zero(T)
  total = T(n.total)
  for i in 1:nresidues(n)
    @inbounds ni = T(n.marginals[i, marginal])
    if ni != 0.0
      H += ni * log(ni/total)
    end
  end
  (-H/total)/log(measure.base)
end

# Kullback-Leibler
# ================

"""
Kullback-Leibler (KL). This `SymmetricMeasure` has two fields.
The first is the base of the logarithm and the second is the backgroud frequency.
"""
immutable KullbackLeibler{T} <: SymmetricMeasure{T}
  base::T
  background::AbstractArray{T}
end

call{T}(::Type{KullbackLeibler{T}}, background::AbstractArray{T}) = KullbackLeibler(T(Base.e), background)

"""
`estimate(KullbackLeibler(base, background), p)`

`p` should be a `ResidueProbability` table, and `background` must have the size of `p`.
The result type is determined by `base` and `background`.
"""
function estimate{B, T, N, UseGap}(measure::KullbackLeibler{B}, p::ResidueProbability{T, N, UseGap})
  q = measure.background
  if size(q) != size(p) || length(q) != length(p)
    throw(ErrorException("p and background should have the same size and length."))
  end
  KL = zero(B)
  @inbounds for i in 1:length(p)
    pi = B(p[i])
    if pi != 0.0
      KL += pi * log(pi/q[i])
    end
  end
  KL/log(measure.base)
end

# Mutual Information
# ==================

"""
Mutual Information (MI)
"""
immutable MutualInformation{T} <: SymmetricMeasure{T}
  base::T
end

call{T}(::Type{MutualInformation{T}}) = MutualInformation(T(Base.e))

@inline _mi{T}(::Type{T}, pij, pi, pj) = ifelse(pij > zero(T) && pi > zero(T), T(pij * log(pij/(pi*pj))), zero(T))

"""
`estimate(MutualInformation(), pxy::ResidueProbability [, base])`

Calculate Mutual Information from `ResidueProbability`. The result type is determined by `base`.
"""
function estimate{B, T, UseGap}(measure::MutualInformation{B}, pxy::ResidueProbability{T, 2,UseGap})
  MI = zero(B)
  marginals = pxy.marginals
  @inbounds for j in 1:nresidues(pxy)
    pj = marginals[j,2]
    if pj > 0.0
      @inbounds @simd for i in 1:nresidues(pxy)
        MI +=  _mi(B, pxy[i,j], marginals[i,1], pj)
      end
    end
  end
  MI/log(measure.base)
end

@inline _mi{T}(N::T, nij, ni, nj) = ifelse(nij > zero(T) && ni > zero(T), T(nij * log((N * nij)/(ni * nj))), zero(T))

"""
`estimate(MutualInformation(), pxy::ResidueCount [, base])`

Calculate Mutual Information from `ResidueCount`. The result type is determined by the `base`.
It's the fastest option (you don't spend time on probability calculations).
"""
function estimate{B, T, UseGap}(measure::MutualInformation{B}, nxy::ResidueCount{T, 2,UseGap})
  MI = zero(B)
  N = B(nxy.total)
  marginals = nxy.marginals
  @inbounds for j in 1:nresidues(nxy)
    nj = marginals[j,2]
    if nj > 0.0
      @inbounds @simd for i in 1:nresidues(nxy)
        MI += _mi(N, nxy[i,j], marginals[i,1], nj)
      end
    end
  end
  (MI/N)/log(measure.base)
end

function estimate{B, T, UseGap}(measure::MutualInformation{B}, pxyz::ResidueContingencyTables{T, 3, UseGap})
  pxy = delete_dimensions(pxyz,3)
  return( estimate_on_marginal(Entropy(measure.base), pxyz, 1) + # H(X)
  estimate_on_marginal(Entropy(measure.base), pxyz ,2) + # H(Y)
  estimate_on_marginal(Entropy(measure.base), pxyz, 3) - # H(Z)
  estimate(Entropy(measure.base), pxy) - # H(X, Y)
  estimate(Entropy(measure.base), delete_dimensions!(pxy, pxyz, 2)) - # H(X, Z)
  estimate(Entropy(measure.base), delete_dimensions!(pxy, pxyz, 1)) + # H(Y, Z)
  estimate(Entropy(measure.base), pxyz) ) # H(X, Y, Z)
end

# Normalized Mutual Information by Entropy
# ========================================

"""
Normalized Mutual Information (nMI) by Entropy.

`nMI(X, Y) = MI(X, Y) / H(X, Y)`
"""
immutable MutualInformationOverEntropy{T} <: SymmetricMeasure{T}
  base::T
end

call{T}(::Type{MutualInformationOverEntropy{T}}) = MutualInformationOverEntropy(T(Base.e))

function estimate{B}(measure::MutualInformationOverEntropy{B}, table)
  H = estimate(Entropy(measure.base), table)
  if H != zero(B)
    MI = estimate(MutualInformation(measure.base), table)
    return(MI/H)
  else
    return(zero(B))
  end
end

# Pairwise Gap Percentage
# =======================

"""
`GapUnionPercentage`
"""
immutable GapUnionPercentage{T} <: SymmetricMeasure{T} end

"""
`GapIntersectionPercentage`
"""
immutable GapIntersectionPercentage{T} <: SymmetricMeasure{T} end

function estimate{B, T}(measure::GapIntersectionPercentage{B}, nxy::ResidueCount{T, 2, true})
  B(100.0) * B(nxy[21, 21]) / B(sum(nxy))
end

function estimate{B, T}(measure::GapUnionPercentage{B}, nxy::ResidueCount{T, 2, true})
  B(100.0) * B(nxy.marginals[21, 1] + nxy.marginals[21, 2] - nxy[21, 21]) / B(sum(nxy))
end


text:
using Test
using TestSetExtensions
using DataLoaders
using DataLoaders:
    BatchViewCollated, BatchDimLast, BatchDimFirst, WorkerPool, RingBuffer, collate, obsslices, _batchsize
using MLDataPattern
using LearnBase
using Random

struct MockDataset
    n::Int
    sz::Any
    inplace::Bool
end

LearnBase.getobs(ds::MockDataset, idx::Int) = (randn(ds.sz...))
LearnBase.getobs(ds::MockDataset, idxs) = [getobs(ds, idx) for idx in idxs]
LearnBase.getobs!(buf, ds::MockDataset, idx::Int) =
    ds.inplace ? fill!(buf, 0.3) : getobs(ds, idx)
LearnBase.nobs(ds::MockDataset) = ds.n

@testset ExtendedTestSet "collate" begin
    @test collate([([1, 2], 3), ([4, 5], 6)]) == ([1 4; 2 5], [3, 6])
    @test collate([(x = [1, 2], y = 3), (x = [4, 5], y = 6)]) ==
          (x = [1 4; 2 5], y = [3, 6])
    @test collate([Dict("x" => [1, 2], "y" => 3), Dict("x" => [4, 5], "y" => 6)]) ==
          Dict("x" => [1 4; 2 5], "y" => [3, 6])
    @test collate([(1, 2), (3, 4)]) == ([1, 3], [2, 4])
end

@testset ExtendedTestSet "obsslices" begin
    batch = rand(10, 10, 16)
    @test_nowarn for obs in obsslices(batch, BatchDimLast())
        @assert size(obs) == (10, 10)
    end

    @test_nowarn for obs in obsslices(batch, BatchDimFirst())
        @assert size(obs) == (10, 16)
    end

    batch2 = Dict(:x => rand(5, 16), :y => rand(1, 16))
    obs = first(obsslices(batch2))
    @test obs isa Dict
    @test size(obs[:x]) == (5,)
    @test size(obs[:y]) == (1,)

    batch3 = (rand(5,16), rand(3,16))
    obs = first(obsslices(batch3))
    @test obs isa Tuple
    @test size(obs[1]) == (5,)
    @test size(obs[2]) == (3,)

    batch4 = (a=rand(5,16), b=rand(3,16))
    obs = first(obsslices(batch4))
    @test obs isa NamedTuple
    @test size(obs.a) == (5,)
    @test size(obs.b) == (3,)
end

@testset ExtendedTestSet "_batchsize" begin
    batch = (rand(10,3,16), rand(10,4,16))
    @test batch isa Tuple
    @test _batchsize(batch, BatchDimFirst()) == 10
    @test _batchsize(batch, BatchDimLast())  == 16

    batch = (a=rand(10,3,16), b=rand(10,4,16))
    @test batch isa NamedTuple
    @test _batchsize(batch, BatchDimFirst()) == 10
    @test _batchsize(batch, BatchDimLast())  == 16
end

@testset ExtendedTestSet "BatchViewCollated" begin
    data = rand(2, 100)
    @testset ExtendedTestSet "basic" begin
        bv = DataLoaders.BatchViewCollated(data, 5)
        @test nobs(bv) == 20
        @test_nowarn for i = 1:nobs(bv)
            getobs(bv, i)
        end
    end

    @testset ExtendedTestSet "getobs!" begin
        bv = DataLoaders.BatchViewCollated(data, 5)
        buf = getobs(bv, 1)
        getobs!(buf, bv, 2)
        @test buf == getobs(bv, 2)
    end

    @testset ExtendedTestSet "droplast" begin
        bv = batchviewcollated(data, 3, partial = true)
        @test nobs(bv) == 34

    end

    @testset ExtendedTestSet "inplace loading for non-supporting datasets" begin
        ds = MockDataset(64, (10,), false)
        bv = BatchViewCollated(ds, 4)
        buf = getobs(bv, 1)
        x = buf[1]
        getobs!(buf, bv, 2)
        @test x != buf[1]
    end
end

@testset ExtendedTestSet "RingBuffer" begin
    @testset ExtendedTestSet "fills correctly" begin
        x = rand(10, 10)
        ringbuffer = RingBuffer(10, x)
        @async begin
            for i ∈ 1:100
                put!(ringbuffer) do buf
                    rand!(buf)
                end
            end
        end
        @test_nowarn for _ ∈ 1:100
            take!(ringbuffer)
        end
    end

    @testset ExtendedTestSet "does mutate" begin
        x = rand(10, 10)
        ringbuffer = RingBuffer(1, copy(x))
        put!(ringbuffer) do buf
            @test x ≈ buf
            copy!(buf, rand(10, 10))
            buf
        end
        x_ = take!(ringbuffer)
        @test !(x ≈ x_)
    end
end

@testset ExtendedTestSet "WorkerPool" begin
    if Threads.nthreads() == 1
        @test_broken false
    else
        xs = falses(10)

        @test !any(xs)
        pool = WorkerPool(collect(zip(1:10))) do i
            xs[i] = true
        end
        DataLoaders.run(pool)
        @test all(xs)
    end

end

@testset ExtendedTestSet "eachobsparallel unbuffered" begin
    make() = eachobsparallel(rand(10, 64), buffered = false, useprimary = true)

    @testset ExtendedTestSet "iterate" begin
        dl = make()
        x, (results, workerpool, idx) = iterate(dl)
        @test idx == 1
        @test_nowarn for obs in dl
        end
    end
end

@testset ExtendedTestSet "eachobsparallel buffered" begin
    make() = eachobsparallel(rand(10, 64), buffered = true, useprimary = true)

    @testset ExtendedTestSet "iterate" begin
        dl = make()
        x, (ringbuffer, workerpool, idx) = iterate(dl)
        @test idx == 1
        @test_nowarn for obs in dl
        end
    end
end

@testset ExtendedTestSet "batchviewcollated" begin
    @testset ExtendedTestSet "basic" begin
        data = rand(10, 64)
        bv = DataLoaders.batchviewcollated(data, 4)
        @test nobs(bv) == 16
        @test_nowarn for _ in obsview(bv) end
    end

    @testset ExtendedTestSet "partial" begin
        data = rand(10, 10)
        bv = DataLoaders.batchviewcollated(data, 3)
        @test nobs(bv) == 4
        buf = getobs(bv, 4)
        buf = getobs!(buf, bv, 4)
    end

end

@testset ExtendedTestSet "DataLoader" begin
    data = MockDataset(256, (10, 5), true)
    bs = 8

    @testset ExtendedTestSet "buffer, collate, parallel" begin
        dl = DataLoader(data, bs)
        @test_nowarn for batch in dl end
    end

    @testset ExtendedTestSet "buffer, collate, parallel, samples" begin
        dl = DataLoader(data, nothing)
        @test_nowarn for batch in dl end
    end

    @testset ExtendedTestSet "collate, parallel" begin
        dl = DataLoader(data, bs, buffered = false)
        @test_nowarn for batch in dl end
    end

    @testset ExtendedTestSet "collate" begin
        dl = DataLoader(data, bs, buffered = false)
        @test_nowarn for batch in dl end
    end

    @testset ExtendedTestSet "buffer, collate" begin
        dl = DataLoader(data, bs, buffered = true)
        @test_nowarn for batch in dl end
    end
end


text:
"""
```
AbstractDSGEModel{T} <: AbstractModel{T}
```

The AbstractDSGEModel is defined as a subtype of AbstractModel
"""
abstract type AbstractDSGEModel{T} <: ModelConstructors.AbstractModel{T} end

function Base.show(io::IO, m::AbstractDSGEModel)
    @printf io "Dynamic Stochastic General Equilibrium Model\n"
    @printf io "no. states:             %i\n" n_states(m)
    @printf io "no. anticipated policy shocks: %i\n" n_mon_anticipated_shocks(m)
    @printf io "data vintage:           %s\n" data_vintage(m)
    @printf io "description:\n %s\n"          description(m)
end

# Retrieve model type
model_type(m::AbstractDSGEModel{T}) where T = T

# Number of anticipated policy shocks
n_anticipated_shocks(m::AbstractDSGEModel) = get_setting(m, :n_mon_anticipated_shocks)
n_anticipated_shocks_padding(m::AbstractDSGEModel) = get_setting(m, :n_mon_anticipated_shocks_padding)
n_mon_anticipated_shocks(m::AbstractDSGEModel) = get_setting(m, :n_mon_anticipated_shocks)
n_mon_anticipated_shocks_padding(m::AbstractDSGEModel) = get_setting(m, :n_mon_anticipated_shocks_padding)
n_anticipated_shocks_padding(m::AbstractDSGEModel) = get_setting(m, :n_mon_anticipated_shocks_padding)
n_z_anticipated_shocks(m::AbstractDSGEModel) = get_setting(m, :n_z_anticipated_shocks)
n_z_anticipated_shocks_padding(m::AbstractDSGEModel) = get_setting(m, :n_z_anticipated_shocks_padding)

# Dates, indices, number of periods for each regime
date_presample_start(m::AbstractDSGEModel) = get_setting(m, :date_presample_start)
date_mainsample_start(m::AbstractDSGEModel) = get_setting(m, :date_mainsample_start)
date_zlb_start(m::AbstractDSGEModel) = get_setting(m, :date_zlb_start)
date_zlb_end(m::AbstractDSGEModel) = get_setting(m, :date_zlb_end)
date_presample_end(m::AbstractDSGEModel) = Dates.lastdayofquarter(get_setting(m, :date_mainsample_start) - Dates.Month(3))
date_prezlb_end(m::AbstractDSGEModel) = Dates.lastdayofquarter(get_setting(m, :date_zlb_start) - Dates.Month(3))
date_mainsample_end(m::AbstractDSGEModel) = Dates.lastdayofquarter(get_setting(m, :date_forecast_start) - Dates.Month(3))
date_conditional_end(m::AbstractDSGEModel) = get_setting(m, :date_conditional_end)

index_presample_start(m::AbstractDSGEModel) = 1
index_mainsample_start(m::AbstractDSGEModel) = subtract_quarters(date_mainsample_start(m), date_presample_start(m)) + 1
index_zlb_start(m::AbstractDSGEModel) = subtract_quarters(date_zlb_start(m), date_presample_start(m)) + 1
index_forecast_start(m::AbstractDSGEModel) = subtract_quarters(date_forecast_start(m), date_presample_start(m)) + 1

"""
```
index_shockdec_start(m::AbstractDSGEModel)
```

Returns the index starting from which the shock decomposition is saved, where 1 is the index corresponding to date_mainsample_start(m).
"""
index_shockdec_start(m::AbstractDSGEModel) = subtract_quarters(date_shockdec_start(m), date_mainsample_start(m)) + 1

"""
```
index_shockdec_end(m::AbstractDSGEModel)
```

Returns the last index for which the shock decomposition is saved, where 1 is the index corresponding to date_mainsample_start(m).
"""
index_shockdec_end(m::AbstractDSGEModel) = subtract_quarters(date_shockdec_end(m), date_mainsample_start(m)) + 1

n_presample_periods(m::AbstractDSGEModel)   = subtract_quarters(date_mainsample_start(m), date_presample_start(m))
n_prezlb_periods(m::AbstractDSGEModel)      = subtract_quarters(date_zlb_start(m), date_mainsample_start(m))
n_zlb_periods(m::AbstractDSGEModel)         = subtract_quarters(date_forecast_start(m), date_zlb_start(m))
n_mainsample_periods(m::AbstractDSGEModel)  = subtract_quarters(date_forecast_start(m), date_mainsample_start(m))
n_conditional_periods(m::AbstractDSGEModel) = subtract_quarters(date_conditional_end(m), date_mainsample_end(m))

inds_presample_periods(m::AbstractDSGEModel)  = collect(index_presample_start(m):(index_mainsample_start(m)-1))
inds_prezlb_periods(m::AbstractDSGEModel)     = collect(index_mainsample_start(m):(index_zlb_start(m)-1))
inds_zlb_periods(m::AbstractDSGEModel)        = collect(index_zlb_start(m):(index_forecast_start(m)-1))
inds_mainsample_periods(m::AbstractDSGEModel) = collect(index_mainsample_start(m):(index_forecast_start(m)-1))

# Convenience functions for working with heterogeneous agent models
# that differentiate between backward looking "state" variables and "jump" variables
n_backward_looking_states(m::AbstractDSGEModel) = get_setting(m, :n_backward_looking_states)
n_jumps(m::AbstractDSGEModel) = get_setting(m, :n_jumps)
n_model_states(m::AbstractDSGEModel) = get_setting(m, :n_model_states)

# The numbers for n_states, and n_jumps assumes normalization
# There is a procedure in klein_solve that normalizes the state variable grids
# by removing an entry from them. Hence the number of states and jumps being
# tracked should be 1 less if normalized
# However, the number of jumps and states unnormalized are required for
# the construction of the Jacobian, hence the reason for these helpers
n_jumps_unnormalized(m::AbstractDSGEModel) = n_jumps(m) + get_setting(m, :jumps_normalization_factor)
n_backward_looking_states_unnormalized(m::AbstractDSGEModel) = n_backward_looking_states(m) + get_setting(m, :backward_looking_states_normalization_factor)
n_model_states_unnormalized(m::AbstractDSGEModel) = n_jumps_unnormalized(m) + n_backward_looking_states_unnormalized(m)
n_model_states_original(m::AbstractDSGEModel) = get_setting(m, :n_model_states_original)

"""
```
AbstractCTModel{T} <: AbstractDSGEModel{T}
```

The AbstractCTModel is defined as a subtype of AbstractDSGEModel to accomodate the
numerical methods and procedures specific to continuous time models.
"""
abstract type AbstractCTModel{T} <: AbstractDSGEModel{T} end

n_states(m::AbstractCTModel) = sum(map(i -> length(collect(m.endogenous_states)[i][2]), 1:length(keys(m.endogenous_states))))
n_shocks_expectational(m::AbstractCTModel) = sum(map(i -> length(collect(m.expected_shocks)[i][2]), 1:length(keys(m.expected_shocks))))

"""
```
AbstractHetModel{T} <: AbstractDSGEModel{T}
```

The AbstractHetModel is defined as a subtype of AbstractDSGEModel to accomodate a bunch of stuff, but for now just different impulse response functions.
"""
abstract type AbstractHetModel{T} <: AbstractDSGEModel{T} end

"""
```
AbstractRepModel{T} <: AbstractDSGEModel{T}
```

The AbstractRepresentativeModel is defined as a subtype of AbstractDSGEModel to accomodate a bunch of stuff, but for now just different impulse response functions.
"""
abstract type AbstractRepModel{T} <: AbstractDSGEModel{T} end

"""
```
get_dict(m, class, index)
```
"""
function get_dict(m::AbstractDSGEModel, class::Symbol)
    if class == :states
        m.endogenous_states
    elseif class == :obs
        m.observables
    elseif class == :pseudo
        m.pseudo_observables
    elseif class in [:shocks, :stdshocks]
        m.exogenous_shocks
    else
        throw(ArgumentError("Invalid class: $class. Must be :states, :obs, :pseudo, :shocks, or :stdshocks"))
    end
end

"""
```
get_key(m, class, index)
```
Returns the name of the state (`class = :states`), observable (`:obs`),
pseudo-observable (`:pseudo`), or shock (`:shocks` or `:stdshocks`)
corresponding to the given `index`.
"""
function get_key(m::AbstractDSGEModel, class::Symbol, index::Int)
    dict = get_dict(m, class)
    out = Base.filter(key -> dict[key] == index, collect(keys(dict)))
    if length(out) == 0
        error("Key corresponding to index $index not found for class: $class")
    elseif length(out) > 1
        error("Multiple keys corresponding to index $index found for class: $class")
    else
        return out[1]
    end
end

# Parse population mnemonic into 2 Nullable{Symbol}s from one
function parse_population_mnemonic(m::AbstractDSGEModel)
    mnemonic = get_setting(m, :population_mnemonic)
    if isnull(mnemonic)
        return [Nullable{Symbol}(), Nullable{Symbol}()]
    else
        return map(s -> Nullable(Symbol(s)), split(string(get(mnemonic)), DSGE_DATASERIES_DELIM))
    end
end

# From an augmented state space with anticipated policy shocks, get indices
# corresponding to pre-ZLB states, shocks, and observables
function inds_states_no_ant(m::AbstractDSGEModel)
    if n_mon_anticipated_shocks(m) > 0
        ind_ant1 = m.endogenous_states[:rm_tl1]
        ind_antn = m.endogenous_states[Symbol("rm_tl$(n_mon_anticipated_shocks(m))")]
        return [1:(ind_ant1-1); (ind_antn+1):n_states_augmented(m)]
    else
        return collect(1:n_states_augmented(m))
    end
end

function inds_shocks_no_ant(m::AbstractDSGEModel)
    if n_mon_anticipated_shocks(m) > 0
        ind_ant1 = m.exogenous_shocks[:rm_shl1]
        ind_antn = m.exogenous_shocks[Symbol("rm_shl$(n_mon_anticipated_shocks(m))")]
        return [1:(ind_ant1-1); (ind_antn+1):n_shocks_exogenous(m)]
    else
        return collect(1:n_shocks_exogenous(m))
    end
end

function inds_obs_no_ant(m::AbstractDSGEModel)
    if n_mon_anticipated_shocks(m) > 0
        ind_ant1 = m.observables[:obs_nominalrate1]
        ind_antn = m.observables[Symbol("obs_nominalrate$(n_mon_anticipated_shocks(m))")]
        return [1:(ind_ant1-1); (ind_antn+1):n_observables(m)]
    else
        return collect(1:n_observables(m))
    end
end

# From an augmented state space with integrated series (e.g. unit root),
# get indices corresponding to stationary states
function inds_states_no_integ_series(m::AbstractDSGEModel)
    if haskey(get_settings(m), :integrated_series)
        inds = map(i -> m.endogenous_states_augmented[i],
                   get_setting(m, :integrated_series))
        return setdiff(1:n_states_augmented(m), inds)
    else
        return collect(1:n_states_augmented(m))
    end
end

# Interface for accessing parameters
get_parameters(m::AbstractDSGEModel) = m.parameters

# Interface for accessing rng
get_rng(m::AbstractDSGEModel) = m.rng

# Interface for accessing settings dictionary
get_settings(m::AbstractDSGEModel) = hasproperty(m, :testing) ? (m.testing ? m.test_settings : m.settings) : m.settings

# Interface for accessing observables dictionary
get_observables(m::AbstractDSGEModel) = m.observables
get_pseudo_observables(m::AbstractDSGEModel) = m.pseudo_observables
get_exogenous_shocks(m::AbstractDSGEModel) = m.exogenous_shocks

# Interface for data
cond_vintage(m::AbstractDSGEModel)    = get_setting(m, :cond_vintage)
cond_id(m::AbstractDSGEModel)         = get_setting(m, :cond_id)
cond_full_names(m::AbstractDSGEModel) = get_setting(m, :cond_full_names)
cond_semi_names(m::AbstractDSGEModel) = get_setting(m, :cond_semi_names)
use_population_forecast(m::AbstractDSGEModel) = get_setting(m, :use_population_forecast)
hpfilter_population(m::AbstractDSGEModel)     = get_setting(m, :hpfilter_population)

# Interface for general computation settings
use_parallel_workers(m::AbstractDSGEModel)    = get_setting(m, :use_parallel_workers)

# Interface for estimation settings
reoptimize(m::AbstractDSGEModel)          = get_setting(m, :reoptimize)
calculate_hessian(m::AbstractDSGEModel) = get_setting(m, :calculate_hessian)
hessian_path(m::AbstractDSGEModel)      = get_setting(m, :hessian_path)
n_hessian_test_params(m::AbstractDSGEModel) = get_setting(m, :n_hessian_test_params)

# Interface for Metropolis-Hastings settings
n_mh_blocks(m::AbstractDSGEModel)       =  get_setting(m, :n_mh_blocks)
n_mh_param_blocks(m::AbstractDSGEModel) =  get_setting(m, :n_mh_param_blocks)
n_mh_simulations(m::AbstractDSGEModel)  =  get_setting(m, :n_mh_simulations)
n_mh_burn(m::AbstractDSGEModel)         =  get_setting(m, :n_mh_burn)
mh_thin(m::AbstractDSGEModel)           =  get_setting(m, :mh_thin)

# Interface for forecast settings
date_forecast_start(m::AbstractDSGEModel)   = get_setting(m, :date_forecast_start)
forecast_block_size(m::AbstractDSGEModel)   = get_setting(m, :forecast_block_size)
forecast_start_block(m::AbstractDSGEModel)  = get_setting(m, :forecast_start_block)
forecast_input_file_overrides(m::AbstractDSGEModel) = get_setting(m, :forecast_input_file_overrides)
forecast_uncertainty_override(m::AbstractDSGEModel) = get_setting(m, :forecast_uncertainty_override)
forecast_smoother(m::AbstractDSGEModel)     = get_setting(m, :forecast_smoother)
forecast_tdist_df_val(m::AbstractDSGEModel) = get_setting(m, :forecast_tdist_df_val)
forecast_tdist_shocks(m::AbstractDSGEModel) = get_setting(m, :forecast_tdist_shocks)
forecast_zlb_value(m::AbstractDSGEModel)    = get_setting(m, :forecast_zlb_value)
impulse_response_horizons(m::AbstractDSGEModel) = get_setting(m, :impulse_response_horizons)
n_shockdec_periods(m::AbstractDSGEModel)    = index_shockdec_end(m) - index_shockdec_start(m) + 1

# Function for getting the alternative policy. Note that, in the case of imperfect awareness,
# i.e. when uncertain_altpolicy is true, this function will return the policy that is actually implemented,
# instead of other potential alternative policies which households believe could occur.
#
# Furthermore, do NOT delete the backup cases when regime_eqcond_info is not defined, namely
# the check for whether :alternative_policy is a Setting. This last check is necessary for
# the scenarios code to continue working.
#
# Finally, the check for :regime_eqcond_info MUST come first, or else the regime-switching code
# will not work as intended, namely that :regime_eqcond_info specifies equilibrium conditions
# which differ from the default ones and thus should determine the "final" permanent policy
# that is believed to hold until "the end of time" (i.e. forever since it is an infinite horizon model).
alternative_policy(m::AbstractDSGEModel) = haskey(get_settings(m), :regime_eqcond_info) && # no check for n_regimes b/c if it is not there, then
    haskey(get_setting(m, :regime_eqcond_info), get_setting(m, :n_regimes)) ?              # it is better to throw an error since regime_eqcond_info
    get_setting(m, :regime_eqcond_info)[get_setting(m, :n_regimes)].alternative_policy :   # will not work in general anyway.
    (haskey(get_settings(m), :alternative_policy) ? get_setting(m, :alternative_policy) : AltPolicy(:historical, eqcond, solve))

# Some additional date settings related to forecasts
function date_forecast_end(m::AbstractDSGEModel)
    if haskey(get_settings(m), :date_forecast_end)
        return get_setting(m, :date_forecast_end)
    else
        date = date_forecast_start(m) + Dates.Month(3 * (forecast_horizons(m)-1))
        return Dates.lastdayofquarter(date)
    end
end

function forecast_horizons(m::AbstractDSGEModel; cond_type::Symbol = :none)
    horizons = get_setting(m, :forecast_horizons)
    if cond_type == :none
        return horizons
    else
        return horizons - n_conditional_periods(m)
    end
end

function date_shockdec_start(m::AbstractDSGEModel)
    startdate = get_setting(m, :shockdec_startdate)
    if (Nullables.isnull(startdate)) | (startdate===nothing)
        return DSGE.date_mainsample_start(m)
    else
        return get(startdate)
    end
end

function date_shockdec_end(m::AbstractDSGEModel)
    enddate =  get_setting(m, :shockdec_enddate)
    if !isnull(enddate)
        return get(enddate)
    else
        return date_forecast_end(m)
    end
end

"""
```
load_parameters_from_file(m::AbstractDSGEModel,path::String)
```
Returns a vector of parameters, read from a file, suitable for updating `m`.
"""
function load_parameters_from_file(m::AbstractDSGEModel, path::String)

    if isfile(path) && splitext(path)[2] == ".h5"
        x  = h5open(path, "r") do file
            try
                read(file, "params")
            catch
                error("$path does not contain variable params")
            end
        end
    else
        error("$path is not a valid HDF5 file.")
    end

    if m.spec=="smets_wouters" && length(x)==42
        x = vcat(x, zeros(7))
    end

    @assert length(x) == length(m.parameters)
    @assert eltype(x) == typeof(m.parameters[1].value)
    return x
end

"""
```
specify_mode!(m::AbstractDSGEModel, mode_file::String=""; verbose=:low)
```

Updates the values of `m.parameters` with the values from
`mode_file`. Sets `reoptimize` setting to `false`.

Usage: should be run before calling `estimate(m)`, e.g.:

    m = Model990()
    specify_mode!(m, modefile)
    estimate(m)
"""
function specify_mode!(m::AbstractDSGEModel, mode_file::String = ""; verbose=:low)

    m <= Setting(:reoptimize, false)

    if mode_file == ""
        mode_file = inpath(m, "user", "paramsmode.h5")
    end
    mode_file = normpath(mode_file)

    DSGE.update!(m,load_parameters_from_file(m,mode_file))

    println(verbose, :low, "Loaded previous mode from $mode_file.")
end

"""
```
specify_hessian!(m::AbstractDSGEModel, path::String=""; verbose=:low)
```

Specify a Hessian matrix calculated at the posterior mode to use in the model estimation. If
no path is provided, will attempt to detect location.
"""
function specify_hessian!(m::AbstractDSGEModel, path::String=""; verbose=:low)
    if isempty(path)
        path = inpath(m, "user", "hessian.h5")
    end

    if isfile(path) && splitext(path)[2] == ".h5"
        m <= Setting(:hessian_path, normpath(abspath(path)))
    else
        error("Invalid input Hessian file: $path",)
    end

    m <= Setting(:calculate_hessian, false)

    println(verbose, :low, "Specified hessian from $path.")
end

"""
```
transform_to_model_space!(m::AbstractDSGEModel, values::Vector{T}; regime_switching::Bool = false) where T<:AbstractFloat
```

Transforms `values` from the real line to the model space, and assigns `values[i]` to
`m.parameters[i].value` for non-steady-state parameters. Recomputes the steady-state
parameter values.

### Arguments
- `m`: the model object
- `values`: the new values to assign to non-steady-state parameters.
- `regime_switching`: set to true if the model's parameters are regime-switching
"""
function transform_to_model_space!(m::AbstractDSGEModel, values::Vector{T};
                                   regime_switching::Bool = false) where {T<:AbstractFloat}
    new_values = transform_to_model_space(m.parameters, values; regime_switching = regime_switching)
    DSGE.update!(m, new_values)
    steadystate!(m)
end

"""
```
update!(m::AbstractDSGEModel, values::Vector{T}) where T<:AbstractFloat
```

Update `m.parameters` with `values`, recomputing the steady-state parameter values.

### Arguments:
- `m`: the model object
- `values`: the new values to assign to non-steady-state parameters.
"""
function update!(m::AbstractDSGEModel, values::AbstractVector{T}) where {T <: Real}
    ModelConstructors.update!(m.parameters, values)
    steadystate!(m)
end

"""
```
update!(m::AbstractDSGEModel, values::ParameterVector{T};
    regime_switching::Bool = false, toggle::Bool = true) where T
```
Update `m.parameters` with `values`, recomputing the steady-state parameter values.

### Arguments
- `m`: the model object
- `values`: the new values to assign to non-steady-state parameters.

### Keyword
- `regime_switching`: if true, then we assume the parameters are regime-switching,
    in which case `update!` assumes the `value` field
    of each parameter in values` holds the parameter value in the first regime, and
    then we update the field `regimes` for each parameter
- `toggle`: if true, we call `ModelConstructors.toggle_regime!(values)` before
    updating any values to ensure the `value` field of the parameters in `values`
    correspond to regime 1 values.
"""
function update!(m::AbstractDSGEModel, values::ParameterVector{T};
                 regime_switching::Bool = false, toggle::Bool = true) where {T <: Real}

    # Update regime-switching if length of `values` exceeds m.parameters
    if regime_switching
        if toggle
            ModelConstructors.toggle_regime!(values, 1)
        end

        # Update first-regime values
        ModelConstructors.update!(m.parameters, [θ.value for θ in values])

        # Update remaining regimes
        for (i, para) in enumerate(m.parameters)
            if !isempty(para.regimes)
                for (ind, val) in para.regimes[:value]
                    if ind == 1
                        ModelConstructors.set_regime_val!(para, 1, para.value)
                    else
                        ModelConstructors.set_regime_val!(para, ind, regime_val(values[i], ind))
                    end
                end
            end
        end
    else
        ModelConstructors.update!(m.parameters, [θ.value for θ in values])
    end

    steadystate!(m)
end

"""
```
mutable struct ShockGroup
```

The `ShockGroup` mutable struct is used in `prepare_means_table_shockdec` and
`plot_shock_decompositions`. When plotting shock decompositions, we usually want
to group the shocks into categories (financial, monetary policy, etc.) so that
the resulting grouped bar plot is legible.

### Fields

- `name::String`
- `shocks::Vector{Symbol}`
- `color::Colorant`
"""
mutable struct ShockGroup
    name::String
    shocks::Vector{Symbol}
    color::Colorant
end

function ShockGroup(name::String, shocks::Vector{Symbol}, color_name::Symbol)
    color = parse(Colorant, color_name)
    return ShockGroup(name, shocks, color)
end

mutable struct SteadyStateConvergenceError <: Exception
    msg::String
end
SteadyStateConvergenceError() = SteadyStateConvergenceError("SteadyState didn't converge")
Base.showerror(io::IO, ex::SteadyStateConvergenceError) = print(io, ex.msg)

"""
```
setup_regime_switching_inds!(m::AbstractDSGEModel; cond_type::Symbol = :none)
```

calculates the indices needed to solve and forecast a model with regime-switching.

### Keywords
- `cond_type`: the correct regime indices for forecasting depend on whether the forecast is conditional or not
"""
function setup_regime_switching_inds!(m::AbstractDSGEModel; cond_type::Symbol = :none)

    if haskey(get_settings(m), :regime_switching) ? !get_setting(m, :regime_switching) : true
        @warn "The setting :regime_switching is either false or is not defined yet. Updating the setting to be true."
        m <= Setting(:regime_switching, true)
    end

    n_hist_regimes = 0
    n_cond_regimes = 0
    n_regimes = length(get_setting(m, :regime_dates))
    post_cond_end = iterate_quarters(date_conditional_end(m), 1) # Period after conditional forecasting ends
    set_reg_forecast_start   = false # These flags are needed to tell whether or not we actually set these regimes.
    set_post_conditional_end = false # W/out these flags, if m already has reg_forecast_start, then it won't be properly set
    reg_forecast_start_str = "Regime in which the forecast starts."
    reg_post_conditional_end_str = "Regime one period after the conditional forecast ends."
    for (key, val) in get_setting(m, :regime_dates)
        if val == date_forecast_start(m)
            m <= Setting(:reg_forecast_start, key, reg_forecast_start_str)
            set_reg_forecast_start = true
        end
        if val == post_cond_end
            m <= Setting(:reg_post_conditional_end, key, reg_post_conditional_end_str)
            set_post_conditional_end = true
        end
        if val < date_forecast_start(m)
            n_hist_regimes += 1
        end
        if date_forecast_start(m) <= val <= date_conditional_end(m)
            n_cond_regimes += 1
        end
    end
    if !set_reg_forecast_start
        # Then the forecast begins in the middle of a regime
        m <= Setting(:reg_forecast_start, findlast(sort!(collect(values(get_setting(m, :regime_dates)))) .< date_forecast_start(m)),
                     reg_forecast_start_str)
    end
    if !set_post_conditional_end
        # Then the conditional forecast ends in the middle of a regime
        m <= Setting(:reg_post_conditional_end,
                     findlast(sort!(collect(values(get_setting(m, :regime_dates)))) .<= date_conditional_end(m)),
                     reg_post_conditional_end_str) # or .< post_cond_end
    end

    # Infer number of regimes in the forecast horizon from reg_forecast_start
    n_fcast_regimes = n_regimes - get_setting(m, :reg_forecast_start) + 1

    m <= Setting(:n_regimes, n_regimes, "Total number of regimes")
    m <= Setting(:n_hist_regimes, n_hist_regimes, "Number of regimes in the history")
    m <= Setting(:n_fcast_regimes, n_fcast_regimes, "Number of regimes in the forecast horizon")
    m <= Setting(:n_cond_regimes, n_cond_regimes, "Number of regime switches during the conditional forecast horizon")

    return m
end

# Dummy function for eqcond, relevant when using gensys2
function eqcond(m::AbstractDSGEModel, Γ0::AbstractMatrix{S}, Γ1::AbstractMatrix{S},
                C::AbstractVector{S}, Ψ::AbstractMatrix{S}, Π::AbstractMatrix{S}) where {S <: Real}
    return Γ0, Γ1, C, Ψ, Π
end

"""
```
setup_param_regimes!(m::AbstractDSGEModel, param_regs::Matrix{Int} = []

Function to set up the model with parameter regime switching
for estimation. param_regs should be a matrix where for each row r,
the i{th} element is the parameter regime for the i{th} model regime
for parameter r (including parameters with no regime-switching: set
the rows for non-regime-switching parameters to all 1s, although any
value will give the same result).
```
"""
function setup_param_regimes!(m::AbstractDSGEModel, param_mat::Array{Int, 2} = Matrix{Int}(undef, 0, 0))

    param_reg   = Dict{Symbol, Dict{Int, Int}}()
    nmodel_regs = haskey(get_settings(m), :n_regimes) ? get_setting(m, :n_regimes) : 1

    if isempty(param_mat)
        param_mat = ones(n_parameters(m), nmodel_regs)
    end

    @assert nmodel_regs == size(param_mat, 2) "The number of columns in `param_mat` must match the number of model regimes"

    for i in 1:n_parameters(m)
        reg_dict = Dict{Int, Int}()
        for reg in 1:nmodel_regs
            reg_dict[reg] = param_mat[i, reg]
        end
        param_reg[m.parameters[i].key] = reg_dict
    end

    m <= Setting(:model2para_regime, param_reg)

    return m
end


text:
function generate_stochiometric_matrix_component(ir_dictionary::Dict{String,Any})::NamedTuple

    # initialize -
    number_of_species = ir_dictionary["number_of_species"]
    system_type = ir_dictionary["system_type"]
    number_of_transcription_models = ir_dictionary["number_of_transcription_models"]
    number_of_translation_models = ir_dictionary["number_of_translation_models"]
    number_of_rates = number_of_transcription_models + number_of_translation_models

    # get the species table -
    model_species_table = ir_dictionary["model_species_table"]

    try

        # what is the overall dimension of the system -
        (number_of_species, _) = size(model_species_table)

        # how many genes do we have?
        number_of_genes = filter(x -> x == :DNA, model_species_table[:,:type]) |> length

        # blocks -
        gene_block = zeros(number_of_genes, number_of_rates)
        id_block = Matrix{Float64}(I, (number_of_species - number_of_genes), (number_of_species - number_of_genes))

        # stm -
        stoichiometric_matrix = [gene_block ; id_block]

        # ok, if we have a CF_ system type, then we have an extra row of zeros -
        if (contains(system_type,"CF_") == true)
            zero_row = zeros(1,number_of_rates)
            stoichiometric_matrix = [stoichiometric_matrix ; zero_row]
        end

        # package -
        stm_program_component = (matrix = stoichiometric_matrix, filename = "Network.dat", component_type = :matrix)

        # return -
        return stm_program_component
    catch error
        rethrow(error)
    end
end

function generate_dilution_degradation_matrix_component(ir_dictionary::Dict{String,Any})::NamedTuple

    # initialize -
    number_of_species = ir_dictionary["number_of_species"]
    system_type = ir_dictionary["system_type"]
    number_of_transcription_models = ir_dictionary["number_of_transcription_models"]
    number_of_translation_models = ir_dictionary["number_of_translation_models"]
    number_of_rates = number_of_transcription_models + number_of_translation_models

    # get the species table -
    model_species_table = ir_dictionary["model_species_table"]

    try

        # what is the overall dimension of the system -
        (number_of_species, _) = size(model_species_table)

        # how many genes do we have?
        number_of_genes = filter(x -> x == :DNA, model_species_table[:,:type]) |> length

        # blocks -
        gene_block = zeros(number_of_genes, number_of_rates)
        id_block = -1 * Matrix{Float64}(I, (number_of_species - number_of_genes), (number_of_species - number_of_genes))

        # stm -
        degradation_matrix = [gene_block ; id_block]

        # ok, if we have a CF_ system type, then we have an extra row of zeros -
        if (contains(system_type,"CF_") == true)
            zero_row = zeros(1,number_of_rates)
            degradation_matrix = [degradation_matrix ; zero_row]
        end

        # package -
        degradation_program_component = (matrix = degradation_matrix, filename = "Degradation.dat", component_type = :matrix)

        # return -
        return degradation_program_component
    catch error
        rethrow(error)
    end
end


text:
using FileIO

function load_tensor(datasetname)
    if datasetname == "AttFace"
        T = load("../data/tensor_att_face.jld2", "tensor_att_face")
    elseif datasetname == "4DLFD"
        T = convert(Array{Float64}, load("../data/tensor_light.jld2", "tensor_light"))
    else
        error("dataset name error")
    end

    return T

end

max_stars_repo_name: gkazunii/Legendre-tucker-rank-reduction

max_issues_repo_name: gkazunii/Legendre-tucker-rank-reduction

max_forks_repo_name: gkazunii/Legendre-tucker-rank-reduction


text:
"""
A 2D Rotation matrix is a mtrix that rotates a vector in a 2D real space by an angle theta.
For more info: https://en.wikipedia.org/wiki/Rotation_matrix

This function takes the angle `theta` in radians as input and returns a 2D Matrix which will rotate the the vector by angle `theta`.
"""
function rotation_matrix(θ)
	rot_mat = Matrix{Float64}(undef,2,2)
	rot_mat[1,1]  = cos(θ)
	rot_mat[1,2] = -sin(θ)
	rot_mat[2,1] = sin(θ)
	rot_mat[2,2] = cos(θ)

	return rot_mat
end

max_stars_repo_name: KohRongSoon/Julia

max_issues_repo_name: KohRongSoon/Julia

max_forks_repo_name: KohRongSoon/Julia


text:
using PtFEM
using Test: @test

data = Dict(
  # Rod(nels, np_types, nip, finite_element(nod, nodof))
  :struc_el => Rod(4, 2, 1, Line(2, 1)),
  :properties => [2.0e3; 1.0e3],
  :etype => [2, 2, 1, 1],
  :x_coords => range(0, stop=1, length=5),
  :support => [
    (1, [0])
    ],
  :fixed_freedoms => [
    (5, 1, 0.05)
    ]
)

@time fem, dis_df, fm_df = p41(data)

println("\nTesting: round.(fem.displacements, digits=7) == [0.0 0.0166667 0.0333333 0.0416667 0.05]' \n")

@test round.(fem.displacements, digits=7) == [0.0 0.0166667 0.0333333 0.0416667 0.05]'


text:
module TestApproxDistinct
    using StreamStats
    using Distributions
    using Base.Test

    for n in rand(1:1_000_000, 100)
        xs = rand(1:n, n)
        stat = StreamStats.ApproxDistinct()
        for x in xs
            update!(stat, x)
        end
        online_d = state(stat)
        online_n = nobs(stat)
        batch_d = length(unique(xs))
        @test abs(online_d - batch_d) / batch_d < 0.03
        @test online_n == n
    end
end


text:
# Useful constants
const DEFAULT_CACHE_DIR = "./__hash_cache__"
const DEFAULT_COMPRESSION = "none"
const DEFAULT_HASHCHAIN_FILENAME = "hashchain.json"
const DEFAULT_HASHCACHE_DIR = "cache"


text:
"""
    Hypergraph{T} <: AbstractMatrix{Union{T, Nothing}}

A hypergraph storing information about vertices and hyperedges.

**Constructors**

    Hypergraph{T}(n::Integer,k::Integer) where {T<:Real}
    Hypergraph{T,V}(n::Integer, k::Integer;
        v_meta=Vector{Union{V,Nothing}}(nothing, n)
        ) where {T<:Real, V}
    Hypergraph{T,V,E}(n::Integer, k::Integer;
        v_meta=Vector{Union{V,Nothing}}(nothing, n),
        he_meta=Vector{Union{E,Nothing}}(nothing, k)
        ) where {T<:Real, V, E}
    Hypergraph{T,V,E,D}(n::Integer, k::Integer,
        v_meta=Vector{Union{V,Nothing}}(nothing, n),
        he_meta=Vector{Union{E,Nothing}}(nothing, k)
        ) where {T<:Real,V,E,D<:AbstractDict{Int,T}}

Construct a hypergraph with a given number of vertices and hyperedges.
Optionally, values of type `V` can be stored at vertices and values of type `E`
can be stored at hyperedges. By default the hypegraph uses a `Dict{Int,T}` for
the internal data storage, however a different dictionary such as `SortedDict`
to ensure result replicability can be used (e.g. when doing stochastic
simulations on hypergraphs).

    Hypergraph(m::AbstractMatrix{Union{T, Nothing}}) where {T<:Real}
    Hypergraph{T, V}(m::AbstractMatrix{Union{T, Nothing}};
        v_meta=Vector{Union{V,Nothing}}(nothing, size(m,1))
        ) where {T<:Real, V}
    Hypergraph{T, V, E}(m::AbstractMatrix{Union{T, Nothing}};
        v_meta=Vector{Union{V,Nothing}}(nothing, size(m,1)),
        he_meta=Vector{Union{E,Nothing}}(nothing, size(m,2))
        ) where {T<:Real, V, E}
    Hypergraph{T, V, E, D}(m::AbstractMatrix{Union{T, Nothing}};
        v_meta=Vector{Union{V,Nothing}}(nothing, size(m,1)),
        he_meta=Vector{Union{E,Nothing}}(nothing, size(m,2))
        ) where {T<:Real, V, E, D<:AbstractDict{Int,T}}

Construct a hypergraph using its matrix representation.
In the matrix representation rows are vertices and columns are hyperedges.
Optionally, values of type `V` can be stored at vertices and values of type `E`
can be stored at hyperedges. By default the hypegraph uses a `Dict{Int,T}` for
the internal data storage, however a different dictionary such as `SortedDict`
to ensure result replicability can be used (e.g. when doing stochastic
simulations on hypergraphs).

    Hypergraph(g::LightGraphs.Graph)

Constructs a hypergraph of degree 2 by making a deep copy of LightGraphs.Graph.
A `SortedDict` will be used for internal data storage of the hypergraph.

**Arguments**

* `T` : type of weight values stored in the hypergraph's adjacency matrix
* `V` : type of values stored in the vertices of the hypergraph
* `E` : type of values stored in the edges of the hypergraph
* `D` : dictionary for storing values the default is `Dict{Int, T}`
* `n` : number of vertices
* `k` : number of hyperedges
* `m` : a matrix representation rows are vertices and columns are hyperedges
"""
struct Hypergraph{T<:Real,V,E,D<:AbstractDict{Int,T}} <: AbstractMatrix{Union{T, Nothing}}
    v2he::Vector{D}
    he2v::Vector{D}
    v_meta::Vector{Union{V,Nothing}}
    he_meta::Vector{Union{E,Nothing}}

    Hypergraph{T,V,E,D}(n::Integer, k::Integer,
        v_meta=Vector{Union{V,Nothing}}(nothing, n),
        he_meta=Vector{Union{E,Nothing}}(nothing, k)
    ) where {T<:Real,V,E,D<:AbstractDict{Int,T}} =
        new{T,V,E,D}(
            [D() for i in 1:n],[D() for i in 1:k],
            v_meta, he_meta)
end

Hypergraph{T,V,E}(n::Integer, k::Integer) where {T<:Real, V, E} = Hypergraph{T,V,E,Dict{Int,T}}(n, k)

Hypergraph{T,V}(n::Integer, k::Integer) where {T<:Real, V, E} = Hypergraph{T,V,Nothing,Dict{Int,T}}(n, k)

Hypergraph{T}(n::Integer, k::Integer) where {T<:Real} =  Hypergraph{T,Nothing,Nothing,Dict{Int,T}}(n, k)

Hypergraph(n::Integer, k::Integer) =  Hypergraph{Bool,Nothing,Nothing,Dict{Int,Bool}}(n, k)

function Hypergraph{T,V,E,D}(m::AbstractMatrix{Union{T, Nothing}};
                        v_meta::Vector{Union{Nothing,V}}=Vector{Union{Nothing,V}}(nothing, size(m,1)),
                        he_meta::Vector{Union{Nothing,E}}=Vector{Union{Nothing,E}}(nothing, size(m,2))
                        ) where {T<:Real,V,E,D<:AbstractDict{Int,T}}
    @assert length(v_meta) == size(m,1)
    @assert length(he_meta) == size(m,2)
    n, k = size(m)
    h = Hypergraph{T,V,E,D}(n, k, v_meta, he_meta)
    h .= m
    h
end

function Hypergraph{T,V,E}(m::AbstractMatrix{Union{T, Nothing}};
                        v_meta::Vector{Union{Nothing,V}}=Vector{Union{Nothing,V}}(nothing, size(m,1)),
                        he_meta::Vector{Union{Nothing,E}}=Vector{Union{Nothing,E}}(nothing, size(m,2))
                        ) where {T<:Real,V,E}
    Hypergraph{T,V,E,Dict{Int,T}}(m;v_meta=v_meta,he_meta=he_meta)
end

function Hypergraph{T,V}(m::AbstractMatrix{Union{T, Nothing}};
                        v_meta::Vector{Union{Nothing,V}}=Vector{Union{Nothing,V}}(nothing, size(m,1))
                        ) where {T<:Real,V}
    Hypergraph{T,V,Nothing,Dict{Int,T}}(m;v_meta=v_meta)
end

function Hypergraph{T}(m::AbstractMatrix{Union{T, Nothing}}) where {T<:Real}
    Hypergraph{T,Nothing,Nothing,Dict{Int,T}}(m)
end

function Hypergraph(m::AbstractMatrix{Union{T, Nothing}}) where {T<:Real}
    Hypergraph{T,Nothing,Nothing,Dict{Int,T}}(m)
end

function Hypergraph(g::LightGraphs.Graph)
    h = Hypergraph{Bool,Nothing,Nothing,SortedDict{Int,Bool}}(maximum(vertices(g)), ne(g))
    e = 0
    for edge in edges(g)
        e+=1
        h[edge.src,e] = true
        h[edge.dst,e] = true
    end
    h
end

"""
    Base.size(h::Hypergraph)

Returns the size of Hypergraph `h`.
The result is a tuple of the number of vertices and the number of hyperedges

"""
Base.size(h::Hypergraph) = (nhv(h), nhe(h))

"""
    Base.getindex(h::Hypergraph, idx::Vararg{Int,2})

Returns a value for a given vertex-hyperedge pair `idx` for a hypergraph `h`.
If a vertex does not belong to a hyperedge `nothing` is returned.

"""
@inline function Base.getindex(h::Hypergraph, idx::Vararg{Int,2})
    @boundscheck checkbounds(h, idx...)
    get(h.v2he[idx[1]], idx[2], nothing)
end

"""
    Base.setindex!(h::Hypergraph, ::Nothing, idx::Vararg{Int,2})

Removes a vertex from a given hyperedge for a hypergraph `h` and a given vertex-hyperedge pair `idx`.
Note that trying to remove a vertex from a hyperedge when it is not present will not throw an error.

"""
@inline function Base.setindex!(h::Hypergraph, ::Nothing, idx::Vararg{Int,2})
    @boundscheck checkbounds(h, idx...)
    pop!(h.v2he[idx[1]], idx[2], nothing)
    pop!(h.he2v[idx[2]], idx[1], nothing)
    h
end

"""
    Base.setindex!(h::Hypergraph, v::Real, idx::Vararg{Int,2})

Adds a vertex to a hyperedge (represented by indices `idx`) and assigns value
`v` to be stored with that assignment.

"""
@inline function Base.setindex!(h::Hypergraph, v::Real, idx::Vararg{Int,2})
    @boundscheck checkbounds(h, idx...)
    h.v2he[idx[1]][idx[2]] = v
    h.he2v[idx[2]][idx[1]] = v
    h
end

"""
    getvertices(h::Hypergraph, he_id::Int)

Returns vertices from a hypergraph `a` for a given hyperedge `he_id`.

"""
@inline getvertices(h::Hypergraph, he_id::Int) = h.he2v[he_id]

"""
    gethyperedges(h::Hypergraph, v_id::Int)

Returns hyperedges for a given vertex `v_id` in a hypergraph `h`.

"""
@inline gethyperedges(h::Hypergraph, v_id::Int) = h.v2he[v_id]

"""
    add_vertex!(h::Hypergraph{T, V, E, D};
                hyperedges::D = D(), v_meta::Union{V,Nothing} = nothing
                ) where {T <: Real, V, E, D <: AbstractDict{Int,T}}

Adds a vertex to a given hypergraph `h`. Optionally, the vertex can be added
to existing hyperedges. The `hyperedges` parameter presents a dictionary
of hyperedge identifiers and values stored at the hyperedges.
Additionally, a value can be stored with the vertex using the `v_meta` keyword
parameter.

"""
function add_vertex!(h::Hypergraph{T, V, E, D};
                     hyperedges::D = D(), v_meta::Union{V,Nothing} = nothing
                    ) where {T <: Real, V, E, D <: AbstractDict{Int,T}}
    @boundscheck (checkbounds(h,1,k) for k in keys(hyperedges))
    push!(h.v2he,hyperedges)
    ix = length(h.v2he)
    for k in keys(hyperedges)
        h[ix,k]=hyperedges[k]
    end
    push!(h.v_meta, v_meta)
    ix
end

"""
    remove_vertex!(h::Hypergraph, v::Int)

Removes the vertex `v` from a given hypergraph `h`.
Note that running this function will cause reordering of vertices in the
hypergraph: the vertex `v` will replaced by the last vertex of the hypergraph
and the list of vertices will be shrunk.
"""
function remove_vertex!(h::Hypergraph, v::Int)
    n = nhv(h)
    if v < n
        h.v2he[v] = h.v2he[n]
        h.v_meta[v] = h.v_meta[n]
    end

    for hv in h.he2v
        if v < n && haskey(hv, n)
            hv[v] = hv[n]
            delete!(hv, n)
        else
            delete!(hv, v)
        end
    end
    resize!(h.v2he, length(h.v2he) - 1)
    h
end

"""
    add_hyperedge!(h::Hypergraph{T, V, E, D};
                   vertices::D = D(), he_meta::Union{E,Nothing}=nothing
                   ) where {T <: Real, V, E, D <: AbstractDict{Int,T}}

Adds a hyperedge to a given hypergraph `h`.
Optionally, existing vertices can be added to the created hyperedge.
The paramater `vertices` represents a dictionary of vertex identifiers and
values stored at the hyperedges. Additionally, a value can be stored with the
hyperedge using the `he_meta` keyword parameter.

"""
function add_hyperedge!(h::Hypergraph{T, V, E, D};
                        vertices::D = D(), he_meta::Union{E,Nothing}=nothing
                        ) where {T <: Real, V, E, D <: AbstractDict{Int,T}}
    @boundscheck (checkbounds(h,k,1) for k in keys(vertices))
    push!(h.he2v,vertices)
    ix = length(h.he2v)
    for k in keys(vertices)
        h[k,ix]=vertices[k]
    end
    push!(h.he_meta, he_meta)
    ix
end

"""
    remove_hyperedge!(h::Hypergraph, e::Int)
Removes the heyperedge `e` from a given hypergraph `h`.
Note that running this function will cause reordering of hyperedges in the
hypergraph: the hyperedge `e` will replaced by the last hyperedge of the hypergraph
and the list of hyperedges will be shrunk.
"""
function remove_hyperedge!(h::Hypergraph, e::Int)
    ne = nhe(h)
	@assert(e <= ne)
	if e < ne
	    h.he2v[e] = h.he2v[ne]
	    h.he_meta[e] = h.he_meta[ne]
	end

    for he in h.v2he
	    if e < ne && haskey(he, ne)
		    he[e] = he[ne]
            delete!(he, ne)
		else
			delete!(he, e)
		end
    end
    resize!(h.he2v, length(h.he2v) - 1)
    h
end

"""
    prune_hypergraph!(h::Hypergraph)

Remove all vertices with degree 0 and all hyperedges of size 0.

"""
function prune_hypergraph!(h::Hypergraph)
	for e in reverse(1:nhe(h))
        length(h.he2v[e]) == 0 && remove_hyperedge!(h,e)
    end
	for v in reverse(1:nhv(h))
    	length(h.v2he[v]) == 0 && 	remove_vertex!(h,v)
    end
	h
end

"""
    prune_hypergraph(h::Hypergraph)

Return a pruned copy of `h`, removing all vertices with degree 0 and
all hyperedges of size 0.

"""
function prune_hypergraph(h::Hypergraph)
    prune_hypergraph!(deepcopy(h))
end

"""
    set_vertex_meta!(h::Hypergraph{T, V, E, D}, new_value::Union{V,Nothing},
        id::Int) where {T <: Real, V, E, D <: AbstractDict{Int,T}}

Sets a new meta value `new_value` for the vertex `id` in the hypegraph `h`.

"""
function set_vertex_meta!(h::Hypergraph{T, V, E, D},
        new_value::Union{V,Nothing}, id::Int
        ) where {T <: Real, V, E, D <: AbstractDict{Int,T}}
    checkbounds(h.v_meta, id)
    h.v_meta[id] = new_value
    h.v_meta
end

"""
    get_vertex_meta(h::Hypergraph{T, V, E, D}, id::Int
                    ) where {T <: Real, V, E, D <: AbstractDict{Int,T}}

Returns a meta value stored at the vertex `id` in the hypergraph `h`.

"""
function get_vertex_meta(h::Hypergraph{T, V, E, D}, id::Int
                         ) where {T <: Real, V, E, D <: AbstractDict{Int,T}}
    checkbounds(h.v_meta, id)
    h.v_meta[id]
end

"""
    set_hyperedge_meta!(h::Hypergraph{T, V, E, D},
        new_value::Union{E,Nothing}, id::Int
        ) where {T <: Real, V, E, D <: AbstractDict{Int,T}}

Sets a new meta value `new_value` for the hyperedge `id` in the hypegraph `h`.

"""
function set_hyperedge_meta!(h::Hypergraph{T, V, E, D},
                            new_value::Union{E,Nothing}, id::Int
                            ) where {T <: Real, V, E, D <: AbstractDict{Int,T}}
    checkbounds(h.he_meta, id)
    h.he_meta[id] = new_value
    h.he_meta
end

"""
    get_hyperedge_meta(h::Hypergraph{T, V, E, D}, id::Int)
        where {T <: Real, V, E, D <: AbstractDict{Int,T}}
Returns a meta value stored at the hyperedge `id` in the hypergraph `h`.

"""
function get_hyperedge_meta(h::Hypergraph{T, V, E, D}, id::Int
                            ) where {T <: Real, V, E, D <: AbstractDict{Int,T}}
    checkbounds(h.he_meta, id)
    h.he_meta[id]
end

"""
    nhe(h::Hypergraph)

Return the number of hyperedges in the hypergraph `h`.
"""
function nhe(h::Hypergraph)
    length(h.he2v)
end

"""
    nhv(h::Hypergraph)

Return the number of vertices in the hypergraph `h`.
"""
function nhv(h::Hypergraph)
    length(h.v2he)
end

function _default_heselect(h::Hypergraph, v::Int)
    hes = gethyperedges(h, v)
    sort!(collect(keys(hes))), ones(length(hes))
end

function _default_vselect(h::Hypergraph, he::Int)
    vs = getvertices(h, he)
    sort!(collect(keys(vs))), ones(length(vs))
end

"""
    random_walk(h::Hypergraph, start::Int; heselect::Function, vselect::Function)

Return a next vertex visited in assuming a random walk starting from vertex `start`.
First a hyperedge is sampled with weights proportional to `heselect` function
(by default each hyperedge is sampled with the same probability).
Next a vertex within hyperedge is with weights proportional to `vselect` function
(by default each vertex, including the source, is sampled with the same probability).

`heselect` and `vselect` functions take two arguments a `Hypergraph` and respectively
a vertex identifier or a hyperedge identifier. The return values of both functions
should be respectively a list of hyperedges or vertices and their weights.
"""
function random_walk(h::Hypergraph, start::Int;
                     heselect::Function=_default_heselect,
                     vselect::Function=_default_vselect)
    1 <= start <= nhv(h) || throw(ArgumentError("invalid start vertex index"))
    hes, hew = heselect(h, start)
    he = sample(hes, Weights(hew))
    ves, vw = vselect(h, he)
    return sample(ves, Weights(vw))
end

"""
    _walk!(h::Hypergraph, s::AbstractVector{Int}, i::Int, visited::AbstractVector{Bool})

Appends the list of neighbors `s` of a given vertex `i` (an auxiliary function for `get_connected_components`).
"""
function _walk!(h::Hypergraph, s::AbstractVector{Int}, i::Int, visited::AbstractVector{Bool})
    visited[i] && return
    visited[i] = true
    push!(s, i)
    for he in keys(gethyperedges(h, i))
        for j in keys(getvertices(h, he))
            _walk!(h, s, j, visited)
        end
    end
end

"""
    get_connected_components(h::Hypergraph)

Return an array of connected components in the hypergraph `h`
(array of vectors of vertices) using recurrence.
"""
function get_connected_components(h::Hypergraph)
    visited = falses(nhv(h))
    cc = Vector{Int}[]
        for i in 1:nhv(h)
            if !visited[i]
                s = Int[]
                _walk!(h, s, i, visited)
                push!(cc, s)
        end
    end
    cc
end

"""
    adjacency_matrix(h::Hypergraph; s::Int=1, weighted::Bool=false)

The sparse weighted `s`-adjacency matrix.

NOTE
The concept of `s`-adjacency matrix has been firstly defined in the
Python library [HyperNetX](https://github.com/pnnl/HyperNetX)

From [HyperNetX](https://pnnl.github.io/HyperNetX/build/classes/classes.html#classes.hypergraph.Hypergraph.adjacency_matrix)
If weighted is `true` each off diagonal cell will equal the number
of edges shared by the nodes indexing the row and column if that number is
greater than `s`, otherwise the cell will equal 0. If weighted is `false`,
the off diagonal cell will equal 1 if the nodes indexed by the row and column
share at least `s` edges and 0 otherwise.

! information about the weight of a vertex in a he will be lost.

"""
function adjacency_matrix(h; s::Int=1, weighted::Bool=true)
    M = Matrix(h)
    _incidence_to_adjacency(M; s=s, weighted=weighted)
end

"""
    edge_adjacency_matrix(h::Hypergraph; s::Int=1, weighted::Bool=false)

The sparse weighted `s`-adjacency matrix for the dual hypergraph.

NOTE
The concept of `s`-adjacency matrix has been firstly defined in the
Python library [HyperNetX](https://github.com/pnnl/HyperNetX)

From [HyperNetX](https://pnnl.github.io/HyperNetX/build/classes/classes.html#classes.hypergraph.Hypergraph.edge_adjacency_matrix)
This is also the adjacency matrix for the line graph.
Two edges are `s`-adjacent if they share at least `s` nodes.

If weighted is `true` each off diagonal cell will equal the number
of nodes shared by the hyperedges indexing the row and column if that number is
greater than `s`, otherwise the cell will equal 0. If weighted is `false`,
the off diagonal cell will equal 1 if the hyperedges indexed by the row and column
share at least `s` nodes and 0 otherwise.

"""
function edge_adjacency_matrix(h; s::Int=1, weighted::Bool=true)
    M = Matrix(h)
	M[M .== nothing] .= 0
    _incidence_to_adjacency(transpose(M); s=s, weighted=weighted)
end

"""
    _incidence_to_adjacency(M; s::Int=1, weighted::Bool=true)

Helper method to obtain adjacency matrix from incidence matrix.

"""
function _incidence_to_adjacency(M; s::Int=1, weighted::Bool=true)
    M[M .== nothing] .= 0
    M[M .> 0] .= 1

    A = *(M, transpose(M))
    A[diagind(A)] .= 0

    if s > 1
        A = A .* (A .>= s)
    end
    if !weighted
        A = (A .> 0) .* 1
    end

    A
end

# TODO find connected components without recurrence
# TODO needs validate_hypergraph!(h::Hypergraph{T})


text:
function grouped_bar_xy_test()
    nx = 3;
    ny = 4;
    dataM = range(1.0, 2.0, length = nx) * range(1.0, 1.5, length = ny)';
    yErrM = 0.1 .* dataM;
    groupLabels = ["Group $j"  for j = 1 : ny];
    legendV = ["x $j" for j = 1 : nx];

    fName = "grouped_bar_xy_test";
    delete_test_figure(fName);
    p = grouped_bar_xy(dataM, groupLabels; xLabel = "x axis", yLabel = "y axis",
        legendV = legendV, yErrorM = yErrM);
    save_test_figure(p, fName);
    @test isfile(test_fig_path(fName));
end

@testset "XY" begin
    grouped_bar_xy_test();
end

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