text:
# CSCI5573Final/schedulers.jl
# Licensed under the MIT License. See LICENSE.md file in the project root for
# full license information.

# -------------------- Begin Shared Functions -----------------------------
function EST(ni, pj, aft)
    available = current_time() + TotalWork(PROCESSORS[pj])
    max_pred = 0

    for parent in parentsof(ni, TASKS)
        cij = COMMS[Edge(ni, parent)]
        if parent in PROCESSORS[pj].queue || parent === PROCESSORS[pj].task
            cij = 0
        end
        new_pred = aft[parent] + cij
        if new_pred > max_pred
            max_pred = new_pred
        end
    end

    return max(available, max_pred)
end

function AFT(ranks, processors)
    aft = Dict{Int64,Float64}()

    for (ID, rank) in ranks
        comm_cost = 0
        if !isempty(TASKS[ID].Children)
            comm_cost = sum(l->COMMS[Edge(ID,l)], TASKS[ID].Children)
        end

        aft[ID] = max(((p->p.Multiplier*rank).(processors))...) + comm_cost |> round
    end

    return aft
end

EFT(t, p, aft) = EST(t, p, aft) + WIJ(t, p)

function WIJ(task, processor)
    return PROCESSORS[processor].Multiplier * TASKS[task].Cost
end

function CommCost(ni, pj)
    comm_cost = 0
    if !isempty(parentsof(ni, TASKS))
        for parent in parentsof(ni, TASKS)
            if parent in PROCESSORS[pj].queue || parent === PROCESSORS[pj].task
                cij = 0
            else
                comm_cost += COMMS[Edge(ni, parent)]
            end
        end
    end
    return comm_cost
end

Cost(ID) = TASKS[ID].Cost
# ------------------------- End Functions ---------------------------------

# ----------------------- Begin List Scheduler -----------------------
@process ProposedScheduler() begin
    global PROCESSORS
    global SCHEDULER = "ProposedScheduler"

    rank = ProposedRank(TASKS)
    aft = AFT(rank, PROCESSORS)

    while(true)
        if isempty(READY_QUEUE)
            wait(CLOCK_CYCLE)
            continue
        end

        ready_list = []
        max_rank   = Inf * -1
        min_EFT    = Inf
        pj         = -1
        ni         = -1

        # Empty the queue to start scheduling
        while(!isempty(READY_QUEUE))
            ID = dequeue!(READY_QUEUE)
            push!(ready_list, ID)
        end

        while(length(ready_list) > 0)
            # Max Rank Task
            for task in ready_list
                if rank[task] > max_rank
                    max_rank = rank[task]
                    ni = task
                end
            end

            # Minimum EFT Processor
            for j in 1:length(PROCESSORS)
                if EFT(ni, j, aft) < min_EFT
                    min_EFT = EFT(ni, j, aft)
                    pj = j
                end
            end

            wik = min(map(k->WIJ(ni, k), 1:length(PROCESSORS))...)

            if WIJ(ni, pj) <= wik
                # Schedule this task to start on that processor
                @schedule now Dispatcher(ni, pj, Cost(ni), CommCost(ni, pj))
            else
                pk = filter(k->WIJ(ni, k) === wik, 1:length(PROCESSORS))[begin]

                wa_numerator    = EFT(ni, pj, aft) - EFT(ni, pk, aft)
                wa_denominator  = EFT(ni, pj, aft) / EFT(ni, pk, aft)

                weight_ni       = WIJ(ni, pj) / TASKS[ni].Cost

                weight_abstract = abs(wa_numerator/wa_denominator)
                cross_threshold = abs(weight_ni/weight_abstract)

                if cross_threshold <= 1-rand()
                    @schedule now Dispatcher(ni, pj, Cost(ni), CommCost(ni, pj))
                else
                    @schedule now Dispatcher(ni, pk, Cost(ni), CommCost(ni, pj))
                end
            end
            filter!(nj->nj!==ni, ready_list)
            work(CLOCK_CYCLE)
        end
    end
end

function ProposedRank(task_graph)
    ranks = Dict{Int64,Float64}()

    RecursiveRankTask(ID, ranks) = begin
        task = task_graph[ID]

        rank = 0.0
        for child in task.Children
            new_rank = RecursiveRankTask(child, ranks) + COMMS[Edge(ID, child)]
            rank = max(new_rank, rank)
        end

        ranks[ID] = task.Cost + rank
        return ranks[ID]
    end

    RecursiveRankTask(0, ranks)

    return ranks
end

# ----------------------- End Proposed -------------------------

# ----------------------------- Begin HEFT -----------------------------
@process HEFTScheduler() begin
    global SCHEDULER = "HEFTScheduler"
    rank = UpwardRank(TASKS, COMMS)
    aft = AFT(rank, PROCESSORS)

    while(true)
        if isempty(READY_QUEUE)
            wait(CLOCK_CYCLE)
            continue
        end

        ready_list = []
        min_EFT    = Inf
        pj         = -1
        ni         = -1

        while(!isempty(READY_QUEUE)) # Empty the queue to start scheduling
            ID = dequeue!(READY_QUEUE)
            push!(ready_list, ID)
        end

        while(length(ready_list) > 0)
            ready_list = sort(ready_list, lt=(a, b) -> rank[a] > rank[b])
            ni = ready_list[1]
            deleteat!(ready_list, 1)

            for j in 1:length(PROCESSORS)
                if EFT(ni, j, aft) < min_EFT
                    min_EFT = EFT(ni, j, aft)
                    pj = j
                end
            end

            comm_cost = 0
            if !isempty(parentsof(ni, TASKS))
                for parent in parentsof(ni, TASKS)
                    if parent in PROCESSORS[pj].queue || parent === PROCESSORS[pj].task
                        cij = 0
                    else
                        comm_cost += COMMS[Edge(ni, parent)]
                    end
                end
            end

            @schedule now Dispatcher(ni, pj, Cost(ni), CommCost(ni, pj))
        end

        work(CLOCK_CYCLE)
    end
end

function UpwardRank(task, communication)
    ranks = Dict{Int64, Int64}()

    RecursiveRankTask(ID, ranks) = begin

        rank = 0.0
        for child in task[ID].Children
            new_rank = communication[Edge(ID, child)] + RecursiveRankTask(child, ranks)

            rank = max(new_rank, rank)
        end
        ranks[ID] = task[ID].Cost + rank
    end

    RecursiveRankTask(0, ranks)

    return ranks
end

# ----------------------------- End HEFT -----------------------------

# ----------------------------- Begin PEFT -----------------------------
@process PEFTScheduler() begin
    global SCHEDULER = "PEFTScheduler"
    oct = OCT(TASKS)
    rank_oct = RankOCT(oct, TASKS)
    aft = AFT(rank_oct, PROCESSORS)

    while(true)
        if isempty(READY_QUEUE)
            wait(CLOCK_CYCLE)
            continue
        end

        ready_list  = []
        min_OEFT    = Inf
        pj          = -1
        ni          = -1

        # Empty the queue to start scheduling
        while(!isempty(READY_QUEUE))
            ID = dequeue!(READY_QUEUE)
            push!(ready_list, ID)
        end

        while(length(ready_list) > 0)
            ready_list = sort(ready_list, lt=(a, b) -> rank_oct[a] > rank_oct[b])
            ni = ready_list[1]
            deleteat!(ready_list, 1)

            OEFT(n, p) = EFT(n, p, aft) + oct[n][p]

            for j in 1:length(PROCESSORS)
                if OEFT(ni, j) < min_OEFT
                    min_OEFT = OEFT(ni, j)
                    pj = j
                end
            end

            @schedule now Dispatcher(ni, pj, Cost(ni), CommCost(ni, pj))
        end
        work(CLOCK_CYCLE)
    end
end
function OCT(tasks) # Paper used recursion but... it blows the stack.
    exit_task = max(unique(collect(keys(tasks)))...)
    oct = Dict{Int64, Array{Int64, 1}}()

    oct[exit_task] = repeat([0], length(PROCESSORS)) # Set Exit

    predecessors = parentsof(exit_task, tasks)
    while length(predecessors) > 0
        ni = pop!(predecessors)

        while !all(map(s->haskey(oct, s), tasks[ni].Children))
            nj = -1
            try
                nj = pop!(predecessors)
            catch
                println("An error occured during OCT calculation")
                sleep(5)
                return
            end
            insert!(predecessors, 1, ni)
            ni = nj
        end

        oct[ni] = repeat([0], length(PROCESSORS))

        for pj in 1:length(PROCESSORS)
            my_oct = -Inf
            for child in tasks[ni].Children
                min_child_oct = Inf
                for pw in 1:length(PROCESSORS)
                    child_oct = oct[child][pw]
                    child_cost = tasks[child].Cost
                    child_comm_cost = pw === pj ? 0 : COMMS[Edge(ni, child)]
                    new_min_oct = child_oct + child_cost + child_comm_cost
                    if new_min_oct < min_child_oct
                        min_child_oct = new_min_oct
                    end
                end
                if min_child_oct > my_oct
                    my_oct = min_child_oct
                end
            end
            oct[ni][pj] = my_oct
        end
        predecessors = append!(filter(f->!(f in predecessors), parentsof(ni, tasks)), predecessors)
    end

    return oct
end

function RankOCT(oct, tasks)
    ranks = Dict{Int64, Int64}()
    for (ti, _) in tasks
        ranks[ti] = sum(map(pj->oct[ti][pj], 1:length(PROCESSORS))) / length(PROCESSORS)
    end
    return ranks
end

# ----------------------------- End PEFT -----------------------------

# ----------------------------- Begin FCFS -----------------------------
@process FCFSScheduler() begin
    global SCHEDULER = "FCFSScheduler"
    processor = 0
    while(true)
        if !isempty(READY_QUEUE)
            ID            = dequeue!(READY_QUEUE)
            parents       = parentsof(ID, TASKS)
            comm_cost     = sum(map(p->COMMS[Edge(p, ID)], parents))

            @schedule now Dispatcher(ID, processor, Cost(ID), CommCost(ID, processor))
        end

        processor = processor === length(PROCESSORS) ? 1 : processor + 1

        work(CLOCK_CYCLE)
    end
end
# ----------------------------- End FCFS -----------------------------


text:
"""
    to_vec(x)

Transform `x` into a `Vector`, and return the vector, and a closure which inverts the
transformation.
"""
function to_vec(x::Real)
    function Real_from_vec(x_vec)
        return first(x_vec)
    end
    return [x], Real_from_vec
end

function to_vec(z::Complex)
    function Complex_from_vec(z_vec)
        return Complex(z_vec[1], z_vec[2])
    end
    return [real(z), imag(z)], Complex_from_vec
end

# Base case -- if x is already a Vector{<:Real} there's no conversion necessary.
to_vec(x::Vector{<:Real}) = (x, identity)

function to_vec(x::AbstractVector)
    x_vecs_and_backs = map(to_vec, x)
    x_vecs, backs = first.(x_vecs_and_backs), last.(x_vecs_and_backs)
    function Vector_from_vec(x_vec)
        sz = cumsum(map(length, x_vecs))
        x_Vec = [backs[n](x_vec[sz[n] - length(x_vecs[n]) + 1:sz[n]]) for n in eachindex(x)]
        return oftype(x, x_Vec)
    end
    return vcat(x_vecs...), Vector_from_vec
end

function to_vec(x::AbstractArray)

    x_vec, from_vec = to_vec(vec(x))

    function Array_from_vec(x_vec)
        return oftype(x, reshape(from_vec(x_vec), size(x)))
    end

    return x_vec, Array_from_vec
end

# Some specific subtypes of AbstractArray.

function to_vec(x::T) where {T<:LinearAlgebra.AbstractTriangular}
    x_vec, back = to_vec(Matrix(x))
    function AbstractTriangular_from_vec(x_vec)
        return T(reshape(back(x_vec), size(x)))
    end
    return x_vec, AbstractTriangular_from_vec
end

function to_vec(x::T) where {T<:LinearAlgebra.HermOrSym}
    x_vec, back = to_vec(Matrix(x))
    function HermOrSym_from_vec(x_vec)
        return T(back(x_vec), x.uplo)
    end
    return x_vec, HermOrSym_from_vec
end

function to_vec(X::Diagonal)
    x_vec, back = to_vec(Matrix(X))
    function Diagonal_from_vec(x_vec)
        return Diagonal(back(x_vec))
    end
    return x_vec, Diagonal_from_vec
end

function to_vec(X::Transpose)
    x_vec, back = to_vec(Matrix(X))
    function Transpose_from_vec(x_vec)
        return Transpose(permutedims(back(x_vec)))
    end
    return x_vec, Transpose_from_vec
end

function to_vec(X::Adjoint)
    x_vec, back = to_vec(Matrix(X))
    function Adjoint_from_vec(x_vec)
        return Adjoint(conj!(permutedims(back(x_vec))))
    end
    return x_vec, Adjoint_from_vec
end

function to_vec(X::PermutedDimsArray{T, D, P, Pinv} where {T, D}) where {P, Pinv}
    x_vec, back = to_vec(collect(X))
    function PermutedDimsArray_from_vec(x_vec)
        X_parent = collect(PermutedDimsArray(back(x_vec), Pinv))
        return PermutedDimsArray(X_parent, P)
    end
    return x_vec, PermutedDimsArray_from_vec
end

# Non-array data structures

function to_vec(x::Tuple)
    x_vecs, x_backs = zip(map(to_vec, x)...)
    sz = cumsum(collect(map(length, x_vecs)))
    function Tuple_from_vec(v)
        return ntuple(length(x)) do n
            return x_backs[n](v[sz[n] - length(x_vecs[n]) + 1:sz[n]])
        end
    end
    return reduce(vcat, x_vecs), Tuple_from_vec
end

function to_vec(x::NamedTuple)
    x_vec, back = to_vec(values(x))
    function NamedTuple_from_vec(v)
        v_vec_vec = back(v)
        return typeof(x)(v_vec_vec)
    end
    return x_vec, NamedTuple_from_vec
end

# Convert to a vector-of-vectors to make use of existing functionality.
function to_vec(d::Dict)
    d_vec, back = to_vec(collect(values(d)))
    function Dict_from_vec(v)
        v_vec_vec = back(v)
        return Dict(key => v_vec_vec[n] for (n, key) in enumerate(keys(d)))
    end
    return d_vec, Dict_from_vec
end

ext: jl


text:
# Julia wrapper for header: /Users/yeesian/Downloads/libspatialindex_c-master/include/spatialindex/capi/sidx_api.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

