text:
include("utils.jl")

@testset "Sphere" begin
    M = Sphere(2)
    @testset "Sphere Basics" begin
        @test representation_size(M) == (3,)
        @test !is_manifold_point(M, [1., 0., 0., 0.])
        @test !is_tangent_vector(M, [1.,0.,0.], [0., 0., 1., 0.])
    end
    types = [Vector{Float64},
             SizedVector{3, Float64},
             MVector{3, Float64},
             Vector{Float32},
             SizedVector{3, Float32},
             MVector{3, Float32},
             Vector{Double64},
             MVector{3, Double64},
             SizedVector{3, Double64}]
    for T in types
        @testset "Type $T" begin
            pts = [convert(T, [1.0, 0.0, 0.0]),
                   convert(T, [0.0, 1.0, 0.0]),
                   convert(T, [0.0, 0.0, 1.0])]
            test_manifold(M,
                          pts,
                          test_reverse_diff = isa(T, Vector),
                          test_project_tangent = true,
                          test_musical_isomorphisms = true,
                          test_vector_transport = true,
                          test_mutating_rand = isa(T, Vector),
                          point_distributions = [Manifolds.uniform_distribution(M, pts[1])],
                          tvector_distributions = [Manifolds.normal_tvector_distribution(M, pts[1], 1.0)])

            @test isapprox(-pts[1], exp(M, pts[1], log(M, pts[1], -pts[1])))
        end
    end

    @testset "Distribution tests" begin
        usd_mvector = Manifolds.uniform_distribution(M, @MVector [1.0, 0.0, 0.0])
        @test isa(rand(usd_mvector), MVector)

        gtsd_mvector = Manifolds.normal_tvector_distribution(M, (@MVector [1.0, 0.0, 0.0]), 1.0)
        @test isa(rand(gtsd_mvector), MVector)
    end

    @testset "log edge case" begin
        n = manifold_dimension(M)
        x = normalize(randn(n + 1))
        v = log(M, x, -x)
        @test norm(v) ≈ π
        @test isapprox(dot(x, v), 0; atol=1e-12)
        vexp = normalize(project_tangent(M, x, [1, zeros(n)...]))
        @test v ≈ π * vexp

        x = [1, zeros(n)...]
        v = log(M, x, -x)
        @test norm(v) ≈ π
        @test isapprox(dot(x, v), 0; atol=1e-12)
        vexp = normalize(project_tangent(M, x, [0, 1, zeros(n - 1)...]))
        @test v ≈ π * vexp
    end

end


text:
using Documenter
using CellFreeModelGenerationKit

makedocs(
    sitename = "CellFreeModelGenerationKit",
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true"
    ),
    modules = [CellFreeModelGenerationKit],
    doctest = true,
    pages = [
        "Home" => "index.md",
        "Quick start" => "Quick.md",
        "Files" => "generated_files.md",
        "VFF" => "vff_format.md",
        "Advanced" => "parser.md"
    ],
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(
    repo = "github.com/varnerlab/CellFreeModelGenerationKit.jl.git",
    devbranch = "main",
    devurl = "dev",
)

ext: jl

max_stars_repo_name: varnerlab/CellFreeModelGenerationKit

max_issues_repo_name: varnerlab/CellFreeModelGenerationKit

max_forks_repo_name: varnerlab/CellFreeModelGenerationKit


text:
@inline function read_type_swap(io, ::Type{T}) where {T<:Tuple}
    return ((read_type_swap(io, T_i) for T_i in T.parameters)...,)
end

@inline function read_type_noswap(io, ::Type{T}) where {T<:Tuple}
    return ((read_type_noswap(io, T_i) for T_i in T.parameters)...,)
end

read_type_noswap(io, ::Type{T}) where {T} = read(io, T)
read_type_swap(io, ::Type{T}) where {T} = bswap(read(io, T))

@inline function read_type(io, ::Type{T}, needswap::Bool) where {T}
    if needswap
        return read_type_swap(io, T)
    else
        return read_type_noswap(io, T)
    end
end

function swap_array(A, needswap::Bool)
    if needswap
        return mappedarray(ntoh, hton, A)
    else
        return A
    end
end

function read_array(
    io,
    ::Type{T},
    sz::Tuple{Vararg{<:Integer,N}},
    needswap::Bool,
    mmap::Bool=true,
    mode::AbstractString="r";
    grow::Bool=true,
    shared::Bool=true
) where {T,N}

    if mmap
        return swap_array(Mmap.mmap(io, Array{T,N}, sz, grow=grow, shared=shared), needswap)
    else
        return swap_array(read!(io, Array{T,N}(undef, sz)), needswap)
    end
end

#=

@inline function read_type_noswap(io, ::Type{T}) where {T<:Tuple}
    return ((read_type_noswap(io, T_i) for T_i in T.parameters)...,)
end
@inline read_type_noswap(io, ::Type{T}) where {T} = read(io, T)

@inline function read_type_swap(io, ::Type{T}) where {T<:Tuple}
    return ((read_type_swap(io, T_i) for T_i in T.parameters)...,)
end
@inline read_type_swap(io, ::Type{T}) where {T} = bswap(read(io, T))
=#


text:
export MaxMetricGoal

"Goal specification where costs are differences in a state-based metric."
struct MaxMetricGoal <: Goal
    terms::Vector{Term} # Goal terms to be satisfied
    metric::Term # Reward metric to be maximized
end

function MaxMetricGoal(problem::Problem)
    goals = flatten_conjs(problem.goal)
    sign, metric = problem.metric
    if sign > 0 metric = Compound(:-, [metric]) end
    return MaxMetricGoal(goals, metric)
end
MaxMetricGoal(goal::Term, metric::Term) =
    MaxMetricGoal(flatten_conjs(goal), metric)

Base.hash(spec::MaxMetricGoal, h::UInt) =
    hash(spec.metric, hash(Set(spec.terms), h))
Base.isequal(s1::MaxMetricGoal, s2::MaxMetricGoal) =
    s1.metric == s2.metric && Set(s1.terms) == Set(s2.terms)

is_goal(spec::MaxMetricGoal, domain::Domain, state::State) =
    satisfy(spec.terms, state, domain)[1]
is_violated(spec::MaxMetricGoal, domain::Domain, state::State) = false
get_cost(spec::MaxMetricGoal, domain::Domain, s1::State, ::Term, s2::State) =
    s1[domain, spec.metric] - s2[domain, spec.metric]
get_reward(spec::MaxMetricGoal, domain::Domain, s1::State, ::Term, s2::State) =
    s2[domain, spec.metric] - s1[domain, spec.metric]
get_goal_terms(spec::MaxMetricGoal) = spec.terms


text:
# These algorithms are from
#   GARSTKA M., CANNON M., GOULART P.:
#   A clique graph based merging strategy for decomposable sdps.
#   IFAC-PapersOnLine 53, 2 (2020), 7355–7361.

import Base.length
length(cg::CliqueGraph) = size(cg.edge_mat, 1)

# Input:
#   cliques : vector of cliques
#   n       : num nodes of graph
# Output:
#   CliqueGraph : rep of relationship bt cliques
"""
    generate_clique_graph(cliques, n::Integer)

Generates CliqueGraph datastructure for an undirected graph from `cliques` and
the number of nodes `n`.

# Reference
* [A clique graph based merging strategy for decomposable SDPs](https://arxiv.org/pdf/1911.05615) by Michael Garstka, Mark Cannon, Paul Goulart
"""
function generate_clique_graph(cliques, n::Integer)
    m = length(cliques)
    membership_mat = spzeros(Bool, n, m)

    # TODO: edge mat as matrix or as list of tuples?
    edge_mat = zeros(Float64, m, m)
    active_cliques = Set(1:m)

    for (idx, clique) in enumerate(cliques)
        for node in clique
            membership_mat[node, idx] = true
        end
    end
    for i in 1:m, j in i+1:m
        ci = @view(membership_mat[:, i])
        cj = @view(membership_mat[:, j])
        l_int = sum(ci .& cj)
        if l_int > 0
            edge_mat[i,j] = edge_mat[j,i] = weight_function(ci, cj)
        end
    end

    return CliqueGraph(membership_mat, edge_mat, active_cliques)
end

# TODO: this should be a user-defined function
"""
    weight_function(ci, cj)

Defines the weight function used to determine if cliques `ci` and `cj` should be
merged. Cliques are merged if this is positive.

# Reference
* [A clique graph based merging strategy for decomposable SDPs](https://arxiv.org/pdf/1911.05615) by Michael Garstka, Mark Cannon, Paul Goulart
"""
function weight_function(ci, cj)
    l_ci = sum(ci)
    l_cj = sum(cj)
    l_union = sum(ci[k] | cj[k] for k in 1:length(ci))
    # l_union = 0.0
    # for k in 1:length(ci)
    #     # l_ci += ci[k]
    #     # l_cj += cj[k]
    #     l_union += ci[k] .| cj[k]
    # end
    return l_ci^3 + l_cj^3 - l_union^3
end

function _merge_cliques!(cg::CliqueGraph, i::Integer, j::Integer)
    # nodes X cliques
    n, m = size(cg.membership_mat)

    ci = @view(cg.membership_mat[:, i])
    cj = @view(cg.membership_mat[:, j])

    n_rem = sum(ci) < sum(cj) ? i : j
    n_keep = sum(ci) < sum(cj) ? j : i
    c_keep = @view(cg.membership_mat[:,n_keep])
    c_rem =  @view(cg.membership_mat[:,n_rem])

    # Update membership matrix
    c_keep .= ci .| cj
    c_rem .= false
    cg.edge_mat[n_keep,n_rem] = cg.edge_mat[n_rem,n_keep] = 0.0

    # TODO: change to only cliques with nonzero intersection
    for k in 1:m
        k == n_keep && continue
        # Remove edges from n_rem
        cg.edge_mat[k,n_rem] = cg.edge_mat[n_rem,k] = 0.0

        # Update edges to n_keep
        wij = weight_function(c_keep, @view(cg.membership_mat[:,k]))
        cg.edge_mat[k,n_keep] = cg.edge_mat[n_keep,k] = wij
    end

    pop!(cg.active_cliques, n_rem)

end

"""
    merge_cliques!(cg::CliqueGraph; verbose=false)

Merges cliques in `cg` in a greedy fashion starting with the pair with the largest
`weight_function`. Stops when all pairs of cliques have a non-positive weight.

# Reference
* [A clique graph based merging strategy for decomposable SDPs](https://arxiv.org/pdf/1911.05615) by Michael Garstka, Mark Cannon, Paul Goulart
"""
function merge_cliques!(cg::CliqueGraph; verbose=false)
    # nodes X cliques
    n, m = size(cg.membership_mat)
    if verbose
        max_clique_size = maximum(nnz, @view(cg.membership_mat[:,i]) for i in 1:m)
        @info "Starting merging; num cliques = $m. Max size is $max_clique_size."
    end

    max_val, max_ind = findmax(cg.edge_mat)
    while max_val > 0
        #TODO: verbose printing
        i = max_ind[1]
        j = max_ind[2]
        _merge_cliques!(cg, i, j)
        max_val, max_ind = findmax(cg.edge_mat)
    end
    dropzeros!(cg.membership_mat)
    if verbose
        max_clique_size = maximum(nnz, @view(cg.membership_mat[:,i]) for i in cg.active_cliques)
        @info "Finished merging; num cliques = $(length(cg.active_cliques)). Max size is $max_clique_size."
    end
end

"""
    get_cliques(cg::CliqueGraph)

Returns a list of the cliques in CliqueGraph `cg`.
"""
function get_cliques(cg::CliqueGraph)
    return [findnz(cg.membership_mat[:,i])[1] for i in sort(collect(cg.active_cliques))]
end


text: quit() # terminates program normally, with its child processes. See also exit(0).

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
@testset "Sets" begin
    ES = ExtensionalSet  # TODO: Come up with better syntax for defining extensional sets

    @testset "Set Matching" begin
        @symbols x y z

        p1 = ES([x])
        p2 = ES([x, ES([y])])
        p3 = ES([x, y])

        A = ES([∅])
        B = ES([∅, ES([∅])])
        C = ES([∅, ES([∅, ES([∅])])])

        # test the trival pattern
        @test first(set_symbol_matches(x, A)[x]) == A
        @test first(set_symbol_matches(x, B)[x]) == B
        @test first(set_symbol_matches(x, C)[x]) == C

        @test first(set_symbol_matches(p1, A)[x]) == ∅
        @test classify_matches(set_symbol_matches(p1, B)) == Deductive.NO_MATCHES
        @test classify_matches(set_symbol_matches(p1, C)) == Deductive.NO_MATCHES

        m1A, m1B = set_symbol_matches(p2, B), set_symbol_matches(p2, C)
        @test first(m1A[x]) == ∅
        @test first(m1A[y]) == ∅
        @test classify_matches(m1B) == Deductive.NO_MATCHES

        m2A, m2B = set_symbol_matches(p3, A), set_symbol_matches(p3, B)
        @test classify_matches(m2A) == Deductive.NO_MATCHES
        @test classify_matches(m2B) == Deductive.PARTIAL_MATCH
        @test length(m2B[x]) == 2
        @test length(m2B[y]) == 2
    end

    @testset "Mixed Set Matching" begin
        # distinguished from test cases above by the mixing of set logic and logical expressions
        #   ends up being used in proof algorithms
        @symbols a b c x y z
        pattern = Set([x, x → y])  # modus ponens without the conclusion :)
        # premises1 = Set([a, a → b])

        # @test set_matches(pattern, premises1; strict=false)
        # TODO: actually write these tests lol
    end
end


text:
using FastGaussQuadrature