function Index_Create(properties::IndexPropertyH)
    ccall((:Index_Create, libspatialindex_c), IndexH, (IndexPropertyH,), properties)
end

function Index_CreateWithStream(properties::IndexPropertyH, readNext)
    ccall((:Index_CreateWithStream, libspatialindex_c), IndexH, (IndexPropertyH, Ptr{Nothing}), properties, readNext)
end

function Index_Destroy(index::IndexH)
    ccall((:Index_Destroy, libspatialindex_c), Nothing, (IndexH,), index)
end

function Index_GetProperties(index::IndexH)
    ccall((:Index_GetProperties, libspatialindex_c), IndexPropertyH, (IndexH,), index)
end

function Index_DeleteData(index::IndexH, id::Integer, pdMin, pdMax, nDimension::Integer)
    ccall((:Index_DeleteData, libspatialindex_c), RTError, (IndexH, Int64, Ptr{Cdouble}, Ptr{Cdouble}, UInt32), index, id, pdMin, pdMax, nDimension)
end

function Index_DeleteTPData(index::IndexH, id::Integer, pdMin, pdMax, pdVMin, pdVMax, tStart::Cdouble, tEnd::Cdouble, nDimension::Integer)
    ccall((:Index_DeleteTPData, libspatialindex_c), RTError, (IndexH, Int64, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Cdouble, Cdouble, UInt32), index, id, pdMin, pdMax, pdVMin, pdVMax, tStart, tEnd, nDimension)
end

function Index_DeleteMVRData(index::IndexH, id::Integer, pdMin, pdMax, tStart::Cdouble, tEnd::Cdouble, nDimension::Integer)
    ccall((:Index_DeleteMVRData, libspatialindex_c), RTError, (IndexH, Int64, Ptr{Cdouble}, Ptr{Cdouble}, Cdouble, Cdouble, UInt32), index, id, pdMin, pdMax, tStart, tEnd, nDimension)
end

function Index_InsertData(index::IndexH, id::Integer, pdMin, pdMax, nDimension::Integer, pData, nDataLength::Integer)
    ccall((:Index_InsertData, libspatialindex_c), RTError, (IndexH, Int64, Ptr{Cdouble}, Ptr{Cdouble}, UInt32, Ptr{UInt8}, Cint), index, id, pdMin, pdMax, nDimension, pData, nDataLength)
end

function Index_InsertTPData(index::IndexH, id::Integer, pdMin, pdMax, pdVMin, pdVMax, tStart::Cdouble, tEnd::Cdouble, nDimension::Integer, pData, nDataLength::Integer)
    ccall((:Index_InsertTPData, libspatialindex_c), RTError, (IndexH, Int64, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Cdouble, Cdouble, UInt32, Ptr{UInt8}, Cint), index, id, pdMin, pdMax, pdVMin, pdVMax, tStart, tEnd, nDimension, pData, nDataLength)
end

function Index_InsertMVRData(index::IndexH, id::Integer, pdMin, pdMax, tStart::Cdouble, tEnd::Cdouble, nDimension::Integer, pData, nDataLength::Integer)
    ccall((:Index_InsertMVRData, libspatialindex_c), RTError, (IndexH, Int64, Ptr{Cdouble}, Ptr{Cdouble}, Cdouble, Cdouble, UInt32, Ptr{UInt8}, Cint), index, id, pdMin, pdMax, tStart, tEnd, nDimension, pData, nDataLength)
end

function Index_IsValid(index::IndexH)
    ccall((:Index_IsValid, libspatialindex_c), UInt32, (IndexH,), index)
end

function Index_TPIntersects_obj(index::IndexH, pdMin, pdMax, pdVMin, pdVMax, tStart::Cdouble, tEnd::Cdouble, nDimension::Integer, items, nResults)
    ccall((:Index_TPIntersects_obj, libspatialindex_c), RTError, (IndexH, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Cdouble, Cdouble, UInt32, Ptr{Ptr{IndexItemH}}, Ptr{UInt64}), index, pdMin, pdMax, pdVMin, pdVMax, tStart, tEnd, nDimension, items, nResults)
end

function Index_MVRIntersects_obj(index::IndexH, pdMin, pdMax, tStart::Cdouble, tEnd::Cdouble, nDimension::Integer, items, nResults)
    ccall((:Index_MVRIntersects_obj, libspatialindex_c), RTError, (IndexH, Ptr{Cdouble}, Ptr{Cdouble}, Cdouble, Cdouble, UInt32, Ptr{Ptr{IndexItemH}}, Ptr{UInt64}), index, pdMin, pdMax, tStart, tEnd, nDimension, items, nResults)
end

function Index_Intersects_obj(index::IndexH, pdMin, pdMax, nDimension::Integer, items, nResults)
    ccall((:Index_Intersects_obj, libspatialindex_c), RTError, (IndexH, Ptr{Cdouble}, Ptr{Cdouble}, UInt32, Ptr{Ptr{IndexItemH}}, Ptr{UInt64}), index, pdMin, pdMax, nDimension, items, nResults)
end

function Index_TPIntersects_id(index::IndexH, pdMin, pdMax, pdVMin, pdVMax, tStart::Cdouble, tEnd::Cdouble, nDimension::Integer, ids, nResults)
    ccall((:Index_TPIntersects_id, libspatialindex_c), RTError, (IndexH, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Cdouble, Cdouble, UInt32, Ptr{Ptr{Int64}}, Ptr{UInt64}), index, pdMin, pdMax, pdVMin, pdVMax, tStart, tEnd, nDimension, ids, nResults)
end

function Index_MVRIntersects_id(index::IndexH, pdMin, pdMax, tStart::Cdouble, tEnd::Cdouble, nDimension::Integer, ids, nResults)
    ccall((:Index_MVRIntersects_id, libspatialindex_c), RTError, (IndexH, Ptr{Cdouble}, Ptr{Cdouble}, Cdouble, Cdouble, UInt32, Ptr{Ptr{Int64}}, Ptr{UInt64}), index, pdMin, pdMax, tStart, tEnd, nDimension, ids, nResults)
end

function Index_Intersects_id(index::IndexH, pdMin, pdMax, nDimension::Integer, items, nResults)
    ccall((:Index_Intersects_id, libspatialindex_c), RTError, (IndexH, Ptr{Cdouble}, Ptr{Cdouble}, UInt32, Ptr{Ptr{Int64}}, Ptr{UInt64}), index, pdMin, pdMax, nDimension, items, nResults)
end

function Index_TPIntersects_count(index::IndexH, pdMin, pdMax, pdVMin, pdVMax, tStart::Cdouble, tEnd::Cdouble, nDimension::Integer, nResults)
    ccall((:Index_TPIntersects_count, libspatialindex_c), RTError, (IndexH, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Cdouble, Cdouble, UInt32, Ptr{UInt64}), index, pdMin, pdMax, pdVMin, pdVMax, tStart, tEnd, nDimension, nResults)
end

function Index_MVRIntersects_count(index::IndexH, pdMin, pdMax, tStart::Cdouble, tEnd::Cdouble, nDimension::Integer, nResults)
    ccall((:Index_MVRIntersects_count, libspatialindex_c), RTError, (IndexH, Ptr{Cdouble}, Ptr{Cdouble}, Cdouble, Cdouble, UInt32, Ptr{UInt64}), index, pdMin, pdMax, tStart, tEnd, nDimension, nResults)
end

function Index_Intersects_count(index::IndexH, pdMin, pdMax, nDimension::Integer, nResults)
    ccall((:Index_Intersects_count, libspatialindex_c), RTError, (IndexH, Ptr{Cdouble}, Ptr{Cdouble}, UInt32, Ptr{UInt64}), index, pdMin, pdMax, nDimension, nResults)
end

function Index_TPNearestNeighbors_obj(index::IndexH, pdMin, pdMax, pdVMin, pdVMax, tStart::Cdouble, tEnd::Cdouble, nDimension::Integer, items, nResults)
    ccall((:Index_TPNearestNeighbors_obj, libspatialindex_c), RTError, (IndexH, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Cdouble, Cdouble, UInt32, Ptr{Ptr{IndexItemH}}, Ptr{UInt64}), index, pdMin, pdMax, pdVMin, pdVMax, tStart, tEnd, nDimension, items, nResults)
end

function Index_MVRNearestNeighbors_obj(index::IndexH, pdMin, pdMax, tStart::Cdouble, tEnd::Cdouble, nDimension::Integer, items, nResults)
    ccall((:Index_MVRNearestNeighbors_obj, libspatialindex_c), RTError, (IndexH, Ptr{Cdouble}, Ptr{Cdouble}, Cdouble, Cdouble, UInt32, Ptr{Ptr{IndexItemH}}, Ptr{UInt64}), index, pdMin, pdMax, tStart, tEnd, nDimension, items, nResults)
end

function Index_NearestNeighbors_obj(index::IndexH, pdMin, pdMax, nDimension::Integer, items, nResults)
    ccall((:Index_NearestNeighbors_obj, libspatialindex_c), RTError, (IndexH, Ptr{Cdouble}, Ptr{Cdouble}, UInt32, Ptr{Ptr{IndexItemH}}, Ptr{UInt64}), index, pdMin, pdMax, nDimension, items, nResults)
end

function Index_TPNearestNeighbors_id(index::IndexH, pdMin, pdMax, pdVMin, pdVMax, tStart::Cdouble, tEnd::Cdouble, nDimension::Integer, ids, nResults)
    ccall((:Index_TPNearestNeighbors_id, libspatialindex_c), RTError, (IndexH, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Cdouble, Cdouble, UInt32, Ptr{Ptr{Int64}}, Ptr{UInt64}), index, pdMin, pdMax, pdVMin, pdVMax, tStart, tEnd, nDimension, ids, nResults)
end

function Index_MVRNearestNeighbors_id(index::IndexH, pdMin, pdMax, tStart::Cdouble, tEnd::Cdouble, nDimension::Integer, ids, nResults)
    ccall((:Index_MVRNearestNeighbors_id, libspatialindex_c), RTError, (IndexH, Ptr{Cdouble}, Ptr{Cdouble}, Cdouble, Cdouble, UInt32, Ptr{Ptr{Int64}}, Ptr{UInt64}), index, pdMin, pdMax, tStart, tEnd, nDimension, ids, nResults)
end

function Index_NearestNeighbors_id(index::IndexH, pdMin, pdMax, nDimension::Integer, items, nResults)
    ccall((:Index_NearestNeighbors_id, libspatialindex_c), RTError, (IndexH, Ptr{Cdouble}, Ptr{Cdouble}, UInt32, Ptr{Ptr{Int64}}, Ptr{UInt64}), index, pdMin, pdMax, nDimension, items, nResults)
end

function Index_GetBounds(index::IndexH, ppdMin, ppdMax, nDimension)
    ccall((:Index_GetBounds, libspatialindex_c), RTError, (IndexH, Ptr{Ptr{Cdouble}}, Ptr{Ptr{Cdouble}}, Ptr{UInt32}), index, ppdMin, ppdMax, nDimension)
end

function Index_GetLeaves(index::IndexH, nLeafNodes, nLeafSizes, nLeafIDs, nLeafChildIDs, pppdMin, pppdMax, nDimension)
    ccall((:Index_GetLeaves, libspatialindex_c), RTError, (IndexH, Ptr{UInt32}, Ptr{Ptr{UInt32}}, Ptr{Ptr{Int64}}, Ptr{Ptr{Ptr{Int64}}}, Ptr{Ptr{Ptr{Cdouble}}}, Ptr{Ptr{Ptr{Cdouble}}}, Ptr{UInt32}), index, nLeafNodes, nLeafSizes, nLeafIDs, nLeafChildIDs, pppdMin, pppdMax, nDimension)
end

function Index_SetResultSetOffset(index::IndexH, value::Integer)
    ccall((:Index_SetResultSetOffset, libspatialindex_c), RTError, (IndexH, Int64), index, value)
end

function Index_GetResultSetOffset(index::IndexH)
    ccall((:Index_GetResultSetOffset, libspatialindex_c), Int64, (IndexH,), index)
end

function Index_SetResultSetLimit(index::IndexH, value::Integer)
    ccall((:Index_SetResultSetLimit, libspatialindex_c), RTError, (IndexH, Int64), index, value)
end

function Index_GetResultSetLimit(index::IndexH)
    ccall((:Index_GetResultSetLimit, libspatialindex_c), Int64, (IndexH,), index)
end

function Index_DestroyObjResults(results, nResults::Integer)
    ccall((:Index_DestroyObjResults, libspatialindex_c), Nothing, (Ptr{IndexItemH}, UInt32), results, nResults)
end

function Index_ClearBuffer(index::IndexH)
    ccall((:Index_ClearBuffer, libspatialindex_c), Nothing, (IndexH,), index)
end

function Index_Free(object)
    ccall((:Index_Free, libspatialindex_c), Nothing, (Ptr{Nothing},), object)
end

function Index_Flush(index::IndexH)
    ccall((:Index_Flush, libspatialindex_c), Nothing, (IndexH,), index)
end

function IndexItem_Destroy(item::IndexItemH)
    ccall((:IndexItem_Destroy, libspatialindex_c), Nothing, (IndexItemH,), item)
end

function IndexItem_GetID(item::IndexItemH)
    ccall((:IndexItem_GetID, libspatialindex_c), Int64, (IndexItemH,), item)
end

function IndexItem_GetData(item::IndexItemH, data, length)
    ccall((:IndexItem_GetData, libspatialindex_c), RTError, (IndexItemH, Ptr{Ptr{UInt8}}, Ptr{UInt64}), item, data, length)
end

function IndexItem_GetBounds(item::IndexItemH, ppdMin, ppdMax, nDimension)
    ccall((:IndexItem_GetBounds, libspatialindex_c), RTError, (IndexItemH, Ptr{Ptr{Cdouble}}, Ptr{Ptr{Cdouble}}, Ptr{UInt32}), item, ppdMin, ppdMax, nDimension)
end

function IndexProperty_Create()
    ccall((:IndexProperty_Create, libspatialindex_c), IndexPropertyH, ())
end

function IndexProperty_Destroy(hProp::IndexPropertyH)
    ccall((:IndexProperty_Destroy, libspatialindex_c), Nothing, (IndexPropertyH,), hProp)
end

function IndexProperty_SetIndexType(iprop::IndexPropertyH, value::RTIndexType)
    ccall((:IndexProperty_SetIndexType, libspatialindex_c), RTError, (IndexPropertyH, RTIndexType), iprop, value)
end