get_nodes(::Val{:gauss_legendre}, n, T=Float64) = convert(Array{T}, gausslegendre(n)[1])
get_nodes(::Val{:lobatto_legendre}, n, T=Float64) = convert(Array{T}, gausslobatto(n)[1])

function get_nodes(::Val{:gauss_chebyshev}, n, T=Float64; k=1)
    x = zeros(T, n)
    if k == 1
        for i in eachindex(x)
            x[i] = sin( π*(n-2i+1) / (2n) )
        end
        return x[end:-1:1]
    elseif k == 2
        for i in eachindex(x)
            x[i] = cos( π*i / (n-1) )
        end
        return x
    end
end

function get_nodes(::Val{:lobatto_chebyshev}, n, T=Float64)
    x = zeros(T, n)
    for i in eachindex(x)
        x[i] = cos( (i-1)*π / (n-1) )
    end
    return x[end:-1:1]
end


text:
using DWARF
using DWARF.CallFrameInfo
using DWARF.CallFrameInfo: EhFrameRef, CIECache
using ObjFileBase
using ObjFileBase: handle, isrelocatable, Sections, mangle_sname
using ELF
using MachO

const InverseSymtab = Vector{UInt32}
const FDETab = Vector{Tuple{Int,UInt}}

type ConcreteModule{T<:ObjFileBase.ObjectHandle, SR<:ObjFileBase.SectionRef}
    base::UInt64
    sz::UInt
    handle::T
    # The .eh_frame (or .debug_frame) section when unwinding using DWARF
    eh_frame::Nullable{SR}
    # Either a DWARF or SEH header/unwind table pair
    ehfr::Nullable{EhFrameRef}
    xpdata::Nullable{XPUnwindRef}
    # If the object containing this module's DWARF info is different,
    # this is the handle to it
    dwarfhandle::T
    # Keeps symbol table indices in ip order as opposed to alphabetically,
    # to accelerate symbol lookup by ip.
    inverse_symtab::Nullable{InverseSymtab}
    FDETab::Nullable{FDETab}
    ciecache::Nullable{CIECache}
    is_jit_dobj::Bool
    is_eh_not_debug::Bool
end

"""
    Entry point for users of this package to provide modules that are not backed
    by object files.
"""
type SyntheticModule
    base::UInt64
    sz::UInt
    # Note: ip is guaranteed to be in this module, for the next two functions
    # (session, RC) -> RC
    unwind_step
    # (session, ip) -> (Bool, symb) - see symbolicate()
    symbolicate
    # (session, ip) -> Range
    get_proc_bounds
end

const Module = Union{ConcreteModule,SyntheticModule}

"""
A way to obtain new modules, e.g. dynamic libraries or jit modules
"""
abstract ModuleSource

immutable ModuleSet{intptrT}
    sources::Vector{ModuleSource}
    modules::Dict{RemotePtr{Void, intptrT},Module}
end
(::ModuleSet{intptrT}){intptrT}(sources::Vector{ModuleSource},x::Dict) =
  ModuleSet{intptrT}(sources, Dict{RemotePtr{Void, intptrT},Module}(x))

type MultiASModules{F,T,S}
    new_as_callback
    modules_by_as::Dict{T, S}  # ASID => modules
end
MultiASModules{F,T,S}(f::F,modules_by_as::Dict{T,S}) =
  MultiASModules{F,T,S}(f, modules_by_as)
function current_asid
end

ObjFileBase.handle(mod::Module) = mod.handle
dhandle(mod::Module) = mod.dwarfhandle
dhandle(h) = h
type LazyJITModules{T}
    modules::T
    nsharedlibs::Int
end
LazyJITModules() = LazyJITModules(Dict{RemotePtr{Void}, Any}(), 0)

function first_actual_segment(h)
    deref(first(Iterators.filter(x->isa(deref(x),MachO.segment_commands)&&(segname(x)!="__PAGEZERO"), LoadCmds(handle(h)))))
end

function get_syms(mod)
    local syms
    h = dhandle(mod)
    sections = Sections(h)
    if isa(h, ELF.ELFHandle)
        secs = collect(Iterators.filter(x->sectionname(x) == ".symtab",sections))
        isempty(secs) && (secs = collect(Iterators.filter(x->sectionname(x) == ".dynsym",sections)))
        syms = ELF.Symbols(secs[1])
    elseif isa(h, MachO.MachOHandle)
        syms = MachO.Symbols(h)
    elseif isa(h, COFF.COFFHandle)
        syms = COFF.Symbols(h)
    end
    syms
end

function make_inverse_symtab(h)
    sects = collect(Sections(h))
    first_ph_vaddr = isa(h, ELF.ELFHandle) ? ELF.compute_first_ph_vaddr(h) : 0
    UInt32[symbolnum(x) for x in sort!(collect(get_syms(h)), by = x->symbolvalue(x,sects,first_ph_vaddr))]
end

function compute_mod_size(h)
    if isrelocatable(handle(h))
        sz = sectionsize(first(Iterators.filter(x->sectionname(x) == mangle_sname(handle(h),"text"), Sections(handle(h)))))
    elseif isa(handle(h), ELF.ELFHandle)
        phs = ELF.ProgramHeaders(handle(h))
        sz = first(Iterators.filter(x->x.p_type == ELF.PT_LOAD, phs)).p_memsz
    elseif isa(handle(h), MachO.MachOHandle)
        sz = first_actual_segment(h).vmsize
    elseif isa(handle(h), COFF.COFFHandle)
        sz = COFF.readoptheader(handle(h)).standard.SizeOfCode
    end
    UInt(sz)
end
compute_mod_size(m::Module) = UInt(m.sz)

@inline function _find_module(modules, theip)
    ip = UInt(theip)
    for (mbase, h) in modules
        base = UInt(mbase)
        sz = compute_mod_size(h)::UInt
        if base <= ip <= base + sz
            return Nullable{Pair{UInt,Any}}(Pair{UInt,Any}(base,h))
        end
    end
    Nullable{Pair{UInt,Any}}()
end

function find_module(session, modules, ip)
    ret = _find_module(module_dict(modules), ip)
    isnull(ret) && error("ip 0x$(hex(UInt(ip),2sizeof(UInt))) not found")
    tuple(get(ret)...)
end

"""
base is the remote address of the text section
"""
function make_fdetab(base, mod, is_eh_not_debug=isa(mod, Module)?mod.is_eh_not_debug:true)
    FDETab = Vector{Tuple{Int,UInt}}()
    eh_frame = find_ehframes(mod)[1]
    for fde in FDEIterator(eh_frame,is_eh_not_debug,ObjFileBase.intptr(handle(mod)))
        # Make this location, relative to base. Note that base means something different,
        # depending on whether we're isrelocatable or not. If we are, base is the load address
        # of the text section. Otherwise, base is the load address of the start of the mapping.
        if isa(handle(mod), ELF.ELFHandle)
            # For relocated ELF Objects and unlrelocated shared libraries,
            # it is an offset from the load address of the FDE (in the shared library case,
            # sh_addr == sectionoffset).
            ip = (Int(deref(eh_frame).sh_addr) - Int(sectionoffset(eh_frame)) +
                Int(initial_loc(fde))) - Int(UInt(base))
        elseif isa(handle(mod), COFF.COFFHandle)
            ip = initial_loc(fde) - COFF.readoptheader(handle(mod)).windows.ImageBase
        elseif isrelocatable(handle(mod))
            # MachO eh_frame section doesn't get relocated, so it's still relative to
            # the file's local address space.
            text = first(Iterators.filter(x->sectionname(x)==mangle_sname(mod,"text"),Sections(mod)))
            ip = initial_loc(fde) - sectionoffset(text)
        else
            ip = initial_loc(fde)
        end
        push!(FDETab,(ip, fde.offset))
    end
    sort!(FDETab, by=x->x[1])
    FDETab
end

function obtain_uuid(h)
    deref(first(Iterators.filter(x->isa(deref(x),MachO.uuid_command), LoadCmds(h)))).uuid
end

@static if is_apple()
    function obtain_dsym(fname, objecth)
        dsympath = string(fname, ".dSYM/Contents/Resources/DWARF/", basename(fname))
        isfile(dsympath) || return objecth
        debugh = readmeta(IOBuffer(open(Base.Mmap.mmap, dsympath)))
        (obtain_uuid(objecth) != obtain_uuid(debugh)) && return debugh
        debugh
    end
end

filter_named_sections(sects::ObjFileBase.Sections, name) = Iterators.filter(x->!isBSS(x) &&
      (sectionname(x)==mangle_sname(handle(sects),name)), sects)
filter_named_sections(h, name) = filter_named_sections(Sections(handle(h)), name)

find_ehframes{T<:Union{ELF.ELFHandle,MachO.MachOHandle}}(sects::ObjFileBase.Sections{T}) =
    collect(filter_named_sections(sects, "eh_frame"))
find_ehframes{T<:COFF.COFFHandle}(sects::ObjFileBase.Sections{T}) =
    collect(filter_named_sections(sects, "debug_frame"))
find_ehframes(h::ObjFileBase.ObjectHandle) = find_ehframes(Sections(h))
function find_ehframes(h::MachO.MachOHandle)
    mapreduce(x->find_ehframes(Sections(x)), vcat,
        Iterators.filter(x->isa(deref(x), MachO.segment_commands), LoadCmds(h)))
end
function find_ehframes{T}(m::ConcreteModule{T})
    (get(m.eh_frame)::SectionRef(T),)
end
_find_eh_frame_hdr(h) = filter_named_sections(h, "eh_frame_hdr")
function find_eh_frame_hdr(h)
    first(_find_eh_frame_hdr(h))
end
find_eh_frame_hdr{T}(m::ConcreteModule{T}) = (get(mod.ehfr).hdr_sec)::SectionRef(T)

find_ehfr(mod::ConcreteModule) = get(mod.ehfr)
find_ehfr(h) = EhFrameRef(find_eh_frame_hdr(h), find_ehframes(h)[1])

@static if is_apple()
    function update_module!(session::LocalSession, modules, base, fname)
        h = readmeta(IOBuffer(open(Base.Mmap.mmap, fname)))
        if isa(h, MachO.FatMachOHandle)
            h = h[findfirst(arch->arch.cputype == MachO.CPU_TYPE_X86_64, h.archs)]
        end
        # Do not record the dynamic linker in our module list
        # Also skip the executable for now
        ft = readheader(h).filetype
        if ft == MachO.MH_DYLINKER
            return false
        end
        # For now, don't record shared libraries for which we only
        # have compact unwind info.
        ehfs = find_ehframes(h)
        length(collect(ehfs)) == 0 && return false
        vmaddr = first_actual_segment(h).vmaddr
        base += vmaddr
        modules.modules[RemotePtr{Void}(base)] =
            ConcreteModule(base, compute_mod_size(h), h, Nullable(ehfs[]),
                Nullable{EhFrameRef}(),
                Nullable{XPUnwindRef}(),
                obtain_dsym(fname, h), Nullable{InverseSymtab}(),
                Nullable{FDETab}(), Nullable{CIECache}(),
                false, true)
        return true
    end

    function update_shlibs!(session::LocalSession, modules::LazyJITModules)
        nactuallibs = ccall(:_dyld_image_count, UInt32, ())
        if modules.nsharedlibs != nactuallibs
            for idx in (modules.nsharedlibs+1):nactuallibs
                idx -= 1
                base = ccall(:_dyld_get_image_vmaddr_slide, UInt, (UInt32,), idx)
                haskey(modules.modules, UInt64(base)) && continue
                fname = unsafe_string(
                    ccall(:_dyld_get_image_name, Ptr{UInt8}, (UInt32,), idx))
                # hooking is weird
                contains(fname, "hooking.dylib") &&
                    continue
                update_module!(session, modules, base, fname)
            end
            modules.nsharedlibs = nactuallibs
            return true
        end
        return false
    end
elseif is_windows()
    function update_shlibs!(session, modules)
        return false
    end
elseif is_linux()
    immutable dl_phdr_info
        # Base address of object
        addr::Ptr{Void}

        # Null-terminated name of object
        name::Ptr{UInt8}

        # Pointer to array of ELF program headers for this object
        phdr::Ptr{Void}

        # Number of program headers for this object
        phnum::Cshort
    end

    const LibArray = Array{Tuple{String,Ptr{Void}},1}

    # This callback function called by dl_iterate_phdr() on Linux
    function dl_phdr_info_callback(di::dl_phdr_info, size::Csize_t, dynamic_libraries::LibArray)
        # Skip over objects without a path (as they represent this own object)
        name = unsafe_string(di.name)
        push!(dynamic_libraries, (name, di.addr))
        return convert(Cint, 0)::Cint
    end

    function update_module!(session::LocalSession, modules, base, fname)
        is_exe = isempty(fname)
        if is_exe
            fname = readlink("/proc/self/exe")
        end
        # hooking is weird
        contains(fname, "hooking.so") &&
            return false
        # linux-vdso is fake
        contains(fname, "linux-vdso") &&
            return false
        h = readmeta(IOBuffer(open(Base.Mmap.mmap,fname)))
        if is_exe
            phs = ELF.ProgramHeaders(h)
            idx = findfirst(p->p.p_offset==0&&p.p_type==ELF.PT_LOAD, phs)
            base += phs[idx].p_vaddr
        end
        base = RemotePtr{Void}(UInt64(base))
        modules.modules[base] = GlibcDyldModules.mod_for_h(h, base, fname)
        return true
    end

    function update_shlibs!(session::LocalSession, modules::LazyJITModules)
        dynamic_libraries = LibArray(0)

        @static if is_linux()
            const callback = cfunction(dl_phdr_info_callback, Cint,
                                       (Ref{dl_phdr_info}, Csize_t, Ref{LibArray} ))
            ccall(:dl_iterate_phdr, Cint, (Ptr{Void}, Ref{LibArray}), callback, dynamic_libraries)
        end

        for (fname,base) in dynamic_libraries[modules.nsharedlibs+1:end]
            haskey(modules.modules, UInt64(base)) && continue
            update_module!(session, modules, base, fname)
        end

        did_update = length(dynamic_libraries) > modules.nsharedlibs
        modules.nsharedlibs = length(dynamic_libraries)
        did_update
    end
end

function jit_mod_for_h(buf, h, sstart)
    sstart = UInt(sstart)
    fdetab = Vector{Tuple{Int,UInt}}()
    ehfr = Nullable{EhFrameRef}()
    xpdata = Nullable{XPUnwindRef}()
    ciecache = CIECache()
    if isrelocatable(h)
      isa(h, ELF.ELFHandle) && ELF.relocate!(buf, h)
      fdetab = make_fdetab(sstart, h, true)
      if isa(h, MachO.MachOHandle)
        LOI = Dict(:__text => sstart,
            :__debug_str=>0) #This one really shouldn't be necessary
        MachO.relocate!(buf, h; LOI=LOI)
      end
    elseif isa(h, ELF.ELFHandle)
        ehfr = Nullable{EhFrameRef}(find_ehfr(h))
    elseif isa(h, COFF.COFFHandle)
        sects = Sections(h)
        pdata = collect(Iterators.filter(x->sectionname(x)==ObjFileBase.mangle_sname(h,"pdata"),sects))[]
        xdata = collect(Iterators.filter(x->sectionname(x)==ObjFileBase.mangle_sname(h,"xdata"),sects))[]
        xpdata = Nullable{XPUnwindRef}(XPUnwindRef(xdata, pdata))
    end
    isa(h, MachO.MachOHandle) && isempty(fdetab) && (fdetab = make_fdetab(sstart, h, true))
    eh_frame = find_ehframes(h)[]
    ConcreteModule(sstart, compute_mod_size(h), h, Nullable(eh_frame),
        ehfr, xpdata, h, Nullable{InverseSymtab}(),
        Nullable(fdetab), Nullable{CIECache}(), true, true)
end

function retrieve_obj_data(s::LocalSession, modules, ip)
    data = ccall(:jl_get_dobj_data, Any, (UInt,), ip)
    @assert data != nothing
    copy(data)
end
retrieve_section_start(s::LocalSession, modules, ip) =
    ccall(:jl_get_section_start, UInt64, (UInt,), ip)
update_shlibs!(session, modules::LazyJITModules) =
    update_shlibs!(session, modules.modules)
update_shlibs!(session, modules::Dict) = false

function update_shlibs!(session, modules::ModuleSet)
    did_change = false
    for source in modules.sources
        did_change |= update_shlibs!(session, modules, source)
    end
    return did_change
end

function find_module(session, modules::LazyJITModules, ip)
    ret = _find_module(module_dict(session, modules), ip)
    return !isnull(ret) ? tuple(get(ret)...) : begin
        if update_shlibs!(session, modules)
            return find_module(session, modules, ip)
        end
        sstart = RemotePtr{Void}(UInt(retrieve_section_start(session, modules, ip-1)))
        sstart == 0 && error("ip 0x$(hex(UInt(ip),2sizeof(UInt))) not found")
        buf = IOBuffer(retrieve_obj_data(session, modules, ip), true, true)
        h = readmeta(buf)
        module_dict(modules)[sstart] = jit_mod_for_h(buf, h, sstart)
        sstart, module_dict(session, modules)[sstart]
    end
end
find_module(modules, ip) = find_module(LocalSession(), modules, ip)

function find_module{F,T,S}(session, modules::MultiASModules{F,T,S}, ip)
    asid = current_asid(session)
    if !haskey(modules.modules_by_as, asid)
        modules.modules_by_as[asid] = modules.new_as_callback(session)::S
    end
    find_module(session, modules.modules_by_as[asid], ip)
end

module_dict(modules) = module_dict(LocalSession(), modules)
module_dict(session, modules::ModuleSet) = modules.modules
module_dict(session, modules) = modules
module_dict(session, modules::LazyJITModules) = module_dict(session, modules.modules)
module_dict(session, modules::MultiASModules) =
    module_dict(session, modules.modules_by_as[current_asid(session)])

function lookup_sym(session, modules, name)
    ret = lookup_syms(session, modules, name)
    length(ret) == 0 && error("Not found")
    ret[]
end

function lookup_syms(session, modules, name, n = typemax(UInt))
    ret = Any[]
    name = string(name)
    for (base, h) in module_dict(session, modules)
      isa(h, SyntheticModule) && continue
      symtab = get_syms(handle(h))
      strtab = ObjFileBase.StrTab(symtab)
      idx = findfirst(x->ObjFileBase.symname(x, strtab = strtab)==name,symtab)
      if idx != 0
          sym = symtab[idx]
          if ObjFileBase.isundef(sym)
              continue
          end
          push!(ret,(h, base, sym))
          length(ret) >= n && return ret
      end
    end
    ret
end

function first_executable_segment(phs)
    idx = findfirst(p->p.p_type==ELF.PT_LOAD &&
                       ((p.p_flags & ELF.PF_X) != 0), phs)
    phs[idx]
end

"""
Load the set of active modules from the GlibC dynamic linker
"""
module GlibcDyldModules
  using ELF
  using NativeDebugger
  using ObjFileBase
  using NativeDebugger: load, mapped_file, XPUnwindRef, InverseSymtab, CIECache,
    FDETab, ModuleSource, ModuleSet, module_dict
  using NativeDebugger.Registers: intptr, getarch
  import NativeDebugger: update_shlibs!
  using DWARF.CallFrameInfo
  using DWARF.CallFrameInfo: EhFrameRef
  using ObjFileBase: mangle_sname, handle

  immutable link_map{intptrT}
    l_addr::RemotePtr{Void,intptrT}
    l_name::RemotePtr{UInt8,intptrT}
    l_ld::RemotePtr{Void,intptrT}
    l_next::RemotePtr{link_map{intptrT},intptrT}
    l_prev::RemotePtr{link_map{intptrT},intptrT}
  end

  immutable r_debug{intptrT}
    r_version::Cint
    link_map::RemotePtr{link_map{intptrT},intptrT}
    r_brk::intptrT
    r_state::Cint
    r_ldbase::RemotePtr{Void,intptrT}
  end

  """
    Computes the object's entry point from the kernel's auxv data. By comparing
    this information to the specified entrypoint in the executable image, one
    can compute the slide of the executable image in memory.
  """
  function compute_entry_ptr(sess, auxv_data::Vector{UInt8})
    auxv = reinterpret(intptr(getarch(sess)), auxv_data)
    entry_idx = findfirst(i->auxv[i] == ELF.AT_ENTRY, 1:2:length(auxv))
    @assert entry_idx != 0
    #  auxv_idx = 1+2(entry_idx-1), we want auxv_idx + 1 == entry_idx
    entry_ptr = auxv[2entry_idx]
  end

  using CRC
  const crc32 = crc(CRC_32)

  """
    Try to obtain a handle to this object's debug object (a second object that
    contains the separated-out debug information). This is done by searching
    following any .gnu_debuglink section that this object may contain.

    References:
    https://access.redhat.com/documentation/en-US/Red_Hat_Enterprise_Linux/4/html/Debugging_with_gdb/separate-debug-files.html
    https://blogs.oracle.com/dbx/entry/gnu_debuglink_or_debugging_system
  """
  function search_debug_object(h, filename)
      sects = Sections(h)
      debuglinks = collect(Iterators.filter(x->sectionname(x)==mangle_sname(handle(sects),"gnu_debuglink"),sects))
      isempty(debuglinks) && return h
      debuglink = debuglinks[]
      # Read the null terminated name of the debug library
      seek(debuglink, 0)
      fname = readuntil(handle(debuglink), '\0')[1:end-1]
      # Align to 4 byte boundary
      skip(handle(debuglink), ((4-position(debuglink)%4)%4))
      # Read checksum
      crc = read(debuglink, UInt32)
      # Now search for this file in the same locations as GDB
      execdir = dirname(filename)
      const globaldir = "/usr/lib/debug"
      for path in [execdir, joinpath(execdir,".debug"),
                   # For /usr/lib/..., /usr/lib/debug/usr/lib/...
                   joinpath(globaldir, execdir[2:end])]
         fp = joinpath(path, fname)
         isfile(fp) || continue
         buf = open(Base.Mmap.mmap, fp)
         crc == crc32(buf) || continue
         return readmeta(IOBuffer(buf))
      end
      return h
  end

  function mod_for_h(h, base, filename)
      eh_frame = NativeDebugger.find_ehframes(h)[]
      eh_frame_hdrs = collect(NativeDebugger._find_eh_frame_hdr(h))
      ehfr = isempty(eh_frame_hdrs) ? Nullable{EhFrameRef}() :
        Nullable{EhFrameRef}(EhFrameRef(eh_frame_hdrs[], eh_frame))
      dh = search_debug_object(h, filename)
      NativeDebugger.ConcreteModule(UInt64(base),
        NativeDebugger.compute_mod_size(h), h, Nullable(eh_frame),
        ehfr, Nullable{XPUnwindRef}(), dh,
        Nullable{InverseSymtab}(),
        Nullable{FDETab}(),
        Nullable{CIECache}(),
        false, true)
  end

  immutable GlibCRemoteSource <: ModuleSource
     dt_debug_addr_addr::RemotePtr
  end

  function update_shlibs!(vm, modules, source::GlibCRemoteSource; current_ip = 0)
      local debug_struct, last_link_map
      #try
          dt_debug_addr = load(vm, source.dt_debug_addr_addr)
          if dt_debug_addr == 0
              # Dynamic linker not yet set up, if we know the current ip, we can
              # try get the library by looking at the memory mapping for this address
              if current_ip != 0
                  base = NativeDebugger.segment_base(vm, current_ip)
                  fn = mapped_file(vm, base)
                  if !isempty(fn)
                      # Don't use the IOStream directly. We do a look of seeking/poking,
                      # so loading the whole thing into memory and using an IOBuffer is
                      # faster
                      buf = IOBuffer(open(Base.Mmap.mmap, fn))
                      h = readmeta(buf)
                      module_dict(modules)[base] = mod_for_h(h, base, fn)
                      return true
                  end
              end
              return false
          end
          debug_struct = load(vm, dt_debug_addr)
          last_link_map = load(vm, debug_struct.link_map)
      #catch err
      #      @show err
    #      return
      #end

      # Now for the shared libraries. We traverse the linked list of modules
      # in the dynamic linker's r_debug structure.
      lm = last_link_map
      did_update = false
      while lm.l_next.ptr != 0
          # Do not add to the module list if the name is empty. This is true for
          # the main executable as well as the dynamic library loader
          if !haskey(modules.modules, lm.l_addr) && load(vm, lm.l_name) != 0
              # name = load(vm, lm.l_name, 255)
              # idx = findfirst(name,0)
              # idx == 0 ? length(name) : idx
              # @show String(name[1:idx])
              # @show lm.l_addr
              # Some libraries (e.g. wine's fake windows libraries) show up with
              # load address 0, ignore them
              if lm.l_addr != 0
                  fn = mapped_file(vm, lm.l_addr)
                  if !isempty(fn)
                      # The vdso doesn't exist on disk. Load it from the remote
                      # address space.
                      if startswith(fn, "linux-vdso.so.1")
                          # 2 pages is the current size of the vdso. If it ever
                          # grows larger, we may have to ask for the actual
                          # mapping size here
                          buf = IOBuffer(load(vm, RemotePtr{UInt8}(lm.l_addr), 2*0x1000))
                      else
                          # Don't use the IOStream directly. We do a lot of seeking/poking,
                          # so loading the whole thing into memory and using an IOBuffer is
                          # faster
                          buf = IOBuffer(open(Base.Mmap.mmap, fn))
                      end
                      h = readmeta(buf)
                      module_dict(modules)[lm.l_addr] = mod_for_h(h, lm.l_addr, fn)
                      did_update = true
                  end
              end
          end
          lm = load(vm, lm.l_next)
      end
      return did_update
  end

  """
    Load the shared library map from address space `vm`.
    `imageh` should be the ObjectHandle of the main executable and `auxv_data`
    should be the session's auxv buffer.

    If the executable was not loaded at it's defined virtual address, you
    may set `image_slide` to the offset of the executable's load address from
    it's intended load address.
  """
  function load_library_map(vm, imageh, image_slide = 0; current_ip = 0)
    # Construct the actual module map
    modules = Dict{RemotePtr{Void},Any}()

    # First the main executable. To do so we need to find the image base.
    phs = ELF.ProgramHeaders(imageh)
    # We want the first loaded segment that's executable
    imagebase = NativeDebugger.first_executable_segment(phs).p_vaddr + image_slide
    modules[RemotePtr{Void}(imagebase)] = mod_for_h(imageh, imagebase, "")

    intptr_t = isa(imageh.file, ELF.ELF64.File) ? UInt64 : UInt32

    # Obtain the target address space's r_debug struct
    secs = collect(Iterators.filter(x->sectionname(x)==".dynamic",ELF.Sections(imageh)))
    isempty(secs) && (return ModuleSet(ModuleSource[], modules))
    dynamic_sec = secs[]
    dynamic = reinterpret(intptr_t,read(dynamic_sec))
    dt_debug_idx = findfirst(i->dynamic[i]==ELF.DT_DEBUG, 1:2:length(dynamic))
    @assert dt_debug_idx != 0

    dynamic_load_addr = deref(dynamic_sec).sh_addr + image_slide
    dt_debug_addr_addr = RemotePtr{RemotePtr{r_debug{intptr_t},intptr_t},intptr_t}(
        dynamic_load_addr + (2*dt_debug_idx-1)*sizeof(intptr_t))
    source = GlibCRemoteSource(dt_debug_addr_addr)
    modules = ModuleSet{intptr_t}(ModuleSource[source], modules)
    update_shlibs!(vm, modules, source; current_ip=current_ip)

    modules
  end