function IndexProperty_GetIndexType(iprop::IndexPropertyH)
    ccall((:IndexProperty_GetIndexType, libspatialindex_c), RTIndexType, (IndexPropertyH,), iprop)
end

function IndexProperty_SetDimension(iprop::IndexPropertyH, value::Integer)
    ccall((:IndexProperty_SetDimension, libspatialindex_c), RTError, (IndexPropertyH, UInt32), iprop, value)
end

function IndexProperty_GetDimension(iprop::IndexPropertyH)
    ccall((:IndexProperty_GetDimension, libspatialindex_c), UInt32, (IndexPropertyH,), iprop)
end

function IndexProperty_SetIndexVariant(iprop::IndexPropertyH, value::RTIndexVariant)
    ccall((:IndexProperty_SetIndexVariant, libspatialindex_c), RTError, (IndexPropertyH, RTIndexVariant), iprop, value)
end

function IndexProperty_GetIndexVariant(iprop::IndexPropertyH)
    ccall((:IndexProperty_GetIndexVariant, libspatialindex_c), RTIndexVariant, (IndexPropertyH,), iprop)
end

function IndexProperty_SetIndexStorage(iprop::IndexPropertyH, value::RTStorageType)
    ccall((:IndexProperty_SetIndexStorage, libspatialindex_c), RTError, (IndexPropertyH, RTStorageType), iprop, value)
end

function IndexProperty_GetIndexStorage(iprop::IndexPropertyH)
    ccall((:IndexProperty_GetIndexStorage, libspatialindex_c), RTStorageType, (IndexPropertyH,), iprop)
end

function IndexProperty_SetPagesize(iprop::IndexPropertyH, value::Integer)
    ccall((:IndexProperty_SetPagesize, libspatialindex_c), RTError, (IndexPropertyH, UInt32), iprop, value)
end

function IndexProperty_GetPagesize(iprop::IndexPropertyH)
    ccall((:IndexProperty_GetPagesize, libspatialindex_c), UInt32, (IndexPropertyH,), iprop)
end

function IndexProperty_SetIndexCapacity(iprop::IndexPropertyH, value::Integer)
    ccall((:IndexProperty_SetIndexCapacity, libspatialindex_c), RTError, (IndexPropertyH, UInt32), iprop, value)
end

function IndexProperty_GetIndexCapacity(iprop::IndexPropertyH)
    ccall((:IndexProperty_GetIndexCapacity, libspatialindex_c), UInt32, (IndexPropertyH,), iprop)
end

function IndexProperty_SetLeafCapacity(iprop::IndexPropertyH, value::Integer)
    ccall((:IndexProperty_SetLeafCapacity, libspatialindex_c), RTError, (IndexPropertyH, UInt32), iprop, value)
end

function IndexProperty_GetLeafCapacity(iprop::IndexPropertyH)
    ccall((:IndexProperty_GetLeafCapacity, libspatialindex_c), UInt32, (IndexPropertyH,), iprop)
end

function IndexProperty_SetLeafPoolCapacity(iprop::IndexPropertyH, value::Integer)
    ccall((:IndexProperty_SetLeafPoolCapacity, libspatialindex_c), RTError, (IndexPropertyH, UInt32), iprop, value)
end

function IndexProperty_GetLeafPoolCapacity(iprop::IndexPropertyH)
    ccall((:IndexProperty_GetLeafPoolCapacity, libspatialindex_c), UInt32, (IndexPropertyH,), iprop)
end

function IndexProperty_SetIndexPoolCapacity(iprop::IndexPropertyH, value::Integer)
    ccall((:IndexProperty_SetIndexPoolCapacity, libspatialindex_c), RTError, (IndexPropertyH, UInt32), iprop, value)
end

function IndexProperty_GetIndexPoolCapacity(iprop::IndexPropertyH)
    ccall((:IndexProperty_GetIndexPoolCapacity, libspatialindex_c), UInt32, (IndexPropertyH,), iprop)
end

function IndexProperty_SetRegionPoolCapacity(iprop::IndexPropertyH, value::Integer)
    ccall((:IndexProperty_SetRegionPoolCapacity, libspatialindex_c), RTError, (IndexPropertyH, UInt32), iprop, value)
end

function IndexProperty_GetRegionPoolCapacity(iprop::IndexPropertyH)
    ccall((:IndexProperty_GetRegionPoolCapacity, libspatialindex_c), UInt32, (IndexPropertyH,), iprop)
end

function IndexProperty_SetPointPoolCapacity(iprop::IndexPropertyH, value::Integer)
    ccall((:IndexProperty_SetPointPoolCapacity, libspatialindex_c), RTError, (IndexPropertyH, UInt32), iprop, value)
end

function IndexProperty_GetPointPoolCapacity(iprop::IndexPropertyH)
    ccall((:IndexProperty_GetPointPoolCapacity, libspatialindex_c), UInt32, (IndexPropertyH,), iprop)
end

function IndexProperty_SetBufferingCapacity(iprop::IndexPropertyH, value::Integer)
    ccall((:IndexProperty_SetBufferingCapacity, libspatialindex_c), RTError, (IndexPropertyH, UInt32), iprop, value)
end

function IndexProperty_GetBufferingCapacity(iprop::IndexPropertyH)
    ccall((:IndexProperty_GetBufferingCapacity, libspatialindex_c), UInt32, (IndexPropertyH,), iprop)
end

function IndexProperty_SetEnsureTightMBRs(iprop::IndexPropertyH, value::Integer)
    ccall((:IndexProperty_SetEnsureTightMBRs, libspatialindex_c), RTError, (IndexPropertyH, UInt32), iprop, value)
end

function IndexProperty_GetEnsureTightMBRs(iprop::IndexPropertyH)
    ccall((:IndexProperty_GetEnsureTightMBRs, libspatialindex_c), UInt32, (IndexPropertyH,), iprop)
end

function IndexProperty_SetOverwrite(iprop::IndexPropertyH, value::Integer)
    ccall((:IndexProperty_SetOverwrite, libspatialindex_c), RTError, (IndexPropertyH, UInt32), iprop, value)
end

function IndexProperty_GetOverwrite(iprop::IndexPropertyH)
    ccall((:IndexProperty_GetOverwrite, libspatialindex_c), UInt32, (IndexPropertyH,), iprop)
end

function IndexProperty_SetNearMinimumOverlapFactor(iprop::IndexPropertyH, value::Integer)
    ccall((:IndexProperty_SetNearMinimumOverlapFactor, libspatialindex_c), RTError, (IndexPropertyH, UInt32), iprop, value)
end

function IndexProperty_GetNearMinimumOverlapFactor(iprop::IndexPropertyH)
    ccall((:IndexProperty_GetNearMinimumOverlapFactor, libspatialindex_c), UInt32, (IndexPropertyH,), iprop)
end

function IndexProperty_SetWriteThrough(iprop::IndexPropertyH, value::Integer)
    ccall((:IndexProperty_SetWriteThrough, libspatialindex_c), RTError, (IndexPropertyH, UInt32), iprop, value)
end

function IndexProperty_GetWriteThrough(iprop::IndexPropertyH)
    ccall((:IndexProperty_GetWriteThrough, libspatialindex_c), UInt32, (IndexPropertyH,), iprop)
end

function IndexProperty_SetFillFactor(iprop::IndexPropertyH, value::Cdouble)
    ccall((:IndexProperty_SetFillFactor, libspatialindex_c), RTError, (IndexPropertyH, Cdouble), iprop, value)
end

function IndexProperty_GetFillFactor(iprop::IndexPropertyH)
    ccall((:IndexProperty_GetFillFactor, libspatialindex_c), Cdouble, (IndexPropertyH,), iprop)
end

function IndexProperty_SetSplitDistributionFactor(iprop::IndexPropertyH, value::Cdouble)
    ccall((:IndexProperty_SetSplitDistributionFactor, libspatialindex_c), RTError, (IndexPropertyH, Cdouble), iprop, value)
end

function IndexProperty_GetSplitDistributionFactor(iprop::IndexPropertyH)
    ccall((:IndexProperty_GetSplitDistributionFactor, libspatialindex_c), Cdouble, (IndexPropertyH,), iprop)
end

function IndexProperty_SetTPRHorizon(iprop::IndexPropertyH, value::Cdouble)
    ccall((:IndexProperty_SetTPRHorizon, libspatialindex_c), RTError, (IndexPropertyH, Cdouble), iprop, value)
end

function IndexProperty_GetTPRHorizon(iprop::IndexPropertyH)
    ccall((:IndexProperty_GetTPRHorizon, libspatialindex_c), Cdouble, (IndexPropertyH,), iprop)
end

function IndexProperty_SetReinsertFactor(iprop::IndexPropertyH, value::Cdouble)
    ccall((:IndexProperty_SetReinsertFactor, libspatialindex_c), RTError, (IndexPropertyH, Cdouble), iprop, value)
end

function IndexProperty_GetReinsertFactor(iprop::IndexPropertyH)
    ccall((:IndexProperty_GetReinsertFactor, libspatialindex_c), Cdouble, (IndexPropertyH,), iprop)
end

function IndexProperty_SetFileName(iprop::IndexPropertyH, value)
    ccall((:IndexProperty_SetFileName, libspatialindex_c), RTError, (IndexPropertyH, Cstring), iprop, value)
end

function IndexProperty_GetFileName(iprop::IndexPropertyH)
    ccall((:IndexProperty_GetFileName, libspatialindex_c), Cstring, (IndexPropertyH,), iprop)
end

function IndexProperty_SetFileNameExtensionDat(iprop::IndexPropertyH, value)
    ccall((:IndexProperty_SetFileNameExtensionDat, libspatialindex_c), RTError, (IndexPropertyH, Cstring), iprop, value)
end

function IndexProperty_GetFileNameExtensionDat(iprop::IndexPropertyH)
    ccall((:IndexProperty_GetFileNameExtensionDat, libspatialindex_c), Cstring, (IndexPropertyH,), iprop)
end

function IndexProperty_SetFileNameExtensionIdx(iprop::IndexPropertyH, value)
    ccall((:IndexProperty_SetFileNameExtensionIdx, libspatialindex_c), RTError, (IndexPropertyH, Cstring), iprop, value)
end

function IndexProperty_GetFileNameExtensionIdx(iprop::IndexPropertyH)
    ccall((:IndexProperty_GetFileNameExtensionIdx, libspatialindex_c), Cstring, (IndexPropertyH,), iprop)
end

function IndexProperty_SetCustomStorageCallbacksSize(iprop::IndexPropertyH, value::Integer)
    ccall((:IndexProperty_SetCustomStorageCallbacksSize, libspatialindex_c), RTError, (IndexPropertyH, UInt32), iprop, value)
end

function IndexProperty_GetCustomStorageCallbacksSize(iprop::IndexPropertyH)
    ccall((:IndexProperty_GetCustomStorageCallbacksSize, libspatialindex_c), UInt32, (IndexPropertyH,), iprop)
end

function IndexProperty_SetCustomStorageCallbacks(iprop::IndexPropertyH, value)
    ccall((:IndexProperty_SetCustomStorageCallbacks, libspatialindex_c), RTError, (IndexPropertyH, Ptr{Nothing}), iprop, value)
end

function IndexProperty_GetCustomStorageCallbacks(iprop::IndexPropertyH)
    ccall((:IndexProperty_GetCustomStorageCallbacks, libspatialindex_c), Ptr{Nothing}, (IndexPropertyH,), iprop)
end

function IndexProperty_SetIndexID(iprop::IndexPropertyH, value::Integer)
    ccall((:IndexProperty_SetIndexID, libspatialindex_c), RTError, (IndexPropertyH, Int64), iprop, value)
end

function IndexProperty_GetIndexID(iprop::IndexPropertyH)
    ccall((:IndexProperty_GetIndexID, libspatialindex_c), Int64, (IndexPropertyH,), iprop)
end

function SIDX_NewBuffer(bytes::Integer)
    ccall((:SIDX_NewBuffer, libspatialindex_c), Ptr{Nothing}, (Cint,), bytes)
end

function SIDX_DeleteBuffer(buffer)
    ccall((:SIDX_DeleteBuffer, libspatialindex_c), Nothing, (Ptr{Nothing},), buffer)
end

function IndexProperty_SetResultSetLimit(iprop::IndexPropertyH, value::UInt64)
    ccall((:IndexProperty_SetResultSetLimit, libspatialindex_c), RTError, (IndexPropertyH, UInt64), iprop, value)
end

function IndexProperty_GetResultSetLimit(iprop::IndexPropertyH)
    ccall((:IndexProperty_GetResultSetLimit, libspatialindex_c), UInt64, (IndexPropertyH,), iprop)
end

function SIDX_Version()
    ccall((:SIDX_Version, libspatialindex_c), Cstring, ())
end

function Error_GetLastErrorMsg()
    ccall((:Error_GetLastErrorMsg, libspatialindex_c), Cstring, ())
end


text:
using HDF5, JLD
using PyPlot

function get_wf(filename, dset_real="psi_real", dset_imag="psi_imag")
    fid = h5open(filename, "r")
    xmin = read(attrs(fid)["xmin"])
    xmax = read(attrs(fid)["xmax"])
    nx = read(attrs(fid)["nx"])
    ix0 = read(attrs(fid)["ix0"])
    sx = size(fid[dset_real])[1]
    dx = (xmax-xmin)/nx
    x_grid = linspace(xmin+dx*ix0,xmin+dx*(ix0+sx-1),sx)
    if exists(fid,dset_imag)
        psi = read(fid[dset_real])+1im* read(fid[dset_imag])
    else
        psi = read(fid[dset_real])
    end
    close(fid)
    return x_grid, psi
end

function plot_wf(filename)
    x_grid, psi = get_wf(filename)
    plot(x_grid, abs(psi).^2)
end

function plot_potential(filename)
    x_grid, pot = get_wf(filename, "potential", "dummy_imag")
    plot(x_grid, pot)
end


text:
#================================================================

Definition of SimHash, an LSH function for hashing on cosine similarity.

================================================================#

#========================
Typedefs
========================#

mutable struct SimHash{T <: Union{Float32,Float64}} <: SymmetricLSHFunction
    # Random coefficients sampled for SimHash. Each column of coefficients
    # corresponds to another hash function.
    coeff :: Matrix{T}

    # Whether or not SimHash should round up to the next power of 2 when
    # resizing its coefficient array.
    #
    # This is useful if you're going to be hashing inputs of many different
    # sizes, since every call to resize! entails allocating a new array and
    # copying the old coefficient array into it. However, if you're only
    # going to be hashing inputs of one size or a few different sizes,
    # then you should set resize_pow2 to false.
    resize_pow2 :: Bool
end

function SimHash{T}(n_hashes::Integer = DEFAULT_N_HASHES;
                    resize_pow2::Bool = DEFAULT_RESIZE_POW2) where T

    coeff = Matrix{T}(undef, 0, n_hashes)
    SimHash{T}(coeff, resize_pow2)
end

SimHash(args...; dtype = DEFAULT_DTYPE, kws...) =
    SimHash{dtype}(args...; kws...)

@doc """
    SimHash(n_hashes::Integer = $(DEFAULT_N_HASHES);
            dtype::Type = $(DEFAULT_DTYPE),
            resize_pow2::Bool = $(DEFAULT_RESIZE_POW2))

Creates a locality-sensitive hash function for cosine similarity.

# Arguments
- $(N_HASHES_DOCSTR())

# Keyword parameters
- $(DTYPE_DOCSTR(SimHash))
- $(RESIZE_POW2_DOCSTR(SimHash))

# Examples
Construct a hash function by calling `SimHash` with the number of hash functions you want to generate:

```jldoctest; setup = :(using LSHFunctions)
julia> hashfn = SimHash(24);

julia> n_hashes(hashfn) == 24 &&
       similarity(hashfn) == cossim
true
```

You can then call `hashfn(x)` in order to compute hashes:

```jldoctest; setup = :(using LSHFunctions)
julia> hashfn = SimHash(32);

julia> x = randn(30);

julia> hashes = hashfn(x);

```

# References

- Moses S. Charikar. *Similarity estimation techniques from rounding algorithms*. In Proceedings of the Thiry-Fourth Annual ACM Symposium on Theory of Computing, STOC '02, page 380–388, New York, NY, USA, 2002. Association for Computing Machinery. 10.1145/509907.509965.

See also: [`cossim`](@ref)
""" SimHash

#========================
SimHash helper functions
========================#

# Resize SimHash to hash vectors of length up to n
function Base.resize!(hashfn::SimHash, n :: Integer)
    # If resize_pow2 is set, then change n to the next power of 2 that
    # is larger than n.
    n = (hashfn.resize_pow2) ? nextpow(2, n) : n

    old_n, n_hashes = size(hashfn.coeff)
    old_coeff = hashfn.coeff
    new_coeff = similar(hashfn.coeff, n, n_hashes)

    # Copy old coefficients into the new_coeff Matrix
    @views new_coeff[1:min(n,old_n),1:end] .= old_coeff

    # Generate new coefficients in the empty entries
    if n > old_n
        @views map!(_ -> randn(),
                    new_coeff[old_n+1:end,1:end],
                    new_coeff[old_n+1:end,1:end])
    end

    hashfn.coeff = new_coeff
end

# The current maximum input size that a SimHash struct can accept based on
# the number of coefficients that we've generated for it. current_ is
# prepended because this implementation of SimHash supports vectors of
# different sizes.
current_max_input_size(hashfn::SimHash) = size(hashfn.coeff, 1)

#========================
LSHFunction and SymmetricLSHFunction API compliance
========================#

hashtype(::SimHash) = Bool
n_hashes(hashfn::SimHash) = size(hashfn.coeff, 2)
similarity(::SimHash) = cossim
single_hash_collision_probability(::SimHash, sim::Real) =
    @. (1 - acos(sim) / π)

### Hash computation

# Perform type conversion to hit BLAS when the input array has a different
# element type than the hashfn coefficients.
(hashfn::SimHash{T})(x::AbstractArray) where T =
	hashfn(T.(x))