end

"""
Load the set of active modules from the Windows dynamic linker
"""
module WinDyldModules
    using COFF
    using ELF
    using NativeDebugger
    using ObjFileBase
    using NativeDebugger: load, mapped_file, make_fdetab, XPUnwindRef, InverseSymtab,
        ModuleSource, ModuleSet, module_dict
    import NativeDebugger: update_shlibs!
    using DWARF.CallFrameInfo
    using DWARF.CallFrameInfo: EhFrameRef, CIECache
    using ..GlibcDyldModules
    import Base: start, next, done

    # Note: These are the wine implementations
    function get_peb_addr(::NativeDebugger.X86_64.X86_64Arch, vm)
        # Get TEB
        teb_ptr = NativeDebugger.get_dwarf(NativeDebugger.getregs(vm), :gs_base)
        peb_ptr_ptr = RemotePtr{UInt64}(teb_ptr + 0x60)
        load(vm, peb_ptr_ptr)
    end

    function get_peb_addr(::NativeDebugger.X86_32.X86_32Arch, vm)
        # Get TEB
        fs_segment = NativeDebugger.get_dwarf(NativeDebugger.getregs(vm), :fs)
        ((fs_segment & 0x4) == 0) || error("TLS segemnt is in LDT, need to add support")
        gdt_entry = fs_segment >> 3
        teb_ptr = NativeDebugger.get_thread_area_base(vm, gdt_entry)
        peb_ptr_ptr = RemotePtr{UInt32}(teb_ptr + 0x30)
        load(vm, peb_ptr_ptr)
    end

    get_peb_addr(vm) = get_peb_addr(NativeDebugger.getarch(vm), vm)

    # See https://msdn.microsoft.com/en-us/library/windows/desktop/aa380518(v=vs.85).aspx
    immutable UNICODE_STRING
        Length::UInt16
        MaximumLength::UInt16
        Buffer::Ptr{UInt16}
    end

    # See https://msdn.microsoft.com/en-us/library/windows/desktop/aa813706(v=vs.85).aspx
    immutable LDR_DATA_TABLE_ENTRY{ptrT}
        Reserved::NTuple{2,RemotePtr{Void,ptrT}}
        Prev::RemotePtr{LDR_DATA_TABLE_ENTRY{ptrT},ptrT}
        Next::RemotePtr{LDR_DATA_TABLE_ENTRY{ptrT},ptrT}
        Reserved2::NTuple{2,RemotePtr{Void,ptrT}}
        DllBase::RemotePtr{Void,ptrT}
        EntryPoint::RemotePtr{Void,ptrT}
        Reserved3::RemotePtr{Void,ptrT}
        FullDllName::UNICODE_STRING
    end

    immutable ModuleIterator
        vm
        head::RemotePtr
    end

    function mod_for_h(dllbase, h)
        sects = Sections(h)
        pdatas = collect(Iterators.filter(x->sectionname(x)==ObjFileBase.mangle_sname(h,"pdata"),sects))
        xdatas = collect(Iterators.filter(x->sectionname(x)==ObjFileBase.mangle_sname(h,"xdata"),sects))
        xpdata = Nullable{XPUnwindRef}()
        ehframe = Nullable{typeof(first(sects))}()
        is_eh_not_debug = true
        if !isempty(pdatas) && !isempty(xdatas)
            xpdata = Nullable{XPUnwindRef}(XPUnwindRef(xdatas[], pdatas[]))
        else
            # Try to find any eh_frame or debug_frame sections (e.g. MinGW compilers)
            eh_frames = collect(Iterators.filter(x->sectionname(x)==ObjFileBase.mangle_sname(h,"eh_frame"),sects))
            debug_frames = collect(Iterators.filter(x->sectionname(x)==ObjFileBase.mangle_sname(h,"debug_frame"),sects))
            ehframe = !isempty(eh_frames) ? Nullable(eh_frames[]) :
                      !isempty(debug_frames) ? Nullable(debug_frames[]) :
                      ehframe
            is_eh_not_debug = !isempty(eh_frames)
        end
        NativeDebugger.ConcreteModule(UInt64(dllbase),
            NativeDebugger.compute_mod_size(h),
            h, ehframe, Nullable{EhFrameRef}(),
            xpdata, h,
            Nullable{InverseSymtab}(),
            Nullable{Vector{Tuple{Int,UInt}}}(),
            Nullable{CIECache}(), false, is_eh_not_debug)
    end

    const PEB_LDR_DATA_IDX      = 3
    InMemoryOrderList_OFFSET(ptrT) = 8 + 3*sizeof(ptrT)
    function ModuleIterator(vm)
        addr = get_peb_addr(vm)
        ptrT = NativeDebugger.intptr(NativeDebugger.getarch(vm))
        peb_ldr_data_addr = load(vm,
            RemotePtr{ptrT}(addr + PEB_LDR_DATA_IDX*sizeof(ptrT)))
        head = peb_ldr_data_addr + InMemoryOrderList_OFFSET(ptrT)
        ModuleIterator(vm, RemotePtr{LDR_DATA_TABLE_ENTRY{ptrT},ptrT}(head))
    end

    function start(m::ModuleIterator)
        ptrT = typeof(m.head).parameters[2]
        load(m.vm, RemotePtr{RemotePtr{LDR_DATA_TABLE_ENTRY{ptrT},ptrT},ptrT}(m.head))
    end
    function next{T}(m::ModuleIterator, s::RemotePtr{LDR_DATA_TABLE_ENTRY{T}})
        entry = load(m.vm, s-2*sizeof(T))
        (entry, entry.Prev)
    end
    done{T}(m::ModuleIterator, s::RemotePtr{LDR_DATA_TABLE_ENTRY{T}}) =
        s == m.head || s == 0
    Base.iteratorsize(::Type{ModuleIterator}) = Base.SizeUnknown()

    immutable WinRemoteSource <: ModuleSource
    end

    function update_shlibs!(vm, modules, source::WinRemoteSource)
        did_change = false
        map(ModuleIterator(vm)) do entry
            entry.DllBase == 0 && return
            haskey(module_dict(modules), entry.DllBase) && return
            # If this is a Wine library, the first 0x1000 are dynamically
            # allocated to create space for the COFF header. Also try to
            # look 0x1000 bytes after to make sure to get the right file.
            fn = mapped_file(vm, entry.DllBase)
            if isempty(fn)
                fn = mapped_file(vm, entry.DllBase + 0x1000)
                if !isempty(fn)
                    base = NativeDebugger.segment_base(vm, entry.DllBase-0x1000)
                    haskey(module_dict(modules), base) && return
                end
            end
            if !isempty(fn)
                h = readmeta(IOBuffer(open(Base.Mmap.mmap, fn)))
                # If this is an ELF library, it is likely one of wine's internal
                # fake DLLs, follow the Linux codepath instead
                if isa(h, ELF.ELFHandle)
                    # We need to find the actual ELF base address for this library
                    # Wine rounds the location of the PE header up to the nearset
                    # 64K. We can't reverse that computiation, we we try to look
                    # at the base address of the previous page and hope that nobody
                    # messsed with the memory layout.
                    base = NativeDebugger.segment_base(vm, entry.DllBase-0x1000)
                    haskey(module_dict(modules), base) && return
                    module_dict(modules)[base] =
                        GlibcDyldModules.mod_for_h(h, entry.DllBase, fn)
                else
                    module_dict(modules)[entry.DllBase] = mod_for_h(entry.DllBase, h)
                end
                did_change = true
            end
        end
        did_change
    end

    function load_library_map(vm)
        modules = ModuleSet(ModuleSource[WinRemoteSource()], Dict{RemotePtr{Void},Any}())
        update_shlibs!(vm, modules)
        modules
    end
end

# Module to represent the linux kernel
immutable LinuxKernelModule
    # Kernel Symbol table
    addrs::Vector{UInt64}
    names::Vector{String}
    mapped_size::UInt64
    function LinuxKernelModule(mapped_size)
        addrs = UInt64[]
        names = String[]
        for line in eachline("/proc/kallsyms")
            pieces = split(line,' ')
            push!(addrs, parse(UInt64,pieces[1],16))
            push!(names, pieces[3][1:end-1])
        end
        new(addrs, names, mapped_size)
    end
end
compute_mod_size(mod::LinuxKernelModule) = mod.mapped_size

function symbolicate(sess, base, kern::LinuxKernelModule, ip)
    id = searchsortedlast(kern.addrs, ip)
    (true, kern.names[id])
end


text:
# generate examples
import Literate
import ClosedLoopReachability: @modelpath

MODELS = [
    joinpath(@__DIR__, "..", "models", "ACC"),
    joinpath(@__DIR__, "..", "models", "Sherlock-Benchmark-9-TORA"),
    joinpath(@__DIR__, "..", "models", "Sherlock-Benchmark-10-Unicycle"),
    joinpath(@__DIR__, "..", "models", "VertCAS"),
    joinpath(@__DIR__, "..", "models", "Single-Pendulum"),
    joinpath(@__DIR__, "..", "models", "Double-Pendulum"),
    joinpath(@__DIR__, "..", "models", "Airplane"),
    joinpath(@__DIR__, "..", "models", "Sherlock-Benchmark-7"),
    joinpath(@__DIR__, "..", "models", "Cart-Pole")
]
GENERATEDDIR = joinpath(@__DIR__, "src", "models")
MODELDIR = joinpath(@__DIR__, "..", "models")
mkpath(GENERATEDDIR)

macro modelpath(model_path::String, name::String)
    return joinpath(MODELDIR, model_path, name)
end

for model in MODELS
    for file in readdir(model)
        if endswith(file, ".jl")
            input = abspath(joinpath(model, file))
            script = Literate.script(input, GENERATEDDIR, credit=false)
            code = strip(read(script, String))
            mdpost(str) = replace(str, "@__CODE__" => code)
            Literate.markdown(input, GENERATEDDIR, postprocess=mdpost, credit=false)
            Literate.notebook(input, GENERATEDDIR, execute = true, credit = false)
        elseif any(endswith.(file, [".png", ".jpg", ".gif"]))
            cp(joinpath(model, file), joinpath(GENERATEDDIR, file); force=true)
        else
            @warn "ignoring $file in $model"
        end
    end
end


text:
module Pathogen

  # Dependencies
  using Distributed,
        DataFrames,
        Distributions,
        RecipesBase,
        Logging,
        StatsBase,
        Statistics,
        ProgressMeter,
        LinearAlgebra,
        OnlineStats

  import StatsBase.mean,
         StatsBase.mode,
         Base.summary,
         ProgressMeter.next!

  # Core
  include("core/EpidemicModel.jl")
  include("core/DiseaseState.jl")
  include("core/Population.jl")
  include("core/Events/Event.jl")
  include("core/Events/Events.jl")
  include("core/Events/EventRates.jl")
  include("core/Events/EventObservations.jl")
  include("core/Risks/AbstractRisk.jl")
  include("core/Risks/RiskFunctions.jl")
  include("core/Risks/RiskParameters.jl")
  include("core/Transmissions/Transmission.jl")
  include("core/Transmissions/AbstractTransmissionNetwork.jl")
  include("core/Transmissions/TransmissionNetwork.jl")
  include("core/Transmissions/TransmissionRates.jl")
  include("core/Transmissions/TransmissionRateCache.jl")
  include("core/initialize.jl")
  include("core/update!.jl")

  # Simulation
  include("simulation/Simulation.jl")
  include("simulation/generate.jl")
  include("simulation/simulate!.jl")
  include("simulation/update!.jl")
  include("simulation/observe.jl")

  # Inference
  include("inference/TransmissionNetworkDistribution.jl")
  include("inference/RiskPriors.jl")
  include("inference/EventExtents.jl")
  include("inference/MarkovChain.jl")
  include("inference/MCMC.jl")
  include("inference/_pathway_from.jl")
  include("inference/_pathway_to.jl")
  include("inference/_accept.jl")
  include("inference/_bounds.jl")
  include("inference/generate.jl")
  include("inference/initialize.jl")
  include("inference/update!.jl")
  include("inference/loglikelihood.jl")
  include("inference/logpriors.jl")
  include("inference/start!.jl")
  include("inference/iterate!.jl")
  include("inference/summary.jl")

  # Visualization
  include("visualization/epidemic_curve.jl")
  include("visualization/observations.jl")
  include("visualization/population.jl")
  include("visualization/transmission_network.jl")
  include("visualization/transmission_network_distribution.jl")
  include("visualization/degree_distribution.jl")
  include("visualization/trace.jl")

  export
    #re-exports from DataFrames.jl
    DataFrame,
    #re-exports from Distributions.jl
    UnivariateDistribution,
    Arcsine,
    Bernoulli,
    Beta,
    BetaBinomial,
    BetaPrime,
    Binomial,
    Biweight,
    Categorical,
    Cauchy,
    Chernoff,
    Chi,
    Chisq,
    Cosine,
    DiagNormal,
    DiagNormalCanon,
    Dirichlet,
    DirichletMultinomial,
    DiscreteUniform,
    DoubleExponential,
    EdgeworthMean,
    EdgeworthSum,
    EdgeworthZ,
    Erlang,
    Epanechnikov,
    Exponential,
    FDist,
    FisherNoncentralHypergeometric,
    Frechet,
    FullNormal,
    FullNormalCanon,
    Gamma,
    DiscreteNonParametric,
    GeneralizedPareto,
    GeneralizedExtremeValue,
    Geometric,
    Gumbel,
    Hypergeometric,
    InverseWishart,
    InverseGamma,
    InverseGaussian,
    IsoNormal,
    IsoNormalCanon,
    Kolmogorov,
    KSDist,
    KSOneSided,
    Laplace,
    Levy,
    LKJ,
    LocationScale,
    Logistic,
    LogNormal,
    LogitNormal,
    MatrixBeta,
    MatrixFDist,
    MatrixNormal,
    MatrixTDist,
    MixtureModel,
    Multinomial,
    MultivariateNormal,
    MvLogNormal,
    MvNormal,
    MvNormalCanon,
    MvNormalKnownCov,
    MvTDist,
    NegativeBinomial,
    NoncentralBeta,
    NoncentralChisq,
    NoncentralF,
    NoncentralHypergeometric,
    NoncentralT,
    Normal,
    NormalCanon,
    NormalInverseGaussian,
    Pareto,
    PGeneralizedGaussian,
    Product,
    Poisson,
    PoissonBinomial,
    QQPair,
    Rayleigh,
    Semicircle,
    Skellam,
    StudentizedRange,
    SymTriangularDist,
    TDist,
    TriangularDist,
    Triweight,
    Truncated,
    Uniform,
    UnivariateGMM,
    VonMises,
    VonMisesFisher,
    WalleniusNoncentralHypergeometric,
    Weibull,
    Wishart,
    ZeroMeanIsoNormal,
    ZeroMeanIsoNormalCanon,
    ZeroMeanDiagNormal,
    ZeroMeanDiagNormalCanon,
    ZeroMeanFullNormal,
    ZeroMeanFullNormalCanon,
    # Pathogen.jl types/functions
    EpidemicModel,
    SEIR, SEI, SIR, SI,
    DiseaseState,
    State_S, State_E, State_I, State_R,
    RiskFunctions, RiskParameters, RiskPriors,
    Population,
    AbstractTransmissionNetwork,
    TransmissionNetwork,
    TransmissionNetworkDistribution, TNDistribution,
    TransmissionNetworkPrior, TNPrior,
    TransmissionNetworkPosterior, TNPosterior,
    EndogenousTransmission, ExogenousTransmission, NoTransmission,
    Simulation,
    next!, simulate!,
    Events, EventObservations, EventExtents,
    observe,
    MCMC, start!, iterate!,
    summary, mean, mode