function (hashfn::SimHash{T})(x::AbstractArray{T}) where T
    n = size(x,1)
    if n > current_max_input_size(hashfn)
        resize!(hashfn, n)
    end

    coeff = @views hashfn.coeff[1:n,1:end]

    (coeff' * x) .≥ 0
end


text:
# low-level wrappers of the MAGMA library
include("common.jl")

# magma_init
function magmaInit()
	ccall((:magma_init, libmagma),Cint,())
end

# magma_finalize
function magmaFinalize()
	ccall((:magma_finalize, libmagma),Cint,())
end

# magma_cgesvd wrapper
function magmaCgesvd(jobu, jobvt, m, n, A, lda, s, U,
	ldu, VT, ldvt, work, lwork, rwork, info)
    ccall((:magma_cgesvd, libmagma),
               Cint,
               (Cint, Cint, Cint, Cint, PtrOrCuPtr{ComplexF32},
			    Cint, PtrOrCuPtr{Float32}, PtrOrCuPtr{ComplexF32},
				Cint, PtrOrCuPtr{ComplexF32}, Cint,
				PtrOrCuPtr{ComplexF32}, Cint, PtrOrCuPtr{Float32},
				PtrOrCuPtr{Cint}),
                jobu, jobvt, m, n, A, lda, s, U, ldu, VT, ldvt,
			    work, lwork, rwork, info)
end

# magma_dgesvd wrapper
function magmaDgesvd(jobu, jobvt, m, n, A, lda, s, U,
	 ldu, VT, ldvt, work, lwork, info)
    ccall((:magma_dgesvd, libmagma),
               Cint,
               (Cint, Cint, Cint, Cint, PtrOrCuPtr{Float64},
			    Cint, PtrOrCuPtr{Float64}, PtrOrCuPtr{Float64},
				Cint, PtrOrCuPtr{Float64}, Cint,
				PtrOrCuPtr{Float64}, Cint, PtrOrCuPtr{Cint}),
                jobu, jobvt, m, n, A, lda, s, U, ldu, VT, ldvt,
			    work, lwork, info)
end

# magma_sgesvd wrapper
function magmaSgesvd(jobu, jobvt, m, n, A, lda, s, U,
	 ldu, VT, ldvt, work, lwork, info)
    ccall((:magma_sgesvd, libmagma),
               Cint,
               (Cint, Cint, Cint, Cint, PtrOrCuPtr{Float32},
			    Cint, PtrOrCuPtr{Float32}, PtrOrCuPtr{Float32},
				Cint, PtrOrCuPtr{Float32}, Cint, PtrOrCuPtr{Float32},
				Cint, Ptr{Cint}),
                jobu, jobvt, m, n, A, lda, s, U, ldu, VT, ldvt,
				work, lwork, info)
end

# magma_cgesvd wrapper
function magmaZgesvd(jobu, jobvt, m, n, A, lda, s, U,
	 ldu, VT, ldvt, work, lwork, rwork, info)
    ccall((:magma_zgesvd, libmagma),
               Cint,
               (Cint, Cint, Cint, Cint, PtrOrCuPtr{ComplexF64},
			    Cint, PtrOrCuPtr{Float64}, PtrOrCuPtr{ComplexF64},
				Cint, PtrOrCuPtr{ComplexF64}, Cint, PtrOrCuPtr{ComplexF64},
				Cint, PtrOrCuPtr{Float64}, PtrOrCuPtr{Cint}),
                jobu, jobvt, m, n, A, lda, s, U, ldu, VT, ldvt,
			    work, lwork, rwork, info)
end

# wrappers of the low-level MAGMA functions
for (function_name, element_type, singular_value_type) in
	((:magmaCgesvd,:ComplexF32,:Float32),
                      (:magmaZgesvd,:ComplexF32,:Float64))
	@eval begin
		# magma_int_t magma_cgesvd 	( 	magma_vec_t  	jobu,
		# 	magma_vec_t  	jobvt,
		# 	magma_int_t  	m,
		# 	magma_int_t  	n,
		# 	magmaFloatComplex *  	A,
		# 	magma_int_t  	lda,
		# 	float *  	s,
		# 	magmaFloatComplex *  	U,
		# 	magma_int_t  	ldu,
		# 	magmaFloatComplex *  	VT,
		# 	magma_int_t  	ldvt,
		# 	magmaFloatComplex *  	work,
		# 	magma_int_t  	lwork,
		# 	float *  	rwork,
		# 	magma_int_t *  	info
		# )
		function gesvd!(jobu::Cint,
                      jobvt::Cint,
                      A::CuMatrix{$element_type},
                      ldu::Cint,
                      ldvt::Cint,
					  lwork::Cint,
					  rwork::PtrOrCuPtr{$singular_value_type})
			convert(Cint,jobu)
			convert(Cint,jobvt)
           m, n = size(A)
           # there should be some certain checking
		   # for the input matrix and arrays
           # before we can calculate further input
		   # and call the lower C function wrappers
           lda = max(1,stride(A,2))
           s = CuArray{$singular_value_type}(0,min(m,n))

		   # there should be some conditions for the size of U and VT
		   # but we will deal with them later
		   U = zeros($element_type,ldu,m)
		   VT = zeros($element_type,ldvt,n)

		   work = zeros($element_type,max(1,lwork))
		   info = zeros(Cint,1)

           $function_name(jobu, jobvt, m, n, A, lda, s, U, ldu, VT,
		    ldvt, work, lwork, rwork, info)
           return U, s, VT, work, info
	   end
    end
end
for (function_name, element_type, singular_value_type) in
	 ((:magmaDgesvd,:Float64,:Float64),
                      (:magmaSgesvd,:Float32,:Float32))
	@eval begin
		# magma_int_t magma_dgesvd 	( 	magma_vec_t  	jobu,
		# 	magma_vec_t  	jobvt,
		# 	magma_int_t  	m,
		# 	magma_int_t  	n,
		# 	magmaFloatComplex *  	A,
		# 	magma_int_t  	lda,
		# 	float *  	s,
		# 	magmaFloatComplex *  	U,
		# 	magma_int_t  	ldu,
		# 	magmaFloatComplex *  	VT,
		# 	magma_int_t  	ldvt,
		# 	magmaFloatComplex *  	work,
		# 	magma_int_t  	lwork,
		# 	magma_int_t *  	info
		# )
		Base.unsafe_convert(::CuPtr{$element_type},x::CuArray{$element_type})=Base.unsafe_convert(CuPtr{$element_type},Base.cconvert(CuPtr{($element_type)}x))
		Base.unsafe_convert(::PtrOrCuPtr{$element_type},x::CuArray{$element_type})=Base.unsafe_convert(PtrOrCuPtr{$element_type},Base.cconvert(PtrOrCuPtr{($element_type)}x))
		function gesvd!(jobu,
                      jobvt,
                      A::CuMatrix{$element_type},
                      ldu,
                      ldvt,
					  lwork)
			convert(Cint,jobu)
			convert(Cint,jobvt)
           m, n = size(A)
           # there should be some certain checking
		   # for the input matrix and arrays
           # before we can calculate further input
		   # and call the lower C function wrappers
           lda = max(1,stride(A,2))
           s = cu(zeros(m,n))

		   # there should be some conditions for the size of U and VT
		   # but we will deal with them later
		   U = cu(zeros($element_type,ldu,m))
		   VT = cu(zeros($element_type,ldvt,n))

		   work = cu(zeros($element_type,max(1,lwork)))
		   info = (zeros(Cint,1))

		   $function_name(jobu, jobvt, m, n, A, lda, s, U, ldu, VT,
		    ldvt, work, lwork, info)
           return U, s, VT, work, info
	   end
    end
end

#print(1+1)
print("Test begins!\n")

jobu = MagmaAllVec
jobvt = MagmaAllVec
A = cu(zeros(2,2))
A[1,1]=2.0
A[1,2]=1.0
A[2,1]=1.0
A[2,2]=2.0
ldu=2
ldvt=2
lwork=134
success=magmaInit()

print("Magma Initialization success=", success)
print('\n')

U, s, VT, work, info = gesvd!(jobu,jobvt,A,ldu,ldvt,lwork)

print("U=")
print(U)
print('\n')

print("s=")
print(s)
print('\n')

print("V**T=")
print(VT)
print('\n')

print("work=")
print(work)
print('\n')

print("info=")
print(info)
print('\n')

print("Test ends.")


text:
function extract_orders(data::AbstractVector{ST}, pipeline::PipelinePlan; recalc::Bool = false, verbose::Bool = false,
      remove_bad_chunks::Bool = true, orders_to_use = RvSpectMLBase.orders_to_use_default(first(data).inst) ) where { ST<:AbstractSpectra }
   if need_to(pipeline,:extract_orders) || recalc
      if verbose println("# Extracting order list timeseries from spectra.")  end
      @assert !need_to(pipeline,:read_spectra)
      order_list_timeseries = make_order_list_timeseries(data, orders_to_use=orders_to_use)
      if remove_bad_chunks
        order_list_timeseries = filter_bad_chunks(order_list_timeseries,verbose=verbose)
      end
      #RvSpectML.normalize_spectra!(order_list_timeseries,data);
      set_cache!(pipeline,:extract_orders,order_list_timeseries)
      dont_need_to!(pipeline,:extract_orders)
    end
    if has_cache(pipeline,:extract_orders) return read_cache(pipeline,:extract_orders)
    else   @error("Invalid pipeline state.")          end
end


text:
using PlanktonIndividuals, Serialization

grid = RegularRectilinearGrid(size = (1, 1, 1), spacing = (32, 32, 32), halo = (2, 2, 2))

phyt_params = deserialize(dirname(pathof(PlanktonIndividuals))*"/../test/param5.bin")

model = PlanktonModel(CPU(), grid;
                      N_species = 5,
                      N_individual = 1024,
                      max_individuals = 1024*10,
                      phyt_params = update_phyt_params(phyt_params, QuotaMode()))

TP = sum((interior(model.nutrients.PO4.data, grid) .+
          interior(model.nutrients.DOP.data, grid) .+
          interior(model.nutrients.POP.data, grid)) .* grid.Δx .* grid.Δy .* grid.Δz)
TP = TP + sum(model.individuals.phytos.sp1.data.Pq .+
              model.individuals.phytos.sp1.data.Bm .* model.individuals.phytos.sp1.p.R_PC) +
          sum(model.individuals.phytos.sp2.data.Pq .+
              model.individuals.phytos.sp2.data.Bm .* model.individuals.phytos.sp2.p.R_PC) +
          sum(model.individuals.phytos.sp3.data.Pq .+
              model.individuals.phytos.sp3.data.Bm .* model.individuals.phytos.sp3.p.R_PC) +
          sum(model.individuals.phytos.sp4.data.Pq .+
              model.individuals.phytos.sp4.data.Bm .* model.individuals.phytos.sp4.p.R_PC) +
          sum(model.individuals.phytos.sp5.data.Pq .+
              model.individuals.phytos.sp5.data.Bm .* model.individuals.phytos.sp5.p.R_PC)

sim = PlanktonSimulation(model, ΔT = 60, iterations = 10)

update!(sim)

TPt = sum((interior(model.nutrients.PO4.data, grid) .+
          interior(model.nutrients.DOP.data, grid) .+
          interior(model.nutrients.POP.data, grid)) .* grid.Δx .* grid.Δy .* grid.Δz)
TPt=TPt + sum(model.individuals.phytos.sp1.data.Pq .+
                model.individuals.phytos.sp1.data.Bm .* model.individuals.phytos.sp1.p.R_PC) +
          sum(model.individuals.phytos.sp2.data.Pq .+
              model.individuals.phytos.sp2.data.Bm .* model.individuals.phytos.sp2.p.R_PC) +
          sum(model.individuals.phytos.sp3.data.Pq .+
              model.individuals.phytos.sp3.data.Bm .* model.individuals.phytos.sp3.p.R_PC) +
          sum(model.individuals.phytos.sp4.data.Pq .+
              model.individuals.phytos.sp4.data.Bm .* model.individuals.phytos.sp4.p.R_PC) +
          sum(model.individuals.phytos.sp5.data.Pq .+
              model.individuals.phytos.sp5.data.Bm .* model.individuals.phytos.sp5.p.R_PC)


text:
clone_dir(repo) = joinpath(homedir(), repo.dir)

clone!(repo::MockRepo) = nothing
clone!(repo::MockFileRepo) = nothing

function clean_repo!(dir)
    paths = readdir(dir; join=true)
    for path in paths
        if !endswith(path, ".git")
            rm(path; force=true, recursive=true)
        end
    end
    return nothing
end

function clone!(repo::GitHubRepo)
    user = string(repo.user)
    token = string(repo.token)
    repository = string(repo.repo)
    branch = string(repo.branch)
    dir = clone_dir(repo)
    url = if token == ""
        "https://github.com/$repository.git"
    else
        "https://$user:$token@github.com/$repository.git"
    end
    run(`git clone $url $dir`)
    cd(dir) do
        try
            run(`git checkout $branch`)
        catch
            run(`git checkout --orphan $branch`)
            clean_repo!(dir)
        end
    end
    return dir
end

pull_or_clone!(repo::MockFileRepo) = nothing
pull_or_clone!(repo::MockRepo) = nothing

function pull_or_clone!(repo::Repo)
    dir = repo.dir
    if isdir(dir) && ".git" in readdir(dir)
        cd(dir) do
            run(`git pull`)
        end
    else
        clone!(repo)
    end
end

commit!(repo::MockRepo) = nothing
commit!(repo::MockFileRepo) = nothing

git_unchanged() = read(`git status --porcelain`, String) == ""

function commit!(repo::Repo)
    dir = clone_dir(repo)
    branch = repo.branch
    cd(dir) do
        if "CI" in keys(ENV)
            run(`git config --global user.email "skansbot@example.com"`)
            run(`git config --global user.name "skansbot"`)
        end
        if git_unchanged()
            println("Nothing changed")
        else
            run(`git add .`)
            run(`git commit -m '[Bot] Update stored pages'`)
            # Required for first push on new branch.
            run(`git push --set-upstream origin $branch`)
        end
    end
    return nothing
end

"""
    diff(dir::AbstractString=pwd())

Return the output of `git diff` inside `dir`.
"""
function diff(dir::AbstractString=pwd())
    cd(dir) do
        cmd = `git diff`
        return read(cmd, String)
    end
end

"""
    startswith_one(s::AbstractString, c::Char)::Bool

Return `true` if and only if `c` is at the first position in `s` and not the second.
"""
function startswith_one(s::AbstractString, c::Char)::Bool
    if length(s) == 0
        return false
    elseif length(s) == 1
        return s[1] == c
    else
        return s[1] == c && s[2] != c
    end
end

struct BeforeAfter
    old::Union{PageScan,Nothing}
    new::PageScan
end

"""
    diff(old::AbstractString, new::AbstractString)

Return a diff comparing `old` to `new`.
"""
function diff(old::AbstractString, new::AbstractString)
    mktempdir() do dir
        old_path = joinpath(dir, "old.txt")
        new_path = joinpath(dir, "new.txt")
        write(old_path, old * '\n')
        write(new_path, new * '\n')
        cmd = `git diff --no-index $old_path $new_path`
        return read(ignorestatus(cmd), String)
    end
end

function diff(ba::BeforeAfter)
    old = isnothing(ba.old) ? "" : ba.old.content
    return diff(old, ba.new.content)
end

"""
    cleandiff(uncleaned::AbstractString)

Return the output of a cleaned up `git diff` inside `dir`.
This keeps only lines starting with `+` or `-` except `+++` or `---` lines.
"""
function cleandiff(uncleaned::AbstractString)
    sep = '\n'
    lines = split(uncleaned, sep)
    filtered = filter(lines) do line
        startswith_one(line, '+') || startswith_one(line, '-')
    end
    threshold = 22
    if threshold < length(filtered)
        filtered = last(filtered, threshold)
        pushfirst!(filtered, "[...]")
    end
    return join(filtered, sep)
end

function code_block(s::AbstractString, class::AbstractString)
    return """
        ```$class
        $s
        ```
        """
end


text:
using Base: Symbol
"""

    plot_posterior_intervals(model, args...; kwargs...)

Draw a plot with a point estimate measure of centrality and quantiled credible intervals for easy interpretation of regression models fit in `Turing.jl` and `Soss.jl`.

Usage:
```julia-repl
plot_posterior_intervals(model)
```

Arguments:

- `model` : The `Turing.jl` or `Soss.jl` model of class `Chains` or `Array` to draw inferences from.
"""
function plot_posterior_intervals(model, args...; kwargs...)

    # Fix seed for reproducibility

    Random.seed!(123)

    # Turn model objects into DataFrames

    if isa(model, Array)

        myarray = DataFrame(model)

        #-------- Wrangle floats --------

        # Reshape

        thefloats = select(myarray, findall(col -> eltype(col) <: Float64, eachcol(myarray)))
        ncols = size(thefloats, 2)
        stackedfloats = stack(thefloats, 1:ncols)

        # Median

        centrefloats = combine(groupby(stackedfloats, :variable), :value => median)
        centrefloats = rename(centrefloats, :value_median => :centre)

        # Lower quantile

        lowerfloats = combine(groupby(stackedfloats, :variable), :value => t -> quantile(t, .025))
        lowerfloats = rename(lowerfloats, :value_function => :lower)

        # Upper quantile

        upperfloats = combine(groupby(stackedfloats, :variable), :value => t -> quantile(t, .975))
        upperfloats = rename(upperfloats, :value_function => :upper)

        # Join together

        finalfloats = leftjoin(centrefloats, lowerfloats, on = :variable)
        finalfloats = leftjoin(finalfloats, upperfloats, on = :variable)

        #-------- Wrangle arrays --------

        # Reshape

        thearrays = select(myarray, findall(col -> eltype(col) <: Array, eachcol(myarray)))
        ncols_arrays = size(thearrays, 2)

        if ncols_arrays > 1
            error("`plot_posterior_intervals` currently only works for models where an array of parameters exists for one type of parameter (e.g. β coefficients stored in the array). Other parameters can be stored as regular Float64 columns.")
        else

            # Parse columns and retain only the separated ones

            parsedarrays = expandcol(thearrays, Symbol(names(thearrays)[1]))

            parsedarraysCols = select(parsedarrays, findall(col -> eltype(col) <: Float64, eachcol(parsedarrays)))

            # Reshape

            ncols_res_arrays = size(parsedarraysCols, 2)
            stackedarraysfloats = stack(parsedarraysCols, 1:ncols_res_arrays)

            # Median

            centrefloatsarrays = combine(groupby(stackedarraysfloats, :variable), :value => median)
            centrefloatsarrays = rename(centrefloatsarrays, :value_median => :centre)

            # Lower quantile

            lowerfloatsarrays = combine(groupby(stackedarraysfloats, :variable), :value => t -> quantile(t, .025))
            lowerfloatsarrays = rename(lowerfloatsarrays, :value_function => :lower)

            # Upper quantile

            upperfloatsarrays = combine(groupby(stackedarraysfloats, :variable), :value => t -> quantile(t, .975))
            upperfloatsarrays = rename(upperfloatsarrays, :value_function => :upper)

            # Join together

            finalfloatsarrays = leftjoin(centrefloatsarrays, lowerfloatsarrays, on = :variable)
            finalfloatsarrays = leftjoin(finalfloatsarrays, upperfloatsarrays, on = :variable)
        end

        #-------- Final outputs ---------

        # Merge together

        finalPost = [finalfloats; finalfloatsarrays]
        finalPost = rename(finalPost, :variable => :parameters)

        # Standardise outputs

        variable = finalPost[!, :parameters]
        variable = string.(variable)
        lower = finalPost[!, :lower]
        centre = finalPost[!, :centre]
        upper = finalPost[!, :upper]

    elseif isa(model, Chains)

        # Extract values

        finalPost = DataFrame(MCMCChains.quantile(model))

        finalPost = select(finalPost, "parameters" => "parameters", "2.5%" => "lower",
                            "50.0%" => "centre", "97.5%" => "upper")

        # Standardise outputs

        variable = finalPost[!, :parameters]
        variable = string.(variable)
        lower = finalPost[!, :lower]
        centre = finalPost[!, :centre]
        upper = finalPost[!, :upper]
    else
        error("`model` must be an object of type `Chains` or `Array`.")
    end

    #------------ Draw the plot -----------------

    gr() # gr backend for graphics

    myPlot = plot(centre, variable, xerror = (centre .- lower, upper .- centre), st = :scatter,
                title = "Posterior medians w/ 95% credible intervals",
                xlabel = "Value",
                ylabel = "Parameter",
                legend = false,
                markersize = 5,
                size = (600, 600))

    return myPlot
end

"""

    plot_posterior_hist(model, plot_legend, args...; kwargs...)

Draw a plot with a binned histogram of sampled parameters for easy interpretation of regression models fit in `Turing.jl` and `Soss.jl`.

Usage:
```julia-repl
plot_posterior_hist(model, plot_legend)
```

Details:

Note that to get the function to work, you may need to call it using a splat, such as `plot(plot_posterior_hist(model, plot_legend)...)`.

Arguments:

- `model` : The `Turing.jl` or `Soss.jl` model of class `Chains` or `Array` to draw inferences from.
- `plot_legend` : Boolean of whether to add a legend to the plot or not.
"""
function plot_posterior_hist(model, plot_legend::Bool, args...; kwargs...)

    if isa(model, Array)

        myarray = DataFrame(model)

        #-------- Wrangle floats --------

        thefloats = select(myarray, findall(col -> eltype(col) <: Float64, eachcol(myarray)))

        #-------- Wrangle arrays --------

        thearrays = select(myarray, findall(col -> eltype(col) <: Array, eachcol(myarray)))
        ncols_arrays = size(thearrays, 2)

        if ncols_arrays > 1
            error("`plot_posterior_intervals` currently only works for models where an array of parameters exists for one type of parameter (e.g. β coefficients stored in the array). Other parameters can be stored as regular Float64 columns.")
        else

            # Parse columns and retain only the separated ones

            parsedarrays = expandcol(thearrays, Symbol(names(thearrays)[1]))

            parsedarraysCols = select(parsedarrays, findall(col -> eltype(col) <: Float64, eachcol(parsedarrays)))
        end

        #-------- Merge together --------

        posteriorDF = [thefloats parsedarraysCols]
        params = propertynames(posteriorDF)

    elseif isa(model, Chains)

        posteriorDF = DataFrame(model)

        # Retrieve just the model parameter names to avoid the sampler specific columns

        paramsData = DataFrame(MCMCChains.summarize(model))
        params = paramsData[!, :parameters]

    else
        error("`model` must be an object of type `Chains` or `Array`.")
    end

    #------------ Draw the plots ----------------

    gr() # gr backend for graphics

    # Draw plot for each parameter

    myPlotArray = [histhelper(posteriorDF, p, plot_legend) for p in params]

    return myPlotArray
end

"""

    plot_posterior_density(model, plot_legend, args...; kwargs...)

Draw a density plot of sampled parameters for easy interpretation of regression models fit in `Turing.jl` and `Soss.jl`.

Usage:
```julia-repl
plot_posterior_density(model, plot_legend)
```

Details:

Note that to get the function to work, you may need to call it using a splat, such as `plot(plot_posterior_density(model, plot_legend)...)`.

Arguments:

- `model` : The `Turing.jl` or `Soss.jl` model of class `Chains` or `Array` to draw inferences from.
- `plot_legend` : Boolean of whether to add a legend to the plot or not.
"""
function plot_posterior_density(model, plot_legend::Bool, args...; kwargs...)

    if isa(model, Array)

        myarray = DataFrame(model)

        #-------- Wrangle floats --------

        thefloats = select(myarray, findall(col -> eltype(col) <: Float64, eachcol(myarray)))

        #-------- Wrangle arrays --------

        thearrays = select(myarray, findall(col -> eltype(col) <: Array, eachcol(myarray)))
        ncols_arrays = size(thearrays, 2)

        if ncols_arrays > 1
            error("`plot_posterior_intervals` currently only works for models where an array of parameters exists for one type of parameter (e.g. β coefficients stored in the array). Other parameters can be stored as regular Float64 columns.")
        else

            # Parse columns and retain only the separated ones

            parsedarrays = expandcol(thearrays, Symbol(names(thearrays)[1]))

            parsedarraysCols = select(parsedarrays, findall(col -> eltype(col) <: Float64, eachcol(parsedarrays)))
        end

        #-------- Merge together --------

        posteriorDF = [thefloats parsedarraysCols]
        params = propertynames(posteriorDF)

    elseif isa(model, Chains)

        posteriorDF = DataFrame(model)

        # Retrieve just the model parameter names to avoid the sampler specific columns

        paramsData = DataFrame(MCMCChains.summarize(model))
        params = paramsData[!, :parameters]

    else
        error("`model` must be an object of type `Chains` or `Array`.")
    end

    #------------ Draw the plots ----------------

    gr() # gr backend for graphics

    # Draw plot for each parameter

    myPlotArray = [denshelper(posteriorDF, p, plot_legend) for p in params]

    return myPlotArray
end


text:
using Plots
using Plots.PlotMeasures
using Colors
using StatsPlots

include("collect_data.jl")

function plot_D()

    mcolor = ARGB(186/255,12/255,47/255, 0.8)
    pcolor = ARGB(39/255, 56/255, 150/255, 0.75)
    ocolor = :darkolivegreen

    mkl = get_mkl()
    oct = get_octavian()
    psi4  = get_psi4()

    nao = [144, 116, 236, 264, 592, 574, 724, 172, 232, 350, 528, 472, 592, 632, 724, 204, 322, 336, 338, 528, 632, 588]
    perm = sortperm(nao)
    nao = nao[perm]

    p1 = groupedbar([oct mkl psi4], c = [ocolor mcolor pcolor], legend=:topright,
                label=["Fermi - Octavian" "Fermi - MKL" "Psi4"], dpi=1200, orientation = :v,
                bar_width=0.8,xticks=([2, 6, 12, 18, 22],["S2","S6", "S12","S18","S22"]), framestyle=:box)

    ylabel!("Time (seconds)")
    xlabel!("Dimer")
    p1
    p2 = scatter(nao, mkl[perm] ./ psi4[perm], label="", yticks=[0.5, 1.0, 1.5], ylims=[0.3,1.8],
                color=:black, ylabel="Relative time\n(Fermi/Psi4)", framestyle=:box, marker=:diamond, markercolor=mcolor,
                xlabel="Number of basis functions", legend=:bottomleft)

    scatter!(p2, nao, oct[perm] ./ psi4[perm], label="", markercolor = ocolor, marker=:diamond)
    hline!(p2, [1], color=:gray, linestyle=:dot, label="")

    plot(p1, p2, layout = grid(2,1, heights=[0.7, 0.3]), leftmargin=0.8cm)
end

function plot_E()

    mcolor = ARGB(186/255,12/255,47/255, 0.8)
    pcolor = ARGB(39/255, 56/255, 150/255, 0.75)
    ocolor = :palegreen4

    mkl   = 1000 .* raw_average(:fermi, "MKL")
    oct   = 1000 .* raw_average(:fermi, "Octavian")
    psi4  = 1000 .* raw_average(:psi4)

    nao = [144, 116, 236, 264, 592, 574, 724, 172, 232, 350, 528, 472, 592, 632, 724, 204, 322, 336, 338, 528, 632, 588]
    perm = sortperm(nao)
    nao = nao[perm]

    p1 = groupedbar([oct mkl psi4], c = [ocolor mcolor pcolor], legend=:topright,
                label=["Fermi - Octavian" "Fermi - MKL" "Psi4"], dpi=1200, orientation = :v,
                yaxis=:log,
                ylim = (10^-0.5,10^6),
                bar_width=0.8,xticks=([2, 6, 12, 18, 22],["S2","S6", "S12","S18","S22"]), framestyle=:box)

    ylabel!("Time (miliseconds)")
    xlabel!("Dimer")
    p1
    p2 = scatter(nao, mkl[perm] ./ psi4[perm], label="MKL", yticks=[0.5, 1.0, 1.5], ylims=[0.3,1.8],
                color=:black, ylabel="Relative time\n(Fermi/Psi4)", framestyle=:box, marker=:diamond, markercolor=mcolor,
                xlabel="Number of basis functions", legend=:bottomleft, legendfontsize=6, legendfonthalign=:left)

    scatter!(p2, nao, oct[perm] ./ psi4[perm], label="Octavian", markercolor = ocolor, marker=:diamond)
    hline!(p2, [1], color=:gray, linestyle=:dot, label="")

    plot(p1, p2, layout = grid(2,1, heights=[0.7, 0.3]), leftmargin=0.8cm, dpi=1200)

    savefig("dfmp2.png")
end

max_stars_repo_name: FermiQC/Fermi

max_issues_repo_name: FermiQC/Fermi

max_forks_repo_name: FermiQC/Fermi


text:
using CartesianGeneticProgramming
using Cambrian
import Cambrian.mutate
import RDatasets

"""
A simple example demonstrating symbolic regression on the iris dataset. For real
application, should be improved with a validation set, data shuffling, and
lexicase selection.
"""

function data_setup()
    iris = RDatasets.dataset("datasets", "iris")
    X = convert(Matrix, iris[:, 1:4])'
    X = X ./ maximum(X; dims=2)
    r = iris[:, 5].refs
    Y = zeros(maximum(r), size(X, 2))
    for i in 1:length(r)
        Y[r[i], i] = 1.0
    end
    X, Y
end

X, Y = data_setup()

function evaluate(ind::CGPInd, X::AbstractArray, Y::AbstractArray)
    accuracy = 0.0
    for i in 1:size(X, 2)
        out = process(ind, X[:, i])
        if argmax(out) == argmax(Y[:, i])
            accuracy += 1
        end
    end
    [accuracy / size(X, 1)]
end

cfg = get_config("cfg/iris.yaml")
fit(i::CGPInd) = evaluate(i, X, Y)
mutate(i::CGPInd) = goldman_mutate(cfg, i)
e = CGPEvolution(cfg, fit)
run!(e)


text:
## Copyright (c) 2013 Miguel Bazdresch
##
## This file is distributed under the 2-clause BSD License.

function gaston_demo()
    closeall()

    t = -2:0.01:2
    st = sin(10pi*t)
    ct = cos(10pi*t)
    et = exp(abs(t/10))

    # simplest figure
    figure(1)
    addcoords(st)
    llplot()

    # with x coordinates
    figure(2)
    addcoords(t,st)
    llplot()

    # plot configuration
    figure(3)
    c = CurveConf()
    c.legend = "Sinusoidal"
    c. plotstyle = "points"
    c.color = "blue"
    c.marker = "fdmd"
    c.pointsize = 2
    c.linewidth = 1.5
    addcoords(t,st,c)
    llplot()

    # figure configuration
    figure(4)
    a = AxesConf()
    a.title = "Example of figure configuration"
    a.xlabel = "Time (s)"
    a.ylabel = "Amplitude"
    a.box = "bottom left"
    a.axis = "semilogx"
    c.plotstyle = "linespoints"
    addcoords(t,st,c)
    addconf(a)
    llplot()

    # multiple plots
    figure(5)
    c = CurveConf()
    c.legend = "Sin"
    c.color = "black"
    addcoords(t,st,c)
    c.legend = "Cos"
    c.color = "magenta"
    c.plotstyle = "impulses"
    c.linewidth = 0.4
    addcoords(t,ct,c)
    c.legend = "Exp"
    c.color = "red"
    c.plotstyle = "linespoints"
    addcoords(t,et,c)
    a = AxesConf()
    a.xlabel = "Time (s)"
    a.ylabel = "Amplitude"
    a.title = "Multiple plots in the same figure"
    a.box = "outside top horizontal box"
    addconf(a)
    llplot()

    # error bars with ydelta
    y = exp(-(1:.1:4.9))
    figure(6)
    c = CurveConf()
    c.legend = "Random"
    c.plotstyle = "errorbars"
    addcoords(1:40,y,c)
    adderror(0.1*rand(40))
    a = AxesConf()
    a.title = "Error bars (ydelta)"
    addconf(a)
    llplot()

    # error bars with ylow, yhigh
    figure(7)
    c = CurveConf()
    c.legend = "Random"
    c.plotstyle = "errorbars"
    ylow = y - 0.05*rand(40);
    yhigh = y + 0.05*rand(40);
    addcoords(1:40,y,c)
    ylow = y - 0.05*rand(40);
    yhigh = y + 0.05*rand(40);
    adderror(ylow,yhigh)
    a = AxesConf()
    a.title = "Error bars (ylow, yhigh)"
    addconf(a)
    llplot()

    # error lines
    figure(8)
    c = CurveConf()
    c.legend = "Random"
    c.plotstyle = "errorlines"
    addcoords(1:40,y,c)
    adderror(0.1*rand(40))
    a = AxesConf()
    a.title = "Error lines (ydelta)"
    addconf(a)
    llplot()

    # plotting columns of matrices
    figure(9)
    Y = hcat(st, ct, et)
    X = hcat(t, t, t)
    addcoords(X,Y)
    a = AxesConf()
    a.title = "Plotting matrix columns"
    addconf(a)
    llplot()

    # dots plotstyle
    figure(10)
    x = 1 + 0.1*randn(200)
    y = 1 + 0.1*randn(200)
    c = CurveConf()
    c.plotstyle = "points"
    c.marker = "*"
    c.pointsize = 3
    c.color = "red"
    addcoords([1],[1],c)
    c = CurveConf()
    c.plotstyle = "dots"
    c.color = "blue"
    addcoords(x,y,c)
    a = AxesConf()
    a.title = "dots demo"
    addconf(a)
    llplot()

    # simple 3-D plot with default config
    figure(11)
    x=[0,1,2,3]
    y=[0,1,2]
    Z=[10 10 10; 10 5 10;10 1 10; 10 0 10]
    addcoords(x,y,Z)
    a = AxesConf()
    a.title = "3D: Valley of the Gnu from gnuplot manual"
    addconf(a)
    llplot()

    # same plot with colored surfaces
    figure(12)
    c = CurveConf()
    c.plotstyle = "pm3d"
    addcoords(x,y,Z,c)
    a = AxesConf()
    a.title = "3D: Valley of the Gnu with pm3d"
    addconf(a)
    llplot()

    # sombrero
    figure(13)
    c = CurveConf()
    c.plotstyle = "pm3d"
    x = -15:0.33:15
    y = -15:0.33:15
    Z = meshgrid(x,y,(x,y)->sin(sqrt(x.*x+y.*y))/sqrt(x.*x+y.*y))
    addcoords(x,y,Z,c)
    a = AxesConf()
    a.title = "3D: Sombrero"
    addconf(a)
    llplot()

    # simple image
    figure(14)
    c = CurveConf()
    c.plotstyle = "image"
    Z = [5 4 3 1 0; 2 2 0 0 1; 0 0 0 1 0; 0 1 2 4 3]
    addcoords([],[],Z,c)
    a = AxesConf()
    a.title = "Image"
    addconf(a)
    llplot()

    # rgb image
    figure(15)
    c = CurveConf()
    c.plotstyle = "rgbimage"
    R = [ x+y for x=0:5:120, y=0:5:120]
    G = [ x+y for x=0:5:120, y=120:-5:0]
    B = [ x+y for x=120:-5:0, y=0:5:120]
    Z = zeros(25,25,3)
    Z[:,:,1] = R
    Z[:,:,2] = G
    Z[:,:,3] = B
    addcoords([],[],Z,c)
    a = AxesConf()
    a.title = "RGB Image"
    addconf(a)
    llplot()

    # image with x,y coordinates
    figure(16)
    c = CurveConf()
    c.plotstyle = "rgbimage"
    addcoords(linspace(-2.0,5.0,25),linspace(4.0,6.0,25),Z,c)
    a = AxesConf()
    a.title = "RGB Image with explicit (x,y) coordinates"
    addconf(a)
    llplot()

    # histograms
    figure(17)
    c = CurveConf()
    c.plotstyle = "boxes"
    c.color = "blue"
    y = [1 2 3 4 5 6 7 8 9 10]
    (x,y) = histdata(y,10)
    addcoords(x,y,c)
    a = AxesConf()
    a.title = "Simple histogram test"
    addconf(a)
    llplot()

    figure(18)
    c = CurveConf()
    c.plotstyle = "boxes"
    c.color = "blue"
    c.legend = "1000 samples"
    (x,y) = histdata(randn(1000),25)
    delta = x[2]-x[1]
    y = y/(delta*sum(y))  # normalization
    addcoords(x,y,c)
    x = -5:0.05:5
    y = 1/sqrt(2pi)*exp((-x.^2)/2)
    c = CurveConf()
    c.plotstyle = "lines"
    c.color = "black"
    c.legend = "Theoretical"
    addcoords(x,y,c)
    a = AxesConf()
    a.title = "Histogram"
    addconf(a)
    llplot()

    # image and curve on the same figure
    figure(19)
    c = CurveConf()
    c.plotstyle = "image"
    Z = [5 4 3 1 0; 2 2 0 0 1; 0 0 0 1 0; 0 1 2 4 3]
    addcoords([],[],Z,c)
    c = CurveConf()
    c.plotstyle = "lines"
    c.color = "blue"
    t = linspace(0,2pi,10)
    addcoords(3*sin(t)+3,c)
    a = AxesConf()
    a.title = "Image and curve on the same figure"
    addconf(a)
    llplot()

    # Two surfaces in one figure, plus low-level use of gnuplot_send
    # for final tweaking
    figure(20)
    gnuplot_send("set view 69,20")
    c = CurveConf()
    c.plotstyle = "pm3d"
    x = -15:0.33:15
    y = -15:0.33:15
    Z = meshgrid(x,y,(x,y)->sin(sqrt(x.*x+y.*y))/sqrt(x.*x+y.*y))
    addcoords(x,y,Z,c)
    Z = meshgrid(x,y,(x,y)->cos(x/2)*sin(y/2)+3)
    addcoords(x,y,Z,c)
    a = AxesConf()
    a.title = "3D: Two surfaces on the same figure"
    addconf(a)
    llplot()

end


text:
"""
    BuckleyLeverette{T,Dim}

A struct representing the Buckley-Leverette equation in `Dim` dimensions using
`T` as scalar type.
"""
struct BuckleyLeverette{T} <: ScalarBalanceLaw{T,1} end

function BuckleyLeverette(T=Float64)
  BuckleyLeverette{T}()
end

function Base.show(io::IO, model::BuckleyLeverette{T}) where {T}
  print(io, "Buckley-Leverette equation {T=", T, "}",
            " with flux f(u) = u^2 / (u^2 + (1-u)^2)")
end

"""
    flux(u, model::BuckleyLeverette)

Compute the flux of `u` for `model`.
"""
@inline flux(u, model::BuckleyLeverette) = u^2 / (u^2 + (1-u)^2)

"""
    speed(u::Real, model::BuckleyLeverette)

Compute the speed f'(`u`) for `model`.
"""
@inline function speed(u, model::BuckleyLeverette)
  2*u*(1-u) / (u^2 + (1-u)^2)^2
end

"""
    max_abs_speed(uₗ, uᵣ, model::BuckleyLeverette)

Compute the maximal absolute value of the speed in the solution of the Riemann
problem with states `uₗ`, `uᵣ` for `model`.
"""
@inline function max_abs_speed(uₗ, uᵣ, model::BuckleyLeverette)
    # the flux is not convex; the maximal speed is 0.5 at u=0.5
    critical_u = one(uₗ) / 2
    if min(uₗ,uᵣ) < critical_u && max(uₗ,uᵣ) > critical_u
        λ = max_abs_speed(critical_u, model)
    else
        λ = max(max_abs_speed(uₗ,model), max_abs_speed(uᵣ,model))
    end
end

"""
    min_max_speed(uₗ, uᵣ, model::BuckleyLeverette)

Compute the minimal and maximal speed in the solution of the Riemann problem with
states `uₗ`, `uᵣ` for `model`.
"""
@inline function min_max_speed(uₗ, uᵣ, model::BuckleyLeverette)
    # the flux is not convex; the maximal speed is 0.5 at u=0.5
    critical_u = one(uₗ) / 2
    λ₋, λ₊ = minmax(speed(uₗ, model), speed(uᵣ, model))
    λ₊ = max(speed(critical_u, model), λ₊)
    λ₋, λ₊
end
#= TODO: Something else?
@inline function min_max_speed(uₗ, uᵣ, model::BuckleyLeverette)
    # the flux is not convex; the maximal speed is 0.5 at u=0.5
    critical_u = one(uₗ) / 2
    umin, umax = minmax(uₗ, uᵣ)
    if umin < critical_u && umax > critical_u
        λ₋ =
        λ₊ = speed(critical_u, model)
    else
        λ₋, λ₊ = minmax(speed(uₗ, model), speed(uᵣ, model))
    end
    λ₋, λ₊
end
=#

################################################################################

function (::GodunovFlux)(uₗ, uᵣ, model::BuckleyLeverette)
  flux(uₗ, model)
end

################################################################################

"""
    BuckleyLeveretteRiemannSolution{T,T1}

The solution of a Riemann problem `prob` for the Buckley-Leverette equation.
"""
struct BuckleyLeveretteRiemannSolution{T,T1} <: ScalarRiemannSolution
  prob::RiemannProblem{BuckleyLeverette{T},T,T1}
  σ⁻::T
  σ⁺::T
end

"""
    minmax_speeds(sol::BuckleyLeveretteRiemannSolution)

Return the minimal and maximal speeds `σ⁻, σ⁺` that appear in the solution `sol`.
"""
function minmax_speeds(sol::BuckleyLeveretteRiemannSolution)
  sol.σ⁻, sol.σ⁺
end

"""
    (sol::BuckleyLeveretteRiemannSolution)(ξ::Real)

Evaluate the solution `sol` at the value `ξ` of the self-similarity variable
`ξ = (x - x₀) / (t - t₀)`.
"""
function (sol::BuckleyLeveretteRiemannSolution)(ξ::Real)
  @unpack σ⁻, σ⁺ = sol
  @unpack uₗ, uᵣ = sol.prob

  if ξ < σ⁻
    uₗ
  elseif ξ < σ⁺
    # rarefaction wave
    if uₗ < uᵣ
      ( 1 - sqrt(-1 + 4/(1+sqrt(1+4ξ))) ) / 2
    else
      ( 1 + sqrt(-1 + 4/(1+sqrt(1+4ξ))) ) / 2
    end
  else
    uᵣ
  end
end

"""
    (sol::BuckleyLeveretteRiemannSolution)(t::Real, x::Real)

Evaluate the solution `sol` at the time and space coordinates `t` and `x`.
"""
function (sol::BuckleyLeveretteRiemannSolution)(t::Real, x::Real)
  @unpack x₀, t₀ = sol.prob

  sol((x-x₀)/(t-t₀))
end

"""
    shockspeed(uₗ, uᵣ, model::BuckleyLeverette)

Compute the speed f a shock with left and right state `uₗ`, `uᵣ` for `model`.
"""
function shockspeed(uₗ, uᵣ, model::BuckleyLeverette)
  (uₗ + uᵣ - 2*uₗ*uᵣ) / ( (1 + 2*uₗ*(uₗ-1)) * (1 + 2*uᵣ*(uᵣ-1)) )
end

"""
    solve(prob::RiemannProblem{BuckleyLeverette{T},T,T1}) where {T,T1}

Compute the solution of the Riemann prolem `prob`.
"""
function solve(prob::RiemannProblem{BuckleyLeverette{T},T,T1}) where {T,T1}
  @unpack uₗ, uᵣ, model = prob
  u_crit = one(T) / 2

  if uₗ < 0 || uₗ > 1 || uᵣ < 0 || uᵣ > 1
    error("The values uₗ, uᵣ are not in the valid interval [0,1].")
  end

  # f(u) = u^2 / (u^2 + (1-u)^2) is convex on [0, 0.5] and concave on [0.5, 1]
  if uₗ ≈ uᵣ
    σ⁻ = σ⁺ = zero(T)
  elseif uₗ < uᵣ
    # find the lower convex envelope
    if uₗ <= u_crit && uᵣ <= u_crit
      # f is convex
      # the solutions contains a single rarefaction wave
      σ⁻ = speed(uₗ, model)
      σ⁺ = speed(uᵣ, model)
    elseif uₗ >= u_crit && uᵣ >= u_crit
      # f is concave
      # the solutions contains a single shock wave
      σ⁻ = σ⁺ = shockspeed(uₗ, uᵣ, model)
    else
      # f has both convex and concave parts
      σ⁻ = speed(uₗ, model)
      umid = 1 / (2uᵣ + sqrt(2 + 4uᵣ*(uᵣ-1)))
      σ⁺ = speed(umid, model)
    end
  else # uₗ > uᵣ
    # find the upper concave envelope
    if uₗ <= u_crit && uᵣ <= u_crit
      # f is convex
      # the solutions contains a single shock
      σ⁻ = σ⁺ = shockspeed(uₗ, uᵣ, model)
    elseif uₗ >= u_crit && uᵣ >= u_crit
      # f is concave
      # the solutions contains a single rarefaction wave
      σ⁻ = speed(uₗ, model)
      σ⁺ = speed(uᵣ, model)
    else
      # f has both convex and concave parts
      σ⁻ = speed(uₗ, model)
      umid = 1 / (2uᵣ + sqrt(2 + 4uᵣ*(uᵣ-1)))
      σ⁺ = speed(umid, model)
    end
  end

  BuckleyLeveretteRiemannSolution(prob, σ⁻, σ⁺)
end


text:
@doc doc"""The environment variable for APP_SUBMIT_TIME. Set in AppMaster environment only""" ->
const APP_SUBMIT_TIME_ENV = "APP_SUBMIT_TIME_ENV"

@doc doc"""The cache file into which container token is written""" ->
const CONTAINER_TOKEN_FILE_ENV_NAME = "HADOOP_TOKEN_FILE_LOCATION"

@doc doc"""# The environmental variable for APPLICATION_WEB_PROXY_BASE. Set in ApplicationMaster's environment only. This states that for all non-relative web URLs in the app masters web UI what base should they have.""" ->
const APPLICATION_WEB_PROXY_BASE_ENV = "APPLICATION_WEB_PROXY_BASE"

@doc doc"""# The temporary environmental variable for container log directory. This should be replaced by real container log directory on container launch.""" ->
const LOG_DIR_EXPANSION_VAR = "<LOG_DIR>"

@doc doc"""# The environment variable for MAX_APP_ATTEMPTS. Set in AppMaster environment only""" ->
const MAX_APP_ATTEMPTS_ENV = "MAX_APP_ATTEMPTS"

# Environment for Applications
const USER = "USER"
const LOGNAME = "LOGNAME"
const HOME = "HOME"
const PWD = "PWD"
const PATH = "PATH"
const SHELL = "SHELL"
const JAVA_HOME = "JAVA_HOME"
const CLASSPATH = "CLASSPATH"
const APP_CLASSPATH = "APP_CLASSPATH"
const LD_LIBRARY_PATH = "LD_LIBRARY_PATH"
const HADOOP_CONF_DIR = "HADOOP_CONF_DIR"
const HADOOP_COMMON_HOME = "HADOOP_COMMON_HOME"
const HADOOP_HDFS_HOME = "HADOOP_HDFS_HOME"
const MALLOC_ARENA_MAX = "MALLOC_ARENA_MAX"
const HADOOP_YARN_HOME = "HADOOP_YARN_HOME"
const CONTAINER_ID = "CONTAINER_ID"
const NM_HOST = "NM_HOST"
const NM_HTTP_PORT = "NM_HTTP_PORT"
const NM_PORT = "NM_PORT"
const LOCAL_DIRS = "LOCAL_DIRS"
const LOG_DIRS = "LOG_DIRS"


text:
## activate project environment
using Pkg
Pkg.activate(@__DIR__)
Pkg.instantiate()

## compute relative performance of inlining
using DelimitedFiles
using Measurements

function read_data(filename)
  header = readline(filename)
  data = readdlm(filename, comments=true)
  polydegs = data[:, 1]
  results = Vector{Measurement{eltype(data)}}[]
  col = 2
  while col + 1 <= size(data, 2)
    push!(results, data[:, col] .± data[:, col+1])
    col += 2
  end

  header, polydegs, results
end

function compute_relative_performance(file_new, file_baseline)
  header, polydegs, results_baseline = read_data(file_baseline)
  header_new, polydegs_new, results_new = read_data(file_new)
  @assert header == header_new
  @assert polydegs == polydegs_new
  @assert size(results_baseline) == size(results_new)

  data = zeros(length(polydegs_new),
               1 + 2 * length(results_baseline))
  data[:, 1] = polydegs
  for i in eachindex(results_baseline, results_new)
    ratio = results_new[i] ./ results_baseline[i]
    data[:, 2 * i] = Measurements.value.(ratio)
    data[:, 2 * i + 1] = Measurements.uncertainty.(ratio)
  end

  open(file_baseline[1:end-4] * "_relative.dat", "w") do io
    println(io, header)
    writedlm(io, data)
  end
end

setup = "2D_flux_shima_etal"
compute_relative_performance(
  joinpath(dirname(@__DIR__), "Cartesian_vs_curved", "pids_$(setup).dat"),
  joinpath(@__DIR__, "pids_$(setup)_notinlined.dat"))

setup = "2D_flux_ranocha"
compute_relative_performance(
  joinpath(dirname(@__DIR__), "Cartesian_vs_curved", "pids_$(setup).dat"),
  joinpath(@__DIR__, "pids_$(setup)_notinlined.dat"))

setup = "3D_flux_shima_etal"
compute_relative_performance(
  joinpath(dirname(@__DIR__), "Cartesian_vs_curved", "pids_$(setup).dat"),
  joinpath(@__DIR__, "pids_$(setup)_notinlined.dat"))

setup = "3D_flux_ranocha"
compute_relative_performance(
  joinpath(dirname(@__DIR__), "Cartesian_vs_curved", "pids_$(setup).dat"),
  joinpath(@__DIR__, "pids_$(setup)_notinlined.dat"))


text:
##########################
# FunctionCell Interface #
##########################
"""
	ppca([;kwargs])

Constructs an untrained cell that when piped data inside, calculates the
probabilistic PCA projection matrix of the input data. We assume that `d`
is the dimesionality of the input dataset, `n` the number of samples

# Keyword arguments (same as in `MultivariateStats`)
  * `method` is the core algorithm used: `:ml`, `:em` and `:bayes` (default `:ml`)
  * `maxoutdim` is the number of output dimensions (default `<number of variables>-1`)
  * `mean` is the mean vector; can be `nothing`, `0` or precomputed mean (default `nothing`)
  * `tol` convergence tolerance (default `1.0e-6`)
  * `tot` maximum number of iterations (default `1000`)

Read the `MultivariateStats.jl` documentation for more information.
"""
ppca(;kwargs...) = FunctionCell(ppca, (), ModelProperties(), kwtitle("PPCA", kwargs);kwargs...)

############################
# DataCell/Array Interface #
############################
"""
	ppca(x, [;kwargs])

Trains a function cell that when piped data into will compute the principal components of
the observations based on the projection matrix extracted from `x`.
"""
# Training
ppca(x::T where T<:CellData; kwargs...) = ppca(getx!(x);kwargs...)
ppca(x::T where T<:AbstractVector; kwargs...) = ppca(mat(x, LearnBase.ObsDim.Constant{2}());kwargs...)
ppca(x::T where T<:AbstractMatrix; kwargs...)  = begin

	ppcadata = fit(MultivariateStats.PPCA, getobs(x); kwargs...)

	# Build model properties
	modelprops = ModelProperties(size(ppcadata.W,1), size(ppcadata.W,2))

	# Returned trained cell
	FunctionCell(ppca, Model(ppcadata, modelprops), kwtitle("PPCA",kwargs))
end

# Execution
ppca(x::T where T<:CellData, model::Model{<:MultivariateStats.PPCA}) =
	datacell(ppca(getx!(x), model), gety(x))
ppca(x::T where T<:AbstractVector, model::Model{<:MultivariateStats.PPCA}) =
	ppca(mat(x, LearnBase.ObsDim.Constant{2}()), model)
ppca(x::T where T<:AbstractMatrix, model::Model{<:MultivariateStats.PPCA}) =
	MultivariateStats.transform(model.data, getobs(x))

"""
	ppcar(M::Cell), where `M=ppca(x)`

Trains a cell based on a previously trained PPCA function cell. When piped data into,
the returned output will be the reconstruction of the observations in the original
space. Obviously, the input data (e.g. data to be reconstructed) is expected to be
principal components extracted with `M`.

Read the `MultivariateStats.jl` documentation for more information.
"""
# Training for the reconstruction tranform (fixed cell)
ppcar(x::T where T<:CellFunT{<:Model{<:MultivariateStats.PPCA}}) = begin

	# Build model properties

	modelprops = ModelProperties(getx(x).properties.odim, getx(x).properties.idim)

	FunctionCell(ppcar, Model(getx(x).data, modelprops), "PPCA: reconstruct")
end

# Execution
ppcar(x::T where T<:CellData, model::Model{<:MultivariateStats.PPCA}) =
	datacell(ppcar(getx!(x), model), gety(x))
ppcar(x::T where T<:AbstractVector, model::Model{<:MultivariateStats.PPCA}) =
	ppcar(mat(x, LearnBase.ObsDim.Constant{2}()), model)
ppcar(x::T where T<:AbstractMatrix, model::Model{<:MultivariateStats.PPCA})::Matrix{Float64} =
	MultivariateStats.reconstruct(model.data, getobs(x))


text:
# ---
# title: JuliaRL\_BasicDQN\_SingleRoomUndirected
# cover: assets/JuliaRL_BasicDQN_SingleRoomUndirected.png
# description: A simple example to demonstrate how to use environments in GridWorlds.jl
# date: 2021-07-27
# author: "[Siddharth Bhatia](https://github.com/Sid-Bhatia-0)"
# ---

#+ tangle=true
using ReinforcementLearning
using GridWorlds
using StableRNGs
using Flux
using Flux.Losses

function RL.Experiment(
    ::Val{:JuliaRL},
    ::Val{:BasicDQN},
    ::Val{:SingleRoomUndirected},
    ::Nothing;
    seed=123,
)
    rng = StableRNG(seed)

    env = GridWorlds.SingleRoomUndirectedModule.SingleRoomUndirected(rng=rng)
    env = GridWorlds.RLBaseEnv(env)
    env = RLEnvs.StateTransformedEnv(env; state_mapping=x -> vec(Float32.(x)))
    env = RewardTransformedEnv(env; reward_mapping = x -> x - convert(typeof(x), 0.01))
    env = MaxTimeoutEnv(env, 240)

    ns, na = length(state(env)), length(action_space(env))
    agent = Agent(
        policy=QBasedPolicy(
            learner=BasicDQNLearner(
                approximator=NeuralNetworkApproximator(
                    model=Chain(
                        Dense(ns, 128, relu; init=glorot_uniform(rng)),
                        Dense(128, 128, relu; init=glorot_uniform(rng)),
                        Dense(128, na; init=glorot_uniform(rng)),
                    ) |> cpu,
                    optimizer=ADAM(),
                ),
                batch_size=32,
                min_replay_history=100,
                loss_func=huber_loss,
                rng=rng,
            ),
            explorer=EpsilonGreedyExplorer(
                kind=:exp,
                ϵ_stable=0.01,
                decay_steps=500,
                rng=rng,
            ),
        ),
        trajectory=CircularArraySARTTrajectory(
            capacity=1000,
            state=Vector{Float32} => (ns,),
        ),
    )

    stop_condition = StopAfterStep(10_000, is_show_progress=!haskey(ENV, "CI"))
    hook = TotalRewardPerEpisode()
    Experiment(agent, env, stop_condition, hook, "")
end

#+ tangle=false
using Plots
pyplot() #hide
ex = E`JuliaRL_BasicDQN_SingleRoomUndirected`
run(ex)
plot(ex.hook.rewards)
savefig("assets/JuliaRL_BasicDQN_SingleRoomUndirected.png") #hide

# ![](assets/JuliaRL_BasicDQN_SingleRoomUndirected.png)


text:
abstract type AbstractLayer{N} end

Base.ndims(::AbstractLayer{N}) where {N} = N
Base.size(layer::AbstractLayer, d::Int) = size(layer)[d]

"""
    flatten(layer, x)

Returns a vectorized version of `x`.
"""
function flatten(layer::AbstractLayer, x::AbstractArray)
    @assert size(layer) == size(x)[1:ndims(layer)]
    if ndims(layer) == ndims(x)
        return vec(x)
    else
        return reshape(x, length(layer), prod(size(x, d) for d in (ndims(layer) + 1):ndims(x)))
    end
end

"""
    energy(layer, x)

Layer energy, reduced over layer dimensions.
"""
function energy(layer::AbstractLayer, x::AbstractArray)
    @assert size(layer) == size(x)[1:ndims(layer)]
    Es = energies(layer, x)
    if ndims(layer) == ndims(x)
        return sum(Es)
    else
        E = sum(Es; dims = 1:ndims(layer))
        return reshape(E, size(x)[(ndims(layer) + 1):end])
    end
end

"""
    free_energy(layer, inputs = 0)

Cumulant generating function of layer, reduced over layer dimensions.
"""
function free_energy(layer::AbstractLayer, inputs::Union{Real,AbstractArray} = 0)
    F = free_energies(layer, inputs)
    if inputs isa Real || ndims(layer) == ndims(inputs)
        return sum(F)
    else
        f = sum(F; dims = 1:ndims(layer))
        return reshape(f, size(inputs)[(ndims(layer) + 1):end])
    end
end

"""
    suffstats(layer, data; [wts])

Computes the sufficient statistics of `layer` extracted from `data`.
"""
function suffstats end

"""
    batchdims(layer, x)

Indices of batch dimensions in `x`, with respect to `layer`.
"""
function batchdims(layer::AbstractLayer, x::AbstractArray)
    @assert size(layer) == size(x)[1:ndims(layer)]
    return (ndims(layer) + 1):ndims(x)
end

"""
    batch_size(layer, x)

Batch sizes of `x`, with respect to `layer`.
"""
function batch_size(layer::AbstractLayer, x::AbstractArray)
    @assert size(layer) == size(x)[1:ndims(layer)]
    return size(x)[batchdims(layer, x)]
end

"""
    batchmean(layer, x; wts = nothing)

Mean of `x` over batch dimensions, weigthed by `wts`.
"""
function batchmean(layer::AbstractLayer, x::AbstractArray; wts = nothing)
    @assert size(layer) == size(x)[1:ndims(layer)]
    if ndims(layer) == ndims(x)
        wts::Nothing
        return x
    else
        μ = wmean(x; wts, dims = batchdims(layer, x))
        return reshape(μ, size(layer))
    end
end

"""
    batchvar(layer, x; wts = nothing, [mean])

Variance of `x` over batch dimensions, weigthed by `wts`.
"""
function batchvar(
    layer::AbstractLayer, x::AbstractArray; wts = nothing, mean = batchmean(layer, x; wts)
)
    @assert size(layer) == size(x)[1:ndims(layer)] == size(mean)
    return batchmean(layer, (x .- mean).^2; wts)
end

"""
    batchstd(layer, x; wts = nothing, [mean])

Standard deviation of `x` over batch dimensions, weigthed by `wts`.
"""
function batchstd(
    layer::AbstractLayer, x::AbstractArray; wts = nothing, mean = batchmean(layer, x; wts)
)
    return sqrt.(batchvar(layer, x; wts, mean))
end

"""
    batchcov(layer, x; wts = nothing, [mean])

Covariance of `x` over batch dimensions, weigthed by `wts`.
"""
function batchcov(
    layer::AbstractLayer, x::AbstractArray; wts = nothing, mean = batchmean(layer, x; wts)
)
    @assert size(layer) == size(x)[1:ndims(layer)] == size(mean)
    ξ = flatten(layer, x .- mean)
    if isnothing(wts)
        C = ξ * ξ' / size(ξ, 2)
    else
        @assert size(wts) == batch_size(layer, x)
        w = Diagonal(vec(wts))
        C = ξ * w * ξ' / sum(w)
    end
    return reshape(C, size(layer)..., size(layer)...)
end

"""
    total_mean_from_inputs(layer, inputs; wts = nothing)

Total mean of unit activations from inputs.
"""
function total_mean_from_inputs(
    layer::AbstractLayer, inputs::Union{Real,AbstractArray}; wts = nothing
)
    h_ave = mean_from_inputs(layer, inputs)
    return batchmean(layer, h_ave; wts)
end

"""
    total_var_from_inputs(layer, inputs; wts = nothing)

Total variance of unit activations from inputs.
"""
function total_var_from_inputs(
    layer::AbstractLayer, inputs::Union{Real,AbstractArray}; wts = nothing
)
    h_ave, h_var = meanvar_from_inputs(layer, inputs)
    ν_int = batchmean(layer, h_var; wts) # intrinsic noise
    ν_ext = batchvar(layer, h_ave; wts) # extrinsic noise
    return ν_int + ν_ext # law of total variance
end

"""
    total_meanvar_from_inputs(layer, inputs; wts = nothing)

Total mean and total variance of unit activations from inputs.
"""
function total_meanvar_from_inputs(
    layer::AbstractLayer, inputs::Union{Real,AbstractArray}; wts = nothing
)
    h_ave, h_var = meanvar_from_inputs(layer, inputs)
    μ = batchmean(layer, h_ave; wts)
    ν_int = batchmean(layer, h_var; wts) # intrinsic noise
    ν_ext = batchvar(layer, h_ave; wts, mean = μ) # extrinsic noise
    ν = ν_int + ν_ext # law of total variance
    return (μ = μ, ν = ν)
end

"""
    ∂energy(layer, data; wts = nothing)
    ∂energy(layer, stats)

Derivative of average energy of `data` with respect to `layer` parameters.
In the second form, `stats` are the pre-computed sufficient statistics
of the layer. See [`suffstats`](@ref).
"""
function ∂energy(layer::AbstractLayer, data::AbstractArray; wts=nothing)
    stats = suffstats(layer, data; wts)
    return ∂energy(layer, stats)
end

"""
    ∂free_energy(layer, inputs = 0; wts = 1)

Unit activation moments, conjugate to layer parameters.
These are obtained by differentiating `free_energies` with respect to the layer parameters.
Averages over configurations (weigthed by `wts`).
"""
function ∂free_energy(layer::AbstractLayer, inputs::Union{Real,AbstractArray} = 0; wts = nothing)
    ∂F = ∂free_energies(layer, inputs)
    return map(∂F) do ∂f
        batchmean(layer, ∂f; wts)
    end
end


text:
import StateSpaceReconstruction:
    DelaunayTriangulation,
    Simplex,
    centroid,
    radius

"""
idxs_potentially_intersecting_simplices(all_pts::Vector{Vector{T}},
        DT::DelaunayTriangulation, idx::Int) where T

Given a delaunay triangulation `DT` of `all_pts[1:(end - 1)]`, find the indices of
the simplices that potentially intersect with the image simplex with index `idx`.
"""
function idxs_potentially_intersecting_simplices(all_pts,
	        DT::DelaunayTriangulation, idx::Int) where T
	# Vector that will hold the indices of the simplices potentially
	# intersecting with the image of simplex #idx
	inds_potential_simplices = Int[]

	n_simplices = length(DT)

	original_pts = all_pts[1:(end - 1)]
	forward_pts = all_pts[2:end]

	simplices = [ Simplex(original_pts[DT[i]]) for i = 1:n_simplices]
	image_simplices = [Simplex(forward_pts[DT[i]]) for i = 1:n_simplices]

	cs = [centroid(simplices[i]) for i = 1:n_simplices]
	cs_im = [centroid(image_simplices[i]) for i = 1:n_simplices]

	rs = [radius(simplices[i]) for i = 1:n_simplices]
	rs_im = [radius(image_simplices[i]) for i = 1:n_simplices]

	@inbounds for i = 1:n_simplices
	    δ = transpose(cs_im[idx] .- cs[i]) * (cs_im[idx] .- cs[i]) .- ((rs_im[idx] + rs[i])^2)
	    if δ[1] < 0
	        push!(inds_potential_simplices, i)
	    end
	end

	return inds_potential_simplices
end


text:
#################################################################################
#               Genera gráfica de las funciones de activación
#################################################################################
# Paquetes
using Plots
using TOML
include("../../OptimizedNeuralNetwork.jl/src/activation_functions.jl")
using .ActivationFunctions

FICHERO_CONFIGURACION = "Experimentos/.config.toml"
config = TOML.parsefile(FICHERO_CONFIGURACION)["visualizacion-funciones-activacion"]
img_path = config["DIRECTORIO_IMAGENES"]

### Funciones que se van a analizar  y auxiliares
p(x) = 2x +1
id(x)=x
# Estructura caracteríticas a mostrar
struct FuncionActivacion
    nombre::String
    funcion_activacion # función propiamente dicho
    cota_inferior :: Real
    cota_superior :: Real
end

# Declaración de funciones
funciones_a_mostrar = [
    FuncionActivacion("Identidad", id, -1,1),
    FuncionActivacion("Threshold de 2x+1", @ThresholdFunction(p,0), -2, 2),
    FuncionActivacion("Cosine CosineSquasher", CosineSquasher, -2pi, 2pi),
    FuncionActivacion("Indicadora de 0", @IndicatorFunction(0), -2, 2),
    FuncionActivacion("Rampa", RampFunction, -2, 2),
    FuncionActivacion("ReLU", ReLU, -2, 2),
    FuncionActivacion("Sigmoidea", Sigmoid, -4,4),
    FuncionActivacion("Tangente hiperbolica", tanh, -4,4),
    FuncionActivacion("Valor absoluto", abs, -2,2),
    FuncionActivacion("Coseno", cos, -2pi,2pi),
    FuncionActivacion("Hardtanh", HardTanh, -2pi,2pi),
    FuncionActivacion("LReLU", @LReLU(0.01), -2pi,2pi),
]

"""
    MuestraYGuarda(f :: FuncionActivacion, img_path:: String)
Muestra por pantalla una función de activación con la estructura FuncionActivacion
y la guarda en el directorio indicado en `img_path`
"""
function MuestraYGuarda(fa :: FuncionActivacion, img_path:: String)
    plot( fa.funcion_activacion, fa.cota_inferior, fa.cota_superior, title="", label="",linewidth=3)
    savefig(img_path*fa.nombre)
end

# Mostramos y guardamos para las funciones definidas
for fa in funciones_a_mostrar
   MuestraYGuarda(fa, img_path)
end


text:
module DiscO

using Random, Distributions, DataFrames
include("fitnesses.jl")

export onemax, zeromax, fractionmax

include("mutations.jl")

export staticBitMutation, flip

include("db.jl")

export db, setup, finish, store, info, record, result, meta, dbEngine

include("ea.jl")

export ea1p1

end # module

max_stars_repo_name: ahrajabi/DiscO

max_issues_repo_name: ahrajabi/DiscO

max_forks_repo_name: ahrajabi/DiscO


text:
macro loop_xyz(i, j, k, grid, expr)
    return esc(
        quote
            @loop for $k in (1:$grid.Nz; (blockIdx().z - 1) * blockDim().z + threadIdx().z)
                @loop for $j in (1:$grid.Ny; (blockIdx().y - 1) * blockDim().y + threadIdx().y)
                    @loop for $i in (1:$grid.Nx; (blockIdx().x - 1) * blockDim().x + threadIdx().x)
                        $expr
                    end
                end
            end
        end)
end

macro loop_xy(i, j, grid, expr)
    return esc(
        quote
            @loop for $j in (1:$grid.Ny; (blockIdx().y - 1) * blockDim().y + threadIdx().y)
                @loop for $i in (1:$grid.Nx; (blockIdx().x - 1) * blockDim().x + threadIdx().x)
                    $expr
                end
            end
        end)
end

macro loop_xz(i, k, grid, expr)
    return esc(
        quote
            @loop for $k in (1:$grid.Nz; (blockIdx().y - 1) * blockDim().y + threadIdx().y)
                @loop for $i in (1:$grid.Nx; (blockIdx().x - 1) * blockDim().x + threadIdx().x)
                    $expr
                end
            end
        end)
end


text:
__precompile__()

module Dummy1AF

using ComputationalResources, Dummy1

Dummy1.foo(::ArrayFireLibs) = true

end