end


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

mutable struct IoK8sApiExtensionsV1beta1DaemonSetStatus <: SwaggerModel
    collisionCount::Any # spec type: Union{ Nothing, Int32 } # spec name: collisionCount
    conditions::Any # spec type: Union{ Nothing, Vector{IoK8sApiExtensionsV1beta1DaemonSetCondition} } # spec name: conditions
    currentNumberScheduled::Any # spec type: Union{ Nothing, Int32 } # spec name: currentNumberScheduled
    desiredNumberScheduled::Any # spec type: Union{ Nothing, Int32 } # spec name: desiredNumberScheduled
    numberAvailable::Any # spec type: Union{ Nothing, Int32 } # spec name: numberAvailable
    numberMisscheduled::Any # spec type: Union{ Nothing, Int32 } # spec name: numberMisscheduled
    numberReady::Any # spec type: Union{ Nothing, Int32 } # spec name: numberReady
    numberUnavailable::Any # spec type: Union{ Nothing, Int32 } # spec name: numberUnavailable
    observedGeneration::Any # spec type: Union{ Nothing, Int64 } # spec name: observedGeneration
    updatedNumberScheduled::Any # spec type: Union{ Nothing, Int32 } # spec name: updatedNumberScheduled

    function IoK8sApiExtensionsV1beta1DaemonSetStatus(;collisionCount=nothing, conditions=nothing, currentNumberScheduled=nothing, desiredNumberScheduled=nothing, numberAvailable=nothing, numberMisscheduled=nothing, numberReady=nothing, numberUnavailable=nothing, observedGeneration=nothing, updatedNumberScheduled=nothing)
        o = new()
        validate_property(IoK8sApiExtensionsV1beta1DaemonSetStatus, Symbol("collisionCount"), collisionCount)
        setfield!(o, Symbol("collisionCount"), collisionCount)
        validate_property(IoK8sApiExtensionsV1beta1DaemonSetStatus, Symbol("conditions"), conditions)
        setfield!(o, Symbol("conditions"), conditions)
        validate_property(IoK8sApiExtensionsV1beta1DaemonSetStatus, Symbol("currentNumberScheduled"), currentNumberScheduled)
        setfield!(o, Symbol("currentNumberScheduled"), currentNumberScheduled)
        validate_property(IoK8sApiExtensionsV1beta1DaemonSetStatus, Symbol("desiredNumberScheduled"), desiredNumberScheduled)
        setfield!(o, Symbol("desiredNumberScheduled"), desiredNumberScheduled)
        validate_property(IoK8sApiExtensionsV1beta1DaemonSetStatus, Symbol("numberAvailable"), numberAvailable)
        setfield!(o, Symbol("numberAvailable"), numberAvailable)
        validate_property(IoK8sApiExtensionsV1beta1DaemonSetStatus, Symbol("numberMisscheduled"), numberMisscheduled)
        setfield!(o, Symbol("numberMisscheduled"), numberMisscheduled)
        validate_property(IoK8sApiExtensionsV1beta1DaemonSetStatus, Symbol("numberReady"), numberReady)
        setfield!(o, Symbol("numberReady"), numberReady)
        validate_property(IoK8sApiExtensionsV1beta1DaemonSetStatus, Symbol("numberUnavailable"), numberUnavailable)
        setfield!(o, Symbol("numberUnavailable"), numberUnavailable)
        validate_property(IoK8sApiExtensionsV1beta1DaemonSetStatus, Symbol("observedGeneration"), observedGeneration)
        setfield!(o, Symbol("observedGeneration"), observedGeneration)
        validate_property(IoK8sApiExtensionsV1beta1DaemonSetStatus, Symbol("updatedNumberScheduled"), updatedNumberScheduled)
        setfield!(o, Symbol("updatedNumberScheduled"), updatedNumberScheduled)
        o
    end
end # type IoK8sApiExtensionsV1beta1DaemonSetStatus

const _property_map_IoK8sApiExtensionsV1beta1DaemonSetStatus = Dict{Symbol,Symbol}(Symbol("collisionCount")=>Symbol("collisionCount"), Symbol("conditions")=>Symbol("conditions"), Symbol("currentNumberScheduled")=>Symbol("currentNumberScheduled"), Symbol("desiredNumberScheduled")=>Symbol("desiredNumberScheduled"), Symbol("numberAvailable")=>Symbol("numberAvailable"), Symbol("numberMisscheduled")=>Symbol("numberMisscheduled"), Symbol("numberReady")=>Symbol("numberReady"), Symbol("numberUnavailable")=>Symbol("numberUnavailable"), Symbol("observedGeneration")=>Symbol("observedGeneration"), Symbol("updatedNumberScheduled")=>Symbol("updatedNumberScheduled"))
const _property_types_IoK8sApiExtensionsV1beta1DaemonSetStatus = Dict{Symbol,String}(Symbol("collisionCount")=>"Int32", Symbol("conditions")=>"Vector{IoK8sApiExtensionsV1beta1DaemonSetCondition}", Symbol("currentNumberScheduled")=>"Int32", Symbol("desiredNumberScheduled")=>"Int32", Symbol("numberAvailable")=>"Int32", Symbol("numberMisscheduled")=>"Int32", Symbol("numberReady")=>"Int32", Symbol("numberUnavailable")=>"Int32", Symbol("observedGeneration")=>"Int64", Symbol("updatedNumberScheduled")=>"Int32")
Base.propertynames(::Type{ IoK8sApiExtensionsV1beta1DaemonSetStatus }) = collect(keys(_property_map_IoK8sApiExtensionsV1beta1DaemonSetStatus))
Swagger.property_type(::Type{ IoK8sApiExtensionsV1beta1DaemonSetStatus }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IoK8sApiExtensionsV1beta1DaemonSetStatus[name]))}
Swagger.field_name(::Type{ IoK8sApiExtensionsV1beta1DaemonSetStatus }, property_name::Symbol) =  _property_map_IoK8sApiExtensionsV1beta1DaemonSetStatus[property_name]

function check_required(o::IoK8sApiExtensionsV1beta1DaemonSetStatus)
    (getproperty(o, Symbol("currentNumberScheduled")) === nothing) && (return false)
    (getproperty(o, Symbol("desiredNumberScheduled")) === nothing) && (return false)
    (getproperty(o, Symbol("numberMisscheduled")) === nothing) && (return false)
    (getproperty(o, Symbol("numberReady")) === nothing) && (return false)
    true
end

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


text:
module BINOD
    print("BINOD")
end


text:
using Compat.Random
import Distributions: ValueSupport, Sampleable
import Base: eltype, rand
using Phylo
using Distributions
using Missings
using IterableTables: getiterator

mutable struct Phylogenetics{T <: AbstractTree} <: ValueSupport end
Base.eltype(::Type{Phylogenetics{T}}) where T <: AbstractTree = T

"""
    Nonultrametric{T <: AbstractTree,
                   RNG <: Sampleable}(n::Int,
                                      rng::RNG = Exponential())
    Nonultrametric{T <: AbstractTree,
                   RNG <: Sampleable}(tiplabels::Vector{String},
                                      rng::RNG = Exponential())

The sampler for non-ultrametric phylogenetic trees of size `n` or with
tip labels `tiplabels`. Generate random trees by calling rand().
Currently only works for `NamedTree`s.
"""
struct Nonultrametric{T <: AbstractTree,
                      RNG <: Sampleable{Univariate, Continuous}} <:
    Sampleable{Univariate, Phylogenetics{T}}
    n::Int
    tiplabels::Vector{String}
    rng::RNG
    leafinfo::Any

    function Nonultrametric{T, RNG}(n::Int, tiplabels::Vector{String}, rng::RNG, leafinfo) where {T, RNG}
        return new{T, RNG}(n, tiplabels, rng, leafinfo)
    end

    function Nonultrametric{T, RNG}(n::Int, rng::RNG) where {T, RNG}
        return new{T, RNG}(n, map(i -> "tip $i", 1:n), rng, missing)
    end

    function Nonultrametric{T, RNG}(tiplabels::Vector{String}, rng::RNG) where {T, RNG}
        return new{T, RNG}(length(tiplabels), tiplabels, rng, missing)
    end
end

function Nonultrametric{T}(n::Int) where T <: AbstractTree
    return Nonultrametric{T, Exponential}(n, Exponential())
end

function Nonultrametric{T}(tiplabels::Vector{String}) where T <: AbstractTree
    return Nonultrametric{T, Exponential}(tiplabels, Exponential())
end

function Nonultrametric{T}(leafinfo) where T <: AbstractTree
    tipnames = unique(collect(map(info -> info[1], getiterator(leafinfo))))
    return Nonultrametric{T, Exponential}(length(tipnames), tipnames,
                                          Exponential(), leafinfo)
end

Nonultrametric(info::LI) where LI =
    Nonultrametric{BinaryTree{LI, Dict{String, Any}}}(info)

Nonultrametric(n::Int) = Nonultrametric{NamedTree}(n)
Nonultrametric(tiplabels::Vector{String}) = Nonultrametric{NamedTree}(tiplabels)

function rand(t::Nonultrametric{T, RNG}) where {T, RNG}
    t.n >= 2 || error("A tree must have at least 2 tips")
    if ismissing(t.leafinfo)
        tree = T(t.tiplabels; rootheight = 0.0)
    else
        tree = T(t.leafinfo; rootheight = 0.0)
    end
    roots = nodenamefilter(isroot, tree)
    while length(roots) > 1
        children = sample(collect(roots), 2, replace=false)
        parent = addnode!(tree)
        addbranch!(tree, parent, children[1], rand(t.rng))
        addbranch!(tree, parent, children[2], rand(t.rng))
        roots = nodenamefilter(isroot, tree)
    end
    return tree
end

"""
    Ultrametric{T <: AbstractTree,
                RNG <: Sampleable}(n::Int,
                                   rng::RNG = Exponential())
    Ultrametric{T <: AbstractTree,
                RNG <: Sampleable}(tiplabels::Vector{String},
                                   rng::RNG = Exponential())

The sampler for ultrametric phylogenetic trees of size `n` or with
tip labels `tiplabels`. Generate random trees by calling rand().
Currently only works for `NamedTree`s.
"""
struct Ultrametric{T <: AbstractTree,
                 RNG <: Sampleable{Univariate, Continuous}} <:
    Sampleable{Univariate, Phylogenetics{T}}
    n::Int
    tiplabels::Vector{String}
    rng::RNG
    leafinfo::Any

    function Ultrametric{T, RNG}(n::Int, tiplabels::Vector{String}, rng::RNG, leafinfo) where {T, RNG}
        return new{T, RNG}(n, tiplabels, rng, leafinfo)
    end

    function Ultrametric{T, RNG}(n::Int, rng::RNG) where {T, RNG}
        return new{T, RNG}(n, map(i -> "tip $i", 1:n), rng, missing)
    end

    function Ultrametric{T, RNG}(tiplabels::Vector{String}, rng::RNG) where {T, RNG}
        return new{T, RNG}(length(tiplabels), tiplabels, rng, missing)
    end
end

function Ultrametric{T}(n::Int) where T <: AbstractTree
    return Ultrametric{T, Exponential}(n, Exponential())
end

function Ultrametric{T}(tiplabels::Vector{String}) where T <: AbstractTree
    return Ultrametric{T, Exponential}(tiplabels, Exponential())
end

function Ultrametric{T}(leafinfo) where T <: AbstractTree
    tipnames = unique(collect(map(info -> info[1], getiterator(leafinfo))))
    return Ultrametric{T, Exponential}(length(tipnames), tipnames,
                                       Exponential(), leafinfo)
end

Ultrametric(info::LI) where LI =
    Ultrametric{BinaryTree{LI, Dict{String, Any}}}(info)

Ultrametric(n::Int) = Ultrametric{NamedTree}(n)
Ultrametric(tiplabels::Vector{String}) = Ultrametric{NamedTree}(tiplabels)

function rand(t::Ultrametric{T, RNG}) where {T, RNG}
    t.n >= 2 || error("A tree must have at least 2 tips")
    if ismissing(t.leafinfo)
        tree = T(t.tiplabels; rootheight = 0.0)
    else
        tree = T(t.leafinfo; rootheight = 0.0)
    end
    roots = nodenamefilter(isroot, tree)
    depth = zero(rand(t.rng))
    leaves = getleafnames(tree)
    while length(roots) > 1
        tocoalesce = collect(roots)
        coalescers = sample(tocoalesce, 2, replace=false)
        parent = addnode!(tree)
        depth += rand(t.rng) * 2.0 / length(tocoalesce)
        c1 = filter(x -> coalescers[1] ∈ nodehistory(tree, x), leaves)
        d1 = map(x -> getheight(tree, x), c1)
        c2 = filter(x -> coalescers[2] ∈ nodehistory(tree, x), leaves)
        d2 = map(x -> getheight(tree, x), c2)
        addbranch!(tree, parent, coalescers[1], depth - d1[1])
        addbranch!(tree, parent, coalescers[2], depth - d2[1])
        roots = nodenamefilter(isroot, tree)
    end
    return tree
end

function rand(s::S, treenames) where
    {TREE <: AbstractTree,
     S <: Sampleable{Univariate, Phylogenetics{TREE}}}
    trees = Dict{eltype(treenames), TREE}()
    for name in treenames
        trees[name] = rand(s)
    end
    return TreeSet(trees, Dict{String, Dict{String, Any}}())
end

function rand(s::S, n::Int) where
    {TREE <: AbstractTree,
     S <: Sampleable{Univariate, Phylogenetics{TREE}}}
    return rand(s, 1:n)
end


text:
@testset "120.triangle.jl" begin
    @test minimum_total([[2], [3, 4], [6, 5, 7], [4, 1, 8, 3]]) == 11
    @test minimum_total([[0]]) == 0
end


text:
#==============================================================================#
# AWSSDB.jl
#
# SDB API. See http://aws.amazon.com/documentation/simpledb/
#
# Copyright OC Technology Pty Ltd 2016 - All rights reserved
#==============================================================================#

__precompile__()

module AWSSDB

export sdb_list_domains, sdb_create_domain, sdb_delete_domain,
       sdb_delete_item, sdb_put, sdb_get, sdb_select,
       SimpleDB

using AWSCore
using XMLDict
using SymDict

# Associative get/setindex! interface...

type SimpleDB <:Associative{AbstractString,AbstractString}
    aws
    domain
end

function Base.get(db::SimpleDB, key::AbstractString, default)
    r = sdb_get(db, key)
    return r == nothing ? default : r
end

Base.setindex!(db::SimpleDB, v, key::AbstractString) = sdb_put(db, key, v)

# SDB API Interface...

function sdb(aws::AWSConfig, action, query::Dict)
    AWSCore.Services.sdb(aws, action, query)
end

sdb(aws::AWSConfig, action, query::SymbolDict) = sdb(aws, action, stringdict(query))

sdb(aws::AWSConfig, action; args...) = sdb(aws, action, stringdict(args))

function sdb_list_domains(aws::AWSConfig)
    r = sdb(aws, "ListDomains")["ListDomainsResult"]
    list = String[]
    if r != ""
        r = r["DomainName"]
        for d in (isa(r, Vector) ? r : [r])
            push!(list, d)
        end
    end
    return list
end

function sdb_create_domain(aws::AWSConfig, name)
    sdb(aws, "CreateDomain", DomainName = name)
end

function sdb_delete_domain(aws::AWSConfig, name)
    sdb(aws, "DeleteDomain", DomainName = name)
end

function sdb_delete_item(aws::AWSConfig, DomainName, ItemName)
    sdb(aws, "DeleteAttributes", @SymDict(DomainName, ItemName))
end

function sdb_put(aws::AWSConfig, DomainName, ItemName, dict; replace=true)

    attrs = StringDict()
    i = 1
    for (a, v) in dict
        if !isa(v, Vector)
            v = [v]
        end
        for v in v
            attrs["Attribute.$i.Name"] = a
            attrs["Attribute.$i.Value"] = v
            attrs["Attribute.$i.Replace"] = string(replace)
            i += 1
        end
    end

    sdb(aws, "PutAttributes", @SymDict(DomainName, ItemName, attrs...))
end

sdb_put(db::SimpleDB, ItemName, dict) = sdb_put(db.aws, db.domain, ItemName, dict)

const SDBAttribute = Pair{String,Union{String,Vector{String}}}
const SDBItem = Union{String,Pair{String,Vector{SDBAttribute}}}

# Convert XML attribute vector "v" to Vector{Pair}...

function sdb_attributes(v, order=nothing)

    v = isa(v, Vector) ? v : [v]

    # Build dict of attribute values...
    d = Dict{String,Union{String,Vector{String}}}()
    for i in v

        # Get Attribute name and value...
        attribute = String(i["Name"])
        value = i["Value"]

        # Decode base64 values...
        if isa(value, XMLDict.XMLDictElement)
            if get(value, :encoding, nothing) == "base64"
                value = base64decode(value[""])
            end
        end
        value = String(value)

        # Build vector for multi-valued attributes...
        if haskey(d, attribute)
            if isa(d[attribute], Vector{String})
                push!(d[attribute], value)
            else
                d[attribute] = String[d[attribute], value]
            end
        else
            d[attribute] = value
        end
    end

    if order == nothing
        return [SDBAttribute(k, v) for (k,v) in d]
    else
        # Return ordered Vector of Pairs...
        return [SDBAttribute(a, get(d, a, String(""))) for a in order]
    end
end

function sdb_get(aws::AWSConfig, DomainName, ItemName; attribute="")

    query = @SymDict(DomainName, ItemName)

    if attribute != ""
        query[:AttributeName] = attribute
    end

    r = sdb(aws, "GetAttributes", query)
    r = r["GetAttributesResult"]
    if r == ""
        return nothing
    end
    r = r["Attribute"]

    return sdb_attributes(r)
end

sdb_get(db::SimpleDB, ItemName) = sdb_get(db.aws, db.domain, ItemName)

type SDBSelectItr
    aws
    SelectExpression
    attributes
    NextToken::String
    result
end

function sdb_select(aws::AWSConfig, SelectExpression; NextToken="")
    return SDBSelectItr(aws, SelectExpression, nothing, NextToken, nothing)
end

function sdb_select(aws::AWSConfig, domain::AbstractString,
                                    attributes::Vector,
                                    condition::AbstractString="";
                                    NextToken="")

    SelectExpression = """
                       select $(join(["`$a`" for a in attributes], ", "))
                       from `$domain`
                       $condition
                       """

    return SDBSelectItr(aws, SelectExpression, attributes, NextToken, nothing)
end

# Iterator interface...

Base.eltype(::Type{SDBSelectItr}) = SDBItem
if VERSION > v"0.5.0-dev+4340"
Base.iteratorsize(::Type{SDBSelectItr}) = Base.SizeUnknown()
end

function Base.start(itr::SDBSelectItr)

    itr.result, itr.NextToken = _sdb_select(itr.aws,
                                            itr.SelectExpression,
                                            itr.attributes)
    return start(itr.result)
end

function Base.next(itr::SDBSelectItr, state)

    if done(itr.result, state)

        itr.result, itr.NextToken =
            _sdb_select(itr.aws,
                        itr.SelectExpression,
                        itr.attributes,
                        itr.NextToken)
        state = start(itr.result)
    end

    return next(itr.result, state)
end

function Base.done(itr::SDBSelectItr, state)

    return done(itr.result, state) && itr.NextToken == ""
end

function _sdb_select(aws::AWSConfig, SelectExpression, attributes, NextToken="")

    query = @SymDict(SelectExpression)

    if NextToken != ""
        query[:NextToken] = NextToken
    end

    r = sdb(aws, "Select", query)
    r = r["SelectResult"]

    if r == ""
        return (SDBItem[], "")
    end

    NextToken = get(r, "NextToken", "")

    r = r["Item"]

    if !isa(r, Vector)
        r = [r]
    end

    if ismatch(r"select itemName()", SelectExpression)
        r = SDBItem[String(i["Name"]) for i in r]
    else
        r = SDBItem[String(i["Name"]) =>
                    sdb_attributes(get(i, "Attribute", []), attributes)
                    for i in r]
    end

    return (r, NextToken)
end

end # module AWSSDB

#==============================================================================#
# End of file.
#==============================================================================#


text:
using LibPQ
using ..PostgresORMUtil, ..Controller

function SchemaInfo.analyse_db_schema(dbconn::LibPQ.Connection)

   result = Dict()

   for schema in SchemaInfo.get_schemas(dbconn)
      result[schema] = Dict()
      for tbl in SchemaInfo.get_tables(schema, dbconn)
         result[schema][tbl] = Dict(
            :pk => SchemaInfo.get_pks(tbl, schema, dbconn),
            :fks => SchemaInfo.get_fks(tbl, schema, dbconn),
            :cols => SchemaInfo.get_columns_types(tbl,schema, dbconn),
            :is_partition => SchemaInfo.check_if_table_is_partition_of_another_table(tbl, schema, dbconn)
         )
      end # ENDOF for SchemaInfo.get_tables
   end # ENDOF for schema

   return result

end

function SchemaInfo.get_schemas(dbconn::LibPQ.Connection)

   querystr = "SELECT n.nspname AS \"Name\",
               pg_catalog.pg_get_userbyid(n.nspowner) AS \"Owner\"
               FROM pg_catalog.pg_namespace n
               WHERE n.nspname !~ '^pg_' AND n.nspname <> 'information_schema'
               ORDER BY 1;"

   queryres = execute_plain_query(querystr, [], dbconn)
   return convert(Vector{String}, queryres.Name)
end

function SchemaInfo.get_tables(schema::String, dbconn::LibPQ.Connection)

   querystr = "SELECT n.nspname as \"Schema\",
                 c.relname as \"Name\",
                 CASE c.relkind WHEN 'r' THEN 'table' WHEN 'v' THEN 'view' WHEN 'm' THEN 'materialized view' WHEN 'i' THEN 'index' WHEN 'S' THEN 'sequence' WHEN 's' THEN 'special' WHEN 'f' THEN 'foreign table' WHEN 'p' THEN 'partitioned table' WHEN 'I' THEN 'partitioned index' END as \"Type\",
                 pg_catalog.pg_get_userbyid(c.relowner) as \"Owner\"
               FROM pg_catalog.pg_class c
                    LEFT JOIN pg_catalog.pg_namespace n ON n.oid = c.relnamespace
               WHERE c.relkind IN ('r','p','s','')
                     AND n.nspname !~ '^pg_toast'
                 AND n.nspname OPERATOR(pg_catalog.~) \$1 COLLATE pg_catalog.default
               ORDER BY 1,2;"
   println(querystr)
   schema_pattern = "^($schema)\$"
   queryres = execute_plain_query(querystr, [schema_pattern], dbconn)
   return convert(Vector{String}, queryres.Name)

end

function SchemaInfo.check_if_table_is_partition_of_another_table(
            table::String,
            schema::String,
            dbconn::LibPQ.Connection)

    #   SELECT c.relchecks, c.relkind, c.relhasindex, c.relhasrules, c.relhastriggers, c.relrowsecurity, c.relforcerowsecurity, false AS relhasoids, c.relispartition, '', c.reltablespace, CASE WHEN c.reloftype = 0 THEN '' ELSE c.reloftype::pg_catalog.regtype::pg_catalog.text END, c.relpersistence, c.relreplident, am.amname
    # FROM pg_catalog.pg_class c
    #  LEFT JOIN pg_catalog.pg_class tc ON (c.reltoastrelid = tc.oid)
    # LEFT JOIN pg_catalog.pg_am am ON (c.relam = am.oid)
    # WHERE c.oid = '116102';

   oid = SchemaInfo.get_table_oid(table, schema, dbconn)
   querystr = "SELECT c.relispartition
       FROM pg_catalog.pg_class c
        LEFT JOIN pg_catalog.pg_class tc ON (c.reltoastrelid = tc.oid)
       LEFT JOIN pg_catalog.pg_am am ON (c.relam = am.oid)
       WHERE c.oid = \$1;"

    queryres = execute_plain_query(querystr, [oid], dbconn)
    return queryres.relispartition[1]

end

function SchemaInfo.check_if_table_or_partition_exists(table::String,
                                                       schema::String,
                                                       dbconn::LibPQ.Connection)
    querystr = "SELECT EXISTS (
                   SELECT FROM information_schema.tables
                   WHERE  table_name  = \$1
                   AND    table_schema = \$2
                   );"
    queryres = execute_plain_query(querystr, [table, schema], dbconn)
    return queryres[1,1]
end

function SchemaInfo.get_table_oid(table::String, schema::String, dbconn::LibPQ.Connection)

   # pagila=# SELECT c.oid,
   # pagila-#   n.nspname,
   # pagila-#   c.relname
   # pagila-# FROM pg_catalog.pg_class c
   # pagila-#      LEFT JOIN pg_catalog.pg_namespace n ON n.oid = c.relnamespace
   # pagila-# WHERE c.relname OPERATOR(pg_catalog.~) '^(film_actor)$' COLLATE pg_catalog.default
   # pagila-#   AND pg_catalog.pg_table_is_visible(c.oid)
   # pagila-# ORDER BY 2, 3;
   #   oid   | nspname |  relname
   # --------+---------+------------
   #  106295 | public  | film_actor
   # (1 row)

   querystr = "SELECT c.oid,
               n.nspname,
               c.relname
               FROM pg_catalog.pg_class c
                     LEFT JOIN pg_catalog.pg_namespace n ON n.oid = c.relnamespace
                WHERE c.relname OPERATOR(pg_catalog.~) \$1 COLLATE pg_catalog.default
                  AND n.nspname OPERATOR(pg_catalog.~) \$2 COLLATE pg_catalog.default
                  ORDER BY 2, 3;"

   table_pattern  = "^($table)\$"
   schema_pattern = "^($schema)\$"

   queryres = execute_plain_query(querystr, [table_pattern,schema_pattern], dbconn)
   return signed(queryres.oid[1])

end

function SchemaInfo.get_pks(table::String, schema::String, dbconn::LibPQ.Connection)

   # pagila=# SELECT c2.relname, i.indisprimary, i.indisunique, i.indisclustered, i.indisvalid, pg_catalog.pg_get_indexdef(i.indexrelid, 0, true),
   # pagila-#   pg_catalog.pg_get_constraintdef(con.oid, true), contype, condeferrable, condeferred, i.indisreplident, c2.reltablespace
   # pagila-# FROM pg_catalog.pg_class c, pg_catalog.pg_class c2, pg_catalog.pg_index i
   # pagila-#   LEFT JOIN pg_catalog.pg_constraint con ON (conrelid = i.indrelid AND conindid = i.indexrelid AND contype IN ('p','u','x'))
   # pagila-# WHERE c.oid = '106295' AND c.oid = i.indrelid AND i.indexrelid = c2.oid
   # pagila-# ORDER BY i.indisprimary DESC, i.indisunique DESC, c2.relname;
   #      relname     | indisprimary | indisunique | indisclustered | indisvalid |                                  pg_get_indexdef                                  |      pg_get_constraintdef       | contype | condeferrable | condeferred | i
   # ndisreplident | reltablespace
   # -----------------+--------------+-------------+----------------+------------+-----------------------------------------------------------------------------------+---------------------------------+---------+---------------+-------------+--
   # --------------+---------------
   #  film_actor_pkey | t            | t           | f              | t          | CREATE UNIQUE INDEX film_actor_pkey ON film_actor USING btree (actor_id, film_id) | PRIMARY KEY (actor_id, film_id) | p       | f             | f           | f
   #               |             0
   #  idx_fk_film_id  | f            | f           | f              | t          | CREATE INDEX idx_fk_film_id ON film_actor USING btree (film_id)                   |                                 |         |               |             | f
   #               |             0
   # (2 rows)

   querystr = "SELECT c2.relname, i.indisprimary, i.indisunique, i.indisclustered, i.indisvalid, pg_catalog.pg_get_indexdef(i.indexrelid, 0, true),
               pg_catalog.pg_get_constraintdef(con.oid, true), contype, condeferrable, condeferred, i.indisreplident, c2.reltablespace
                FROM pg_catalog.pg_class c, pg_catalog.pg_class c2, pg_catalog.pg_index i
                  LEFT JOIN pg_catalog.pg_constraint con ON (conrelid = i.indrelid AND conindid = i.indexrelid AND contype IN ('p','u','x'))
                WHERE c.oid = \$1 AND c.oid = i.indrelid AND i.indexrelid = c2.oid
                ORDER BY i.indisprimary DESC, i.indisunique DESC, c2.relname;"

   table_oid = SchemaInfo.get_table_oid(table, schema, dbconn)

   queryres = execute_plain_query(querystr, [table_oid], dbconn)
   pg_get_constraintdef = filter(x -> begin
                                          !ismissing(x) && startswith(x,"PRIMARY KEY")
                                      end,
                                 queryres.pg_get_constraintdef)

   pg_get_constraintdef =  convert(Vector{String}, pg_get_constraintdef)
   # pg_get_constraintdef = map(x -> replace(x, " " => ""), pg_get_constraintdef)
   # return pg_get_constraintdef

   result = Vector{String}()
   for str in pg_get_constraintdef
      pks = match(r"^PRIMARY KEY \((.*)\)$",str)
      pks =  convert(Vector{String}, pks.captures)
      tmp_res = String[]
      for pk in pks
         push!(tmp_res,
               remove_spaces_and_split(pk)...)
      end
      sort!(tmp_res) # IMPORTANT! Sort the PK columns, we'll do the same with
                     #            the FK columns
      push!(result, tmp_res...)
   end

   return result

end

function SchemaInfo.get_fks(table::String, schema::String, dbconn::LibPQ.Connection)

   # pagila=# SELECT true as sametable, conname,
   # pagila-#   pg_catalog.pg_get_constraintdef(r.oid, true) as condef,
   # pagila-#   conrelid::pg_catalog.regclass AS ontable
   # pagila-# FROM pg_catalog.pg_constraint r
   # pagila-# WHERE r.conrelid = '106295' AND r.contype = 'f'
   # pagila-#      AND conparentid = 0
   # pagila-# ORDER BY conname
   # pagila-#
   # pagila-# ;
   #  sametable |         conname          |                                         condef                                         |  ontable
   # -----------+--------------------------+----------------------------------------------------------------------------------------+------------
   #  t         | film_actor_actor_id_fkey | FOREIGN KEY (actor_id) REFERENCES actor(actor_id) ON UPDATE CASCADE ON DELETE RESTRICT | film_actor
   #  t         | film_actor_film_release_year_fkey  | FOREIGN KEY (film_id, film_release_year) REFERENCES film(film_id, release_year)    | film_actor
   # (2 rows)

   querystr = "SELECT true as sametable, conname,
               pg_catalog.pg_get_constraintdef(r.oid, true) as condef,
                  conrelid::pg_catalog.regclass AS ontable
                FROM pg_catalog.pg_constraint r
                WHERE r.conrelid = \$1 AND r.contype = 'f'
                     AND conparentid = 0
                ORDER BY conname"

   table_oid = SchemaInfo.get_table_oid(table, schema, dbconn)

   queryres = execute_plain_query(querystr, [table_oid], dbconn)
   queryres = filter(x -> begin
                           !ismissing(x.condef) && startswith(x.condef,"FOREIGN KEY")
                       end,
                    queryres)

   result = Dict()
   for r in eachrow(queryres)
      fks = match(r"^FOREIGN KEY \((.*)\) REFERENCES (.*)\((.*)\)",r.condef)
      # @info fks.captures

      # Referencing columns
      referencing_cols = remove_spaces_and_split(string(fks.captures[1]))
      referencing_cols = string.(referencing_cols)

      # Referenced table
      referenced_table = remove_spaces_and_split(string(fks.captures[2]))[1]
      referenced_schema = "public"
      if occursin(".",referenced_table)
         referenced_table_arr = split(referenced_table,'.')
         referenced_table = referenced_table_arr[2]
         referenced_schema = referenced_table_arr[1]
      end

      # Referenced columns
      referenced_cols = remove_spaces_and_split(string(fks.captures[3]))
      referenced_cols = string.(referenced_cols)

      # Reorder the pairs (PK column, FK column) according to the PK columns
      ordered_cols_along_pk_cols =
         collect(zip(referenced_cols, referencing_cols)) |>
         n -> sort(n, by = x -> x[1])
      referenced_cols = map(x -> x[1],ordered_cols_along_pk_cols)
      referencing_cols = map(x -> x[2],ordered_cols_along_pk_cols)

      result[r.conname] = Dict(
         :referencing_cols => referencing_cols,
         :referenced_table => Dict(:table => string(referenced_table),
                                   :schema => string(referenced_schema)),
         :referenced_cols => referenced_cols
      )

   end

   return result

end

function SchemaInfo.get_columns_types(table::String,schema::String, dbconn::LibPQ.Connection)
   query_string = "SELECT column_name,
                          data_type AS column_type,
                          udt_name AS element_type
                   from information_schema.columns
                   WHERE table_schema = \$1 AND table_name = \$2"

   cols = execute_plain_query(query_string,
                              [schema,table],
                              dbconn)

    result = Dict()

   for c in eachrow(cols)
      colname = c.column_name
      coltype = c.column_type
      elttype = c.element_type
      result[colname] = Dict(
         :type => coltype,
         :elttype_if_array => elttype
      )
   end

   return result

end

function SchemaInfo.get_custom_types(dbconn::LibPQ.Connection)

  #  schema |    name    | internal_name | size |  elements  | description
  # --------+------------+---------------+------+------------+-------------
  #  public | value_type | value_type    | 4    | ordinal   +|
  #         |            |               |      | continuous+|
  #         |            |               |      | category  +|
  #         |            |               |      | bool      +|
  #         |            |               |      | text       |

   query_string = "SELECT n.nspname AS schema,
        pg_catalog.format_type ( t.oid, NULL ) AS name,
        t.typname AS internal_name,
        CASE
            WHEN t.typrelid != 0
            THEN CAST ( 'tuple' AS pg_catalog.text )
            WHEN t.typlen < 0
            THEN CAST ( 'var' AS pg_catalog.text )
            ELSE CAST ( t.typlen AS pg_catalog.text )
        END AS size,
        pg_catalog.array_to_string (
            ARRAY( SELECT e.enumlabel
                    FROM pg_catalog.pg_enum e
                    WHERE e.enumtypid = t.oid
                    ORDER BY e.oid ), E'\n'
            ) AS elements,
        pg_catalog.obj_description ( t.oid, 'pg_type' ) AS description
    FROM pg_catalog.pg_type t
    LEFT JOIN pg_catalog.pg_namespace n
        ON n.oid = t.typnamespace
    WHERE ( t.typrelid = 0
            OR ( SELECT c.relkind = 'c'
                    FROM pg_catalog.pg_class c
                    WHERE c.oid = t.typrelid
                )
        )
        AND NOT EXISTS
            ( SELECT 1
                FROM pg_catalog.pg_type el
                WHERE el.oid = t.typelem
                    AND el.typarray = t.oid
            )
        AND n.nspname <> 'pg_catalog'
        AND n.nspname <> 'information_schema'
        AND pg_catalog.pg_type_is_visible ( t.oid )
    ORDER BY 1, 2;"

   cols = execute_plain_query(query_string,
                              [],
                              dbconn)

   result = Dict()

   for c in eachrow(cols)
      name_ = c.internal_name
      possible_values = string.(split(c.elements, "\n"))
      # Cleaning
      filter!(x-> length(x) > 0,possible_values)
      result[name_] = Dict(
         :name => name_,
         :possible_values => possible_values
      )
   end

   return result

end


text:
"""
    LegendreDerivativeOperator{T<:Real}

A derivative operator on a nonperiodic Lobatto-Legendre grid with scalar type
`T` computing the first derivative using a Legendre expansion.
"""
struct LegendreDerivativeOperator{T<:Real} <: AbstractNonperiodicDerivativeOperator{T}
    jac::T
    Δx::T
    grid::Vector{T} # N nodes, including the left and the right boundary
    basis::LobattoLegendre{T}

    function LegendreDerivativeOperator(xmin::T, xmax::T, basis::LobattoLegendre{T}) where {T<:Real}
        grid = map_from_canonical.(basis.nodes, xmin, xmax, basis)
        jac = 2 / (xmax - xmin)
        Δx = inv(jac)

        new{T}(jac, Δx, grid, basis)
    end
end

"""
    LegendreDerivativeOperator(xmin::T, xmax::T, N::Int) where {T<:Real}

Construct the `LegendreDerivativeOperator` on a uniform grid between `xmin` and
`xmax` using `N` nodes and `N-1` Legendre modes.
"""
function LegendreDerivativeOperator(xmin::T, xmax::T, N::Int) where {T<:Real}
    @argcheck N >= 2

    basis = LobattoLegendre(N-1, T)

    LegendreDerivativeOperator(xmin, xmax, basis)
end

"""
    legendre_derivative_operator(xmin::Real, xmax::Real, N::Integer)
    legendre_derivative_operator(; xmin::Real, xmax::Real, N::Integer)

Construct the `LegendreDerivativeOperator` on a uniform grid between `xmin` and
`xmax` using `N` nodes and `N-1` Legendre modes.
"""
function legendre_derivative_operator(xmin::Real, xmax::Real, N::Integer)
    LegendreDerivativeOperator(promote(xmin, xmax)..., N)
end

function legendre_derivative_operator(; xmin::Real, xmax::Real, N::Integer)
    legendre_derivative_operator(xmin, xmax, N)
end

derivative_order(D::LegendreDerivativeOperator) = 1
LinearAlgebra.issymmetric(D::LegendreDerivativeOperator) = false

"""
    LegendreSecondDerivativeOperator{T<:Real}

A derivative operator on a nonperiodic Lobatto-Legendre grid with scalar type
`T` computing the second derivative using a Legendre expansion.
"""
struct LegendreSecondDerivativeOperator{T<:Real} <: AbstractNonperiodicDerivativeOperator{T}
    D2::Matrix{T}
    jac::T
    Δx::T
    grid::Vector{T} # N nodes, including the left and the right boundary
    basis::LobattoLegendre{T}

    function LegendreSecondDerivativeOperator(xmin::T, xmax::T, N::Int) where {T<:Real}
        @argcheck N >= 2

        basis = LobattoLegendre(N-1, T)

        D2 = basis.D^2
        grid = map_from_canonical.(basis.nodes, xmin, xmax, basis)
        jac = 2 / (xmax - xmin)
        Δx = inv(jac)

        new{T}(D2, jac, Δx, grid, basis)
    end
end

"""
    legendre_second_derivative_operator(xmin::Real, xmax::Real, N::Integer)
    legendre_second_derivative_operator(; xmin::Real, xmax::Real, N::Integer)

Construct the `LegendreDerivativeOperator` on a uniform grid between `xmin` and
`xmax` using `N` nodes and `N-1` Legendre modes.
"""
function legendre_second_derivative_operator(xmin::Real, xmax::Real, N::Integer)
  LegendreSecondDerivativeOperator(promote(xmin, xmax)..., N)
end

function legendre_second_derivative_operator(; xmin::Real, xmax::Real, N::Integer)
    legendre_second_derivative_operator(xmin, xmax, N)
end

derivative_order(D::LegendreSecondDerivativeOperator) = 2

integrate(func, u, D::Union{LegendreDerivativeOperator,LegendreSecondDerivativeOperator}) = D.Δx*integrate(func, u, D.basis)
function evaluate_coefficients(u, D::Union{LegendreDerivativeOperator,LegendreSecondDerivativeOperator},
                               npoints=2*size(D,2)+1)
    evaluate_coefficients(u, D.basis, npoints)
end

mass_matrix(D::Union{LegendreDerivativeOperator,LegendreSecondDerivativeOperator}) = Diagonal(D.Δx * D.basis.weights)

Base.eltype(D::Union{LegendreDerivativeOperator{T},LegendreSecondDerivativeOperator{T}}) where {T} = T

function scale_by_mass_matrix!(u::AbstractVector, D::Union{LegendreDerivativeOperator,LegendreSecondDerivativeOperator}, factor=true)
    @unpack Δx, basis = D
    N, _ = size(D)
    @boundscheck begin
        @argcheck N == length(u)
    end

    @inbounds @simd for i in eachindex(u, basis.weights)
        u[i] = factor * u[i] * (Δx * basis.weights[i])
    end

    u
end

function scale_by_inverse_mass_matrix!(u::AbstractVector, D::Union{LegendreDerivativeOperator,LegendreSecondDerivativeOperator}, factor=true)
    @unpack Δx, basis = D
    N, _ = size(D)
    @boundscheck begin
        @argcheck N == length(u)
    end

    @inbounds @simd for i in eachindex(u, basis.weights)
        u[i] = factor * u[i] / (Δx * basis.weights[i])
    end

    u
end

function get_weight(D::Union{LegendreDerivativeOperator,LegendreSecondDerivativeOperator}, i::Int)
    @unpack Δx, basis = D
    @unpack weights = basis
    N, _ = size(D)
    @boundscheck begin
        @argcheck 1 <= i <= N
    end
    @inbounds ω = Δx * weights[i]
    ω
end

function Base.show(io::IO, D::LegendreDerivativeOperator)
    if get(io, :compact, false)
        summary(io, D)
    else
        x = grid(D)
        print(io, "First derivative operator {T=", eltype(D), "}")
        print(io, " on ", length(x), " Lobatto Legendre nodes in [", first(x), ", ", last(x), "]")
    end
end

function Base.show(io::IO, D::LegendreSecondDerivativeOperator)
    if get(io, :compact, false)
        summary(io, D)
    else
        x = grid(D)
        print(io, "Second derivative operator {T=", eltype(D), "}")
        print(io, " on ", length(x), " Lobatto Legendre nodes in [", first(x), ", ", last(x), "]")
    end
end

function mul!(dest::AbstractVector, D::LegendreDerivativeOperator, u::AbstractVector, α=true, β=false)
    @unpack jac, basis = D
    N, _ = size(D)
    @boundscheck begin
        @argcheck N == length(u)
        @argcheck N == length(dest)
    end

    mul!(dest, basis.D, u, α*jac, β)
end

function mul!(dest::AbstractVector, D::LegendreSecondDerivativeOperator, u::AbstractVector, α=true, β=false)
    @unpack jac, D2 = D
    N, _ = size(D2)
    @boundscheck begin
        @argcheck N == length(u)
        @argcheck N == length(dest)
    end

    mul!(dest, D2, u, α*jac^2, β)
end

function derivative_left(D::LegendreSecondDerivativeOperator, u::AbstractVector, ::Val{1})
    @unpack jac, basis = D
    N, _ = size(D)
    @boundscheck begin
        @argcheck N == length(u)
    end

    jac * dot(view(basis.D, 1, :), u)
end

function derivative_right(D::LegendreSecondDerivativeOperator, u::AbstractVector, ::Val{1})
    @unpack jac, basis = D
    N, _ = size(D)
    @boundscheck begin
        @argcheck N == length(u)
    end

    jac * dot(view(basis.D, size(D,1), :), u)
end

function lower_bandwidth(D::Union{LegendreDerivativeOperator,LegendreSecondDerivativeOperator})
    size(D, 1) - 1
end

function upper_bandwidth(D::Union{LegendreDerivativeOperator,LegendreSecondDerivativeOperator})
    size(D, 1) - 1
end

function accuracy_order(D::Union{LegendreDerivativeOperator,LegendreSecondDerivativeOperator})
    size(D, 1) - 1
end

function left_boundary_weight(D::Union{LegendreDerivativeOperator,LegendreSecondDerivativeOperator})
    @inbounds retval = D.Δx * D.basis.weights[1]
    retval
end

function right_boundary_weight(D::Union{LegendreDerivativeOperator,LegendreSecondDerivativeOperator})
    @inbounds retval = D.Δx * D.basis.weights[end]
    retval
end

function ConstantFilter(basis::LobattoLegendre{T}, filter, tmp::Array) where {T}
    Np1 = length(basis.nodes)
    coefficients = Array{T}(undef, Np1)
    set_filter_coefficients!(coefficients, filter)
    modal2nodal = legendre_vandermonde(basis)
    nodal2modal = FactorisationWrapper(factorize(modal2nodal))

    ConstantFilter(coefficients, nodal2modal, modal2nodal, tmp, filter)
end

function ConstantFilter(basis::LobattoLegendre{T}, filter, TmpEltype::DataType=T) where {T}
    Np1 = length(basis.nodes)
    tmp = Array{TmpEltype}(undef, Np1)

    ConstantFilter(basis, filter, tmp)
end

"""
    ConstantFilter(D::LegendreDerivativeOperator, filter, TmpEltype=T)

Create a modal filter with constant parameters adapted to the Legendre
derivative operator `D` with parameters given by the filter function `filter`.
"""
function ConstantFilter(D::LegendreDerivativeOperator{T}, filter, TmpEltype=T) where {T}
    ConstantFilter(D.basis, filter, TmpEltype)
end

# TODO: LegendreSuperSpectralViscosity


text:
module TestSimulate
  using Test
  using BioEnergeticFoodWebs

  food_chain = [0 0 0;0 0 1; 1 0 0]
  p = model_parameters(food_chain)
  # Fail if the biomass vector has bad length
  @test_throws AssertionError simulate(p, rand(5))
  # Fail if stop < start
  @test_throws AssertionError simulate(p, rand(3), start=0, stop=-1)
end

module TestSimulateHandChecked
  using Test
  using BioEnergeticFoodWebs

  A = [0 1 0; 0 0 0; 0 1 0]
  p = model_parameters(A)
  b0 = vec([0.2 0.4 0.1])
  der = similar(b0)
  BioEnergeticFoodWebs.dBdt(der, b0, p, 0.0)
  # 0.1604888888888889,-0.504296,0.08024444444444445
  @test der[1] ≈ 0.160 atol = 0.01
  @test der[2] ≈ -0.5 atol = 0.01
  @test der[3] ≈ 0.080 atol = 0.01
end

module TestSimulateSanityCheck
  using Test
  using BioEnergeticFoodWebs
  using Statistics

  # A producer with no predation reaches K
  free_producers = [0 1 0 0; 0 0 1 0; 0 0 0 0; 0 0 0 0]
  parameters = free_producers |> model_parameters
  bm = rand(4)

  s = simulate(parameters, bm, start=0, stop=25)
  @test s[:B][end,4] ≈ parameters[:K] atol=0.002

  # Using system-wide regulation, producers with no consumption reach K / n
  A = zeros(Int64, (4, 4))
  parameters = model_parameters(A, productivity=:system)
  bm = rand(4)
  s = simulate(parameters, bm, start=0, stop=15)
  @test s[:B][end,4] ≈ parameters[:K]/4 atol=0.001

  # A consumer with a resource with 0 biomass goes extinct
  A = [0 1; 0 0]
  parameters = A |> model_parameters
  bm = vec([1.0, 0.0])

  # s = simulate(p, n, start=0, stop=50, use=:stiff)
  # @test_approx_eq_eps s[:B][end,1] 0.0 0.001
  # @test_approx_eq_eps s[:B][end,2] 0.0 0.001

  s = simulate(parameters, bm, start=0, stop=50)
  @test s[:B][end,1] ≈ .0 atol=0.001
  @test s[:B][end,2] ≈ 0.0 atol=0.001

  #test that stiff and nonstiff methods yield similar results
  A = [0 0 1 ; 0 0 1 ; 0 0 0]
  parameters = model_parameters(A)
  b = [.5, .5, .5]
  s_stiff = simulate(parameters, b, use = :stiff)
  meanB_stiff = mean(s_stiff[:B], dims = 1)
  s_nonstiff = simulate(parameters, b, use = :nonstiff)
  meanB_nonstiff = mean(s_stiff[:B], dims = 1)
  @test meanB_stiff == meanB_nonstiff
end

max_stars_repo_name: SimmonsBI/BioEnergeticFoodWebs

max_issues_repo_name: SimmonsBI/BioEnergeticFoodWebs

max_forks_repo_name: SimmonsBI/BioEnergeticFoodWebs


text:
__precompile__()

module GetGene

using HTTP, LazyJSON, DataFrames

export
    getgenes,
    getgeneinfo

include("getgenes.jl")

end


text:
export freqCont;

"""
	function freqCont

	Frequency continuation procedure for running FWI.
	This function runs GaussNewton on misfit functions defined by the continuation division array contDiv.

	Input:
		mc    		- current model
		pInv		- Inverse param
		pMis 		- misfit params (remote)
		contDiv		- continuation division. Assumes that pMis is contiunous with respecto to the division.
					  If the tasks (frequencies) are divided in pMis {1,2} {3,4} {5,6} then contDiv = [1,3,5,7] (similarly to the ptr array in SparseMatrixCSC)
		windowSize  - How many frequencies to treat at once at the most.
		resultsFilename - a filename for saving the intermediate results according to the GN and continuation (FC) iterations (done in dumpFun())
		dumpFun     - a function for plotting, saving and doing all the things with the intermidiate results.
		mode        - either "1stInit" or anything else.
					  1stInit will use the first group of misfits as an initialization and will not include it together with the next group of misfits.
		startFrom   - a start index for the continuation. Usefull when our run broke in the middle of the iterations.
		cycle       - just for identifying different runs when saving the results.
		method      - "projGN" or "barrierGN"

"""
function freqCont(mc, pInv::InverseParam, pMis::Array{RemoteChannel},contDiv::Array{Int64}, windowSize::Int64,
			resultsFilename::String,dumpFun::Function,mode::String="",startFrom::Int64 = 1,cycle::Int64=0,method::String="projGN")
Dc = 0;
flag = -1;
HIS = [];
for freqIdx = startFrom:(length(contDiv)-1)
	if mode=="1stInit"
		reqIdx1 = freqIdx;
		if freqIdx > 1
			reqIdx1 = max(2,freqIdx-windowSize+1);
		end
		reqIdx2 = freqIdx;
	else
		reqIdx1 = freqIdx;
		if freqIdx > 1
			reqIdx1 = max(1,freqIdx-windowSize+1);
		end
		reqIdx2 = freqIdx;
	end
	currentProblems = contDiv[reqIdx1]:contDiv[reqIdx2+1]-1;
	println("\n======= New Continuation Stage: selecting continuation batches: ",reqIdx1," to ",reqIdx2,"=======\n");
	pMisTemp = pMis[currentProblems];
	pInv.mref = mc[:];

	if resultsFilename == ""
		filename = "";
		hisMatFileName = "";
	else
		Temp = splitext(resultsFilename);
		if cycle==0
			filename = string(Temp[1],"_FC",freqIdx,"_GN",Temp[2]);
			hisMatFileName  = string(Temp[1],"_FC",freqIdx);
		else
			filename = string(Temp[1],"_Cyc",cycle,"_FC",freqIdx,"_GN",Temp[2]);
			hisMatFileName  =  string(Temp[1],"_Cyc",cycle,"_FC",freqIdx);
		end
	end

	# Here we set a dump function for GN for this iteracion of FC
	function dumpGN(mc,Dc,iter,pInv,PF)
		dumpFun(mc,Dc,iter,pInv,PF,filename);
	end

	if method == "projGN"
		mc,Dc,flag,His = projGNCG(mc,pInv,pMisTemp,dumpResults = dumpGN);
	elseif method == "barrierGN"
		mc,Dc,flag,His = barrierGNCG(mc,pInv,pMisTemp,rho=1.0,dumpResults = dumpGN);
	end

	if hisMatFileName != ""
		file = matopen(string(hisMatFileName,"_HisGN.mat"), "w");
		### PUT THE CONTENT OF "HIS" INTO THE MAT FILE ###
		His.Dc = [];
		write(file,"His",His);
		close(file);
	end
	His.Dc = []
	push!(HIS,His)

	clear!(pMisTemp);
end
return mc,Dc,flag,HIS;
end