text:
default_user() = LibGit2.getconfig("github.user", "")
default_version() = VersionNumber(VERSION.major)
default_plugins() = [
    CompatHelper(),
    ProjectFile(),
    SrcDir(),
    Git(),
    License(),
    Readme(),
    Tests(),
    TagBot(),
]

function default_authors()
    name = LibGit2.getconfig("user.name", "")
    isempty(name) && return "contributors"
    email = LibGit2.getconfig("user.email", "")
    authors = isempty(email) ? name : "$name <$email>"
    return "$authors and contributors"
end

"""
    Template(; kwargs...)

A configuration used to generate packages.

## Keyword Arguments

### User Options
- `user::AbstractString="$(default_user())"`: GitHub (or other code hosting service) username.
  The default value comes from the global Git config (`github.user`).
  If no value is obtained, many plugins that use this value will not work.
- `authors::Union{AbstractString, Vector{<:AbstractString}}="$(default_authors())"`: Package authors.
  Like `user`, it takes its default value from the global Git config
  (`user.name` and `user.email`).

### Package Options
- `dir::AbstractString="$(contractuser(Pkg.devdir()))"`: Directory to place packages in.
- `host::AbstractString="github.com"`: URL to the code hosting service where packages will reside.
- `julia::VersionNumber=$(repr(default_version()))`: Minimum allowed Julia version.

### Template Plugins
- `plugins::Vector{<:Plugin}=Plugin[]`: A list of [`Plugin`](@ref)s used by the template.
  The default plugins are [`ProjectFile`](@ref), [`SrcDir`](@ref), [`Tests`](@ref),
  [`Readme`](@ref), [`License`](@ref), [`Git`](@ref), [`CompatHelper`](@ref), and
  [`TagBot`](@ref).
  To disable a default plugin, pass in the negated type: `!PluginType`.
  To override a default plugin instead of disabling it, pass in your own instance.

### Interactive Mode
- `interactive::Bool=false`: In addition to specifying the template options with keywords,
  you can also build up a template by following a set of prompts.
  To create a template interactively, set this keyword to `true`.
  See also the similar [`generate`](@ref) function.

---

To create a package from a `Template`, use the following syntax:

```julia
julia> t = Template();

julia> t("PkgName")
```
"""
struct Template
    authors::Vector{String}
    dir::String
    host::String
    julia::VersionNumber
    plugins::Vector{<:Plugin}
    user::String
end

Template(; interactive::Bool=false, kwargs...) = Template(Val(interactive); kwargs...)
Template(::Val{true}; kwargs...) = interactive(Template; kwargs...)

function Template(::Val{false}; kwargs...)
    kwargs = Dict(kwargs)

    user = getkw!(kwargs, :user)
    dir = abspath(expanduser(getkw!(kwargs, :dir)))
    host = replace(getkw!(kwargs, :host), r".*://" => "")
    julia = getkw!(kwargs, :julia)

    authors = getkw!(kwargs, :authors)
    authors isa Vector || (authors = map(strip, split(authors, ",")))

    # User-supplied plugins come first, so that deduping the list will remove the defaults.
    plugins = Vector{Any}(collect(getkw!(kwargs, :plugins)))
    disabled = map(d -> first(typeof(d).parameters), filter(p -> p isa Disabled, plugins))
    filter!(p -> p isa Plugin, plugins)
    append!(plugins, filter(p -> !(typeof(p) in disabled), default_plugins()))
    plugins = Vector{Plugin}(sort(unique(typeof, plugins); by=string))

    if isempty(user)
        foreach(plugins) do p
            if needs_username(p)
                T = nameof(typeof(p))
                s = """$T: Git hosting service username is required, set one with keyword `user="<username>"`"""
                throw(ArgumentError(s))
            end
        end
    end

    if !isempty(kwargs)
        @warn "Unrecognized keywords were supplied, see the documentation for help" kwargs
    end

    t = Template(authors, dir, host, julia, plugins, user)
    foreach(p -> validate(p, t), t.plugins)
    return t
end

"""
    (::Template)(pkg::AbstractString)

Generate a package named `pkg` from a [`Template`](@ref).
"""
function (t::Template)(pkg::AbstractString)
    endswith(pkg, ".jl") && (pkg = pkg[1:end-3])
    pkg_dir = joinpath(t.dir, pkg)
    ispath(pkg_dir) && throw(ArgumentError("$pkg_dir already exists"))
    mkpath(pkg_dir)

    try
        foreach((prehook, hook, posthook)) do h
            @info "Running $(nameof(h))s"
            foreach(sort(t.plugins; by=p -> priority(p, h), rev=true)) do p
                h(p, t, pkg_dir)
            end
        end
    catch
        rm(pkg_dir; recursive=true, force=true)
        rethrow()
    end

    @info "New package is at $pkg_dir"
end

function Base.:(==)(a::Template, b::Template)
    return a.authors == b.authors &&
        a.dir == b.dir &&
        a.host == b.host &&
        a.julia == b.julia &&
        a.user == b.user &&
        all(map(==, a.plugins, b.plugins))
end

# Does the template have a plugin that satisfies some predicate?
hasplugin(t::Template, f::Function) = any(f, t.plugins)
hasplugin(t::Template, ::Type{T}) where T <: Plugin = hasplugin(t, p -> p isa T)

"""
    getplugin(t::Template, ::Type{T<:Plugin}) -> Union{T, Nothing}

Get the plugin of type `T` from the template `t`, if it's present.
"""
function getplugin(t::Template, ::Type{T}) where T <: Plugin
    i = findfirst(p -> p isa T, t.plugins)
    return i === nothing ? nothing : t.plugins[i]
end

# Get a keyword or a default value.
getkw!(kwargs, k) = pop!(kwargs, k, defaultkw(Template, k))

# Default Template keyword values.
defaultkw(::Type{T}, s::Symbol) where T = defaultkw(T, Val(s))
defaultkw(::Type{Template}, ::Val{:authors}) = default_authors()
defaultkw(::Type{Template}, ::Val{:dir}) = contractuser(Pkg.devdir())
defaultkw(::Type{Template}, ::Val{:host}) = "github.com"
defaultkw(::Type{Template}, ::Val{:julia}) = default_version()
defaultkw(::Type{Template}, ::Val{:plugins}) = Plugin[]
defaultkw(::Type{Template}, ::Val{:user}) = default_user()

function interactive(::Type{Template}; kwargs...)
    # If the user supplied any keywords themselves, don't prompt for them.
    kwargs = Dict{Symbol, Any}(kwargs)
    options = [:user, :authors, :dir, :host, :julia, :plugins]
    customizable = setdiff(options, keys(kwargs))

    # Make sure we don't try to show a menu with < 2 options.
    isempty(customizable) && return Template(; kwargs...)
    just_one = length(customizable) == 1
    just_one && push(customizable, "None")

    return try
        println("Template keywords to customize:")
        menu = MultiSelectMenu(map(string, customizable); pagesize=length(customizable))
        customize = customizable[sort!(collect(request(menu)))]
        just_one && lastindex(customizable) in customize && return Template(; kwargs...)

        # Prompt for each keyword.
        foreach(customize) do k
            kwargs[k] = prompt(Template, fieldtype(Template, k), k)
        end

        Template(; kwargs...)
    catch e
        e isa InterruptException || rethrow()
        println()
        @info "Cancelled"
        nothing
    end
end

prompt(::Type{Template}, ::Type, ::Val{:pkg}) = Base.prompt("Package name")

function prompt(::Type{Template}, ::Type, ::Val{:host})
    hosts = ["github.com", "gitlab.com", "bitbucket.org", "Other"]
    menu = RadioMenu(hosts; pagesize=length(hosts))
    println("Select Git repository hosting service:")
    idx = request(menu)
    return if idx == lastindex(hosts)
        fallback_prompt(String, :host)
    else
        hosts[idx]
    end
end

function prompt(::Type{Template}, ::Type, ::Val{:julia})
    versions = map(format_version, VersionNumber.(1, 0:VERSION.minor))
    push!(versions, "Other")
    menu = RadioMenu(map(string, versions); pagesize=length(versions))
    println("Select minimum Julia version:")
    idx = request(menu)
    return if idx == lastindex(versions)
        fallback_prompt(VersionNumber, :julia)
    else
        VersionNumber(versions[idx])
    end
end

const CR = "\r"
const DOWN = "\eOB"

function prompt(::Type{Template}, ::Type, ::Val{:plugins})
    defaults = map(typeof, default_plugins())
    ndefaults = length(defaults)
    # Put the defaults first.
    options = unique!([defaults; concretes(Plugin)])
    menu = MultiSelectMenu(map(T -> string(nameof(T)), options); pagesize=length(options))
    println("Select plugins:")
    # Pre-select the default plugins and move the cursor to the first non-default.
    # To make this better, we need julia#30043.
    print(stdin.buffer, (CR * DOWN)^ndefaults)
    types = sort!(collect(request(menu)))
    plugins = Vector{Any}(map(interactive, options[types]))
    # Find any defaults that were disabled.
    foreach(i -> i in types || push!(plugins, !defaults[i]), 1:ndefaults)
    return plugins
end

# Call the default prompt method even if a specialized one exists.
fallback_prompt(T::Type, name::Symbol) = prompt(Template, T, Val(name), nothing)


text:
# include("Configuration.jl")
# include("Lattice.jl")

###############################################################
#################### Free fermion orbitals ####################
###############################################################

immutable FreeFermionOrbitals{T}
    all::Matrix{T}
    filled::Matrix{T}

    energies::Vector{Float64}
    Fermi_energy::Float64

    no_degeneracies_at_FS::Bool
    degeneracy_tolerance::Float64

    function FreeFermionOrbitals(Hami::Matrix{T}, Nfilled::Int, deg_tol::Float64 = 1e-12)
        @assert ishermitian(Hami)
        @assert Nfilled <= size(Hami)[1]
        eig_data = eigfact(Hami)
        if Nfilled == size(Hami)[1]
            no_deg = true
            warn("Filling up all possible states. Have fun.")
        else
            # printlnf(eig_data[:values][Nfilled+1] - eig_data[:values][Nfilled] > deg_tol)
            if eig_data[:values][Nfilled+1] - eig_data[:values][Nfilled] > deg_tol
                no_deg = true
            else
                no_deg = false
                warn("Degeneracies exist at the Fermi surface. Proceed with caution!")
            end
        end
        return new(eig_data[:vectors],
                   eig_data[:vectors][:, 1:Nfilled],
                   eig_data[:values],
                   eig_data[:values][Nfilled],
                   no_deg, deg_tol)
    end
end

FreeFermionOrbitals(Hami::Matrix, Nfilled::Int) = FreeFermionOrbitals{eltype(Hami)}(Hami, Nfilled)

get_Nsites(orbs::FreeFermionOrbitals) = size(orbs.filled)[1]
get_Nfilled(orbs::FreeFermionOrbitals) = size(orbs.filled)[2]
no_degeneracies_at_FS(orbs::FreeFermionOrbitals) = orbs.no_degeneracies_at_FS

function fill_detmat(orbs::FreeFermionOrbitals, config::LatticeParticlesConfiguration, species::Int = 1)
    @assert get_Nsites(orbs) == get_Nsites(config)
    @assert species <= get_Nspecies(config)
    @assert get_Nfilled(orbs) == get_Nparticles(config)[species]

    detmat = transpose(orbs.filled[config.particle_positions[species], :])

    return detmat
end

function fill_detmats{T}(orbs::FreeFermionOrbitals{T}, config::LatticeParticlesConfiguration)
    detmats = Matrix{T}[]
    for sp in 1:get_Nspecies(config)
        push!(detmats, fill_detmat(orbs, config, sp))
    end
    return detmats
end

function columns_for_detmat(orbs::FreeFermionOrbitals, positions::Vector{Int})
    @assert all(positions .<= get_Nsites(orbs))
    return transpose(orbs.filled[positions, :])
end

columns_for_detmat(orbs::FreeFermionOrbitals, position::Int) = columns_for_detmat(orbs, [position])

# construct free fermion Green's function as the single-particle density matrix corresponding to a filled Fermi sea
function single_particle_density_matrix{T}(orbs::FreeFermionOrbitals{T})
    # spdm = zeros(T, (get_Nsites(orbs), get_Nsites(orbs)))
    # for i in 1:get_Nfilled(orbs)
    #     spdm = spdm + orbs.filled[:, i] * orbs.filled[:, i]'
    # end
    # @assert ishermitian(spdm)
    # NB: the above lines are slow, but emphasize how the SPDM is defined
    spdm = orbs.filled * orbs.filled' # this is obviously much faster
    return transpose(spdm) # NB: transpose because spdm_{ij} = <c_j^dagger c_i>; we want to return <c_i^dagger c_j>
end

# local density in real space
local_density(orbs::FreeFermionOrbitals) = real(diag(single_particle_density_matrix(orbs)))

# compute density-density correlation function <n_i n_j> using Wick's theorem
function density_density_correlation_function(orbs::FreeFermionOrbitals)
    spdm = single_particle_density_matrix(orbs)
    ni = local_density(orbs)
    return real(-abs(spdm).^2 + diagm(ni) + ni * ni.')
end

# core function in AnalysisTools.jl
function local_particle_number_variance(orbs::FreeFermionOrbitals, sites_per_region::Int)
    ninj = density_density_correlation_function(orbs)
    ni = local_density(orbs)
    return local_particle_number_variance(ninj, ni, sites_per_region)
end

# core function in AnalysisTools.jl
function bipartite_number_fluctuations(orbs::FreeFermionOrbitals, sites_added_per_cut::Int)
    ninj = density_density_correlation_function(orbs)
    ni = local_density(orbs)
    return bipartite_number_fluctuations(ninj, ni, sites_added_per_cut)
end

# calculate von Neumann entanglement entropy using correlation matrix technique for many sequential cuts through a system
function bipartite_vonNeumann_entanglement_entropy(orbs::FreeFermionOrbitals, sites_added_per_cut::Int) # = garbage
    @assert sites_added_per_cut < get_Nsites(orbs)

    spdm = single_particle_density_matrix(orbs)

    if mod(get_Nsites(orbs), sites_added_per_cut) == 0
        numcuts = convert(Int, get_Nsites(orbs)/sites_added_per_cut - 1)
    else
        numcuts = convert(Int, floor(get_Nsites(orbs)/sites_added_per_cut))
    end

    S1 = zeros(Float64, (numcuts, 2))
    for x in 1:numcuts
        sites_in_partition = sites_added_per_cut * x
        all_evals = eigfact(spdm[1:sites_in_partition, 1:sites_in_partition])[:values]
        good_evals = all_evals[find(0. .< all_evals .< 1.)]
        S1[x, 1] = sites_in_partition
        S1[x, 2] = -(good_evals.' * log(good_evals) + (1. - good_evals.') * log(1. - good_evals))[1]
    end

    return S1
end

function Renyi_formula(cs::Vector{Float64}, alpha::Int)
    @assert minimum(cs) > 0. && maximum(cs) < 1.

    if alpha == 1
        # S = -((1. - cs.') * log(1. - cs) + cs.' * log(cs))[1]
        S = -((1. - cs.') * log1p(-cs) + cs.' * log(cs))[1]
    elseif alpha == 2
        S = -sum(log((1. - cs).^2 + cs.^2))
        # S = -sum(log1p(2 * cs .* (cs - 1.)))
    elseif alpha >= 3
        S = 1./(1. - alpha) * sum(log((1. - cs).^alpha + cs.^alpha))
    else
        error("Renyi_formula requires alpha > 0.")
    end

    return S
end

function bipartite_entanglement_entropy(orbs::FreeFermionOrbitals, subsystem::Vector{Int}, alpha::Int)
    @assert 0 <= length(subsystem) <= get_Nsites(orbs)
    if length(subsystem) == 0 || length(subsystem) == get_Nsites(orbs)
        return 0.
    end
    @assert minimum(subsystem) >= 1 && maximum(subsystem) <= get_Nsites(orbs)

    C = single_particle_density_matrix(orbs)[subsystem, subsystem] # = correlation matrix
    cs = eigfact(C)[:values]
    cs = cs[find(0. .< cs .< 1.)] # good evals

    return Renyi_formula(cs, alpha)
end

function bipartite_entanglement_entropy(spdm::Matrix, subsystem::Vector{Int}, alpha::Int)
    @assert ishermitian(spdm)
    Nsites = size(spdm)[1]
    @assert 0 <= length(subsystem) <= Nsites
    if length(subsystem) == 0 || length(subsystem) == Nsites
        return 0.
    end
    @assert minimum(subsystem) >= 1 && maximum(subsystem) <= Nsites

    C = spdm[subsystem, subsystem]
    cs = eigfact(C)[:values]
    cs = cs[find(0. .< cs .< 1.)] # good evals

    return Renyi_formula(cs, alpha)
end

##########################################################
#################### Spinful BCS data ####################
##########################################################

immutable SpinfulBCS{T}
    evecs::Matrix{T}
    U::Matrix{T}
    V::Matrix{T}

    energies::Vector{Float64}

    function SpinfulBCS(t::Matrix, Delta::Matrix)
        # both t and Delta should include on-site terms, e.g., t_{ij} = ttilde_{ij} = t_{ij} + mu_i * krondelta_{ij} from notes
        @assert ishermitian(t)
        @assert issym(Delta)
        @assert size(t)[1] == size(Delta)[1]

        Nsites = size(t)[1]

        H = zeros(T, (2Nsites, 2Nsites))
        H[1:Nsites, 1:Nsites] = -t
        H[Nsites+1:2Nsites, Nsites+1:2Nsites] = conj(t) # = t.'
        H[1:Nsites, Nsites+1:2Nsites] = Delta
        H[Nsites+1:2Nsites, 1:Nsites] = conj(Delta)
        @assert ishermitian(H)

        eig_data = eigfact(H)
        energies = eig_data[:values]

        # check that energies come in plus-minus pairs
        @assert maxabs(flipdim(energies[1:Nsites], 1) + energies[Nsites+1:2Nsites]) < 1e-12

        evecs = eig_data[:vectors]

        # U and V are defined by the positive-energy eigenvectors (see notes)
        U = evecs[1:Nsites, Nsites+1:2Nsites]
        V = evecs[Nsites+1:2Nsites, Nsites+1:2Nsites]

        # check that (v, -u)^* are indeed the correct negative-energy eigenvectors
        for j in 1:Nsites
            evec_test = vcat(conj(V[:, j]), -conj(U[:, j]))
            @assert (maxabs(H * evec_test - energies[Nsites-j+1] * evec_test) < 1e-12) "$(j) is bad"
        end

        return new(evecs, U, V, energies)
    end
end

SpinfulBCS(t::Matrix, Delta::Matrix) = SpinfulBCS{promote_type(eltype(t), eltype(Delta))}(t, Delta)

get_Nsites(bcs::SpinfulBCS) = size(bcs.U)[1]

function fill_detmat_for_SL{T}(bcs::SpinfulBCS{T}, config::LatticeParticlesConfiguration)
    @assert get_Nsites(bcs) == get_Nsites(config)
    Nsites = get_Nsites(bcs)
    @assert Nsites/2 == get_Nparticles(config)[1] # hard-core bosons are at half-filling

    detmat = zeros(T, (Nsites, Nsites))
    detmat[:, 1:2:Nsites-1] = transpose(bcs.evecs[config.particle_positions[1], 1:Nsites])
    detmat[:, 2:2:Nsites] = transpose(bcs.evecs[Nsites + config.particle_positions[1], 1:Nsites])

    return detmat
end

function columns_for_detmat_for_SL{T}(bcs::SpinfulBCS{T}, positions::Vector{Int})
    @assert all(positions .<= get_Nsites(bcs))
    Nsites = get_Nsites(bcs)

    cols = zeros(T, (Nsites, 2*length(positions)))
    cols[:, 1:2:2*length(positions)-1] = transpose(bcs.evecs[positions, 1:Nsites])
    cols[:, 2:2:2*length(positions)] = transpose(bcs.evecs[Nsites + positions, 1:Nsites])

    return cols
end

columns_for_detmat_for_SL(bcs::SpinfulBCS, position::Int) = columns_for_detmat_for_SL(bcs, [position])

#######################################################
#################### Hopping Hamis ####################
#######################################################

function zigzag_strip_hopping_Hami(t1::Complex128, t2::Complex128, mu::Float64, BCphase::Complex128, Lx::Int)
    xhat = unitx(Lx, 1)
    Nsites = Lx
    tmat = zeros(Complex128, (Nsites, Nsites))

    # bulk sites t1
    for i in 1:Nsites-1
        tmat[i, (Site2D(i, Lx, 1) + xhat).index] = t1
    end
    # end site t1
    tmat[Nsites, (Site2D(Nsites, Lx, 1) + xhat).index] = BCphase * t1

    # bulk sites t2
    for i in 1:Nsites-2
        tmat[i, (Site2D(i, Lx, 1) + 2xhat).index] = t2
    end
    # end sites t2
    for i in Nsites-1:Nsites
        tmat[i, (Site2D(i, Lx, 1) + 2xhat).index] = BCphase * t2
    end

    tmat = tmat + tmat'

    for i in 1:Nsites
        tmat[i, i] = mu
    end

    H = -tmat

    # convert to a real matrix if we can
    return imag(H) == zeros(size(H)...) ? real(H) : H
end

# the argument n::Int determines the flux per triangle consistent with the kx preserving gauge
function uniform_flux_triangular_ladder_kx_gauge_hopping_Hami(tx::Float64, ty::Float64, td::Float64, n::Int, BCxphase::Number, Lx::Int, Ly::Int)
    # BCxphase should either be of unit modulus or zero (for OBCs)
    # I think having BCxphase::Number in the argument list is good (so it can be either real or complex)
    # @assert abs(BCxphase) == 1 || BCxphase == 0
    # Eh, maybe not.  Perhaps BCxphase should just be renamed BCxfac or something..

    xhat = unitx(Lx, Ly)
    yhat = unity(Lx, Ly)
    Nsites = Ly * Lx
    phi = pi * n / Ly
    tmat = zeros(Complex128, (Nsites, Nsites))

    for i in 1:Nsites
        if Site2D(i, Lx, Ly).x < Lx-1
            tmat[i, (Site2D(i, Lx, Ly) + xhat).index] = tx * exp(-im * 2 * Site2D(i, Lx, Ly).y * phi)
            tmat[i, (Site2D(i, Lx, Ly) + xhat + yhat).index] = td * exp(-im * (2 * Site2D(i, Lx, Ly).y + 1) * phi)
            tmat[i, (Site2D(i, Lx, Ly) + yhat).index] = ty
        elseif Site2D(i, Lx, Ly).x == Lx-1
            # terms corresponding to the next two lines cross the vertical boundary at the end of the sample
            # NB: collectively, they don't change the flux piercing the corresponding plaquettes
            tmat[i, (Site2D(i, Lx, Ly) + xhat).index] = BCxphase * tx * exp(-im * 2 * Site2D(i, Lx, Ly).y * phi)
            tmat[i, (Site2D(i, Lx, Ly) + xhat + yhat).index] = BCxphase * td * exp(-im * (2 * Site2D(i, Lx, Ly).y + 1) * phi)
            tmat[i, (Site2D(i, Lx, Ly) + yhat).index] = ty
        else
            @assert false
        end
    end

    tmat = tmat + tmat'
    H = -tmat

    # convert to a real matrix if we can
    return imag(H) == zeros(size(H)...) ? real(H) : H
end

# allows for different BCphase factors in both the x and y directions
function uniform_flux_triangular_lattice_kx_gauge_hopping_Hami(tx::Float64, ty::Float64, td::Float64, n::Int,
                                                               BCxphase::Number, BCyphase::Number, Lx::Int, Ly::Int)
# FIXME: build it
end

function kagome_strip_hopping_Hami(tleg::Float64, tcross::Float64, muleg::Float64, mumid::Float64, phi::Float64, BCxfac::Number, Lx::Int)
    @assert abs(BCxfac) == 1 || BCxfac == 0

    Ly = 3
    xhat = unitx(Lx, Ly)
    yhat = unity(Lx, Ly)
    Nsites = Ly * Lx
    tmat = zeros(Complex128, (Nsites, Nsites))

    for i in 1:Nsites
        if Site2D(i, Lx, Ly).x < Lx-1
            if Site2D(i, Lx, Ly).y == 0
                tmat[i, (Site2D(i, Lx, Ly) + xhat).index] = tleg * exp(-im * phi)
                tmat[i, (Site2D(i, Lx, Ly) + xhat + yhat).index] = tcross
                tmat[i, (Site2D(i, Lx, Ly) + yhat).index] = tcross * exp(im * phi)
            elseif Site2D(i, Lx, Ly).y == 1
                tmat[i, (Site2D(i, Lx, Ly) + yhat).index] = tcross
            elseif Site2D(i, Lx, Ly).y == 2
                tmat[i, (Site2D(i, Lx, Ly) + xhat).index] = tleg
                tmat[i, (Site2D(i, Lx, Ly) + xhat - yhat).index] = tcross
            else
                @assert false
            end
        elseif Site2D(i, Lx, Ly).x == Lx-1
            if Site2D(i, Lx, Ly).y == 0
                tmat[i, (Site2D(i, Lx, Ly) + xhat).index] = BCxfac * tleg * exp(-im * phi)
                tmat[i, (Site2D(i, Lx, Ly) + xhat + yhat).index] = BCxfac * tcross
                tmat[i, (Site2D(i, Lx, Ly) + yhat).index] = tcross * exp(im * phi)
            elseif Site2D(i, Lx, Ly).y == 1
                tmat[i, (Site2D(i, Lx, Ly) + yhat).index] = tcross
            elseif Site2D(i, Lx, Ly).y == 2
                tmat[i, (Site2D(i, Lx, Ly) + xhat).index] = BCxfac * tleg
                tmat[i, (Site2D(i, Lx, Ly) + xhat - yhat).index] = BCxfac * tcross
            else
                @assert false
            end
        else
            @assert false
        end
    end

    tmat = tmat + tmat'

    for i in 1:Nsites
        if Site2D(i, Lx, Ly).y == 0
            tmat[i, i] = muleg
        elseif Site2D(i, Lx, Ly).y == 1
            tmat[i, i] = mumid
        elseif Site2D(i, Lx, Ly).y == 2
            tmat[i, i] = muleg
        else
            @assert false
        end
    end

    H = -tmat

    # convert to a real matrix if we can
    return imag(H) == zeros(size(H)...) ? real(H) : H
end

function no_flux_square_lattice_hopping_Hami(tx::Float64, ty::Float64, BCxfac::Number, BCyfac::Number, Lx::Int, Ly::Int)
    @assert abs(BCxfac) == 1 || BCxfac == 0
    @assert abs(BCyfac) == 1 || BCyfac == 0

    xhat = unitx(Lx, Ly)
    yhat = unity(Lx, Ly)
    Nsites = Ly * Lx
    tmat = zeros(Complex128, (Nsites, Nsites))

    for i in 1:Nsites
        if Site2D(i, Lx, Ly).x < Lx-1 && Site2D(i, Lx, Ly).y < Ly-1
            tmat[i, (Site2D(i, Lx, Ly) + xhat).index] = tx
            tmat[i, (Site2D(i, Lx, Ly) + yhat).index] = ty
        elseif Site2D(i, Lx, Ly).x < Lx-1 && Site2D(i, Lx, Ly).y == Ly-1
            tmat[i, (Site2D(i, Lx, Ly) + xhat).index] = tx
            tmat[i, (Site2D(i, Lx, Ly) + yhat).index] = BCyfac * ty
        elseif Site2D(i, Lx, Ly).x == Lx-1 && Site2D(i, Lx, Ly).y < Ly-1
            tmat[i, (Site2D(i, Lx, Ly) + xhat).index] = BCxfac * tx
            tmat[i, (Site2D(i, Lx, Ly) + yhat).index] = ty
        elseif Site2D(i, Lx, Ly).x == Lx-1 && Site2D(i, Lx, Ly).y == Ly-1
            tmat[i, (Site2D(i, Lx, Ly) + xhat).index] = BCxfac * tx
            tmat[i, (Site2D(i, Lx, Ly) + yhat).index] = BCyfac * ty
        else
            @assert false
        end
    end

    tmat = tmat + tmat'
    H = -tmat

    # convert to a real matrix if we can
    return imag(H) == zeros(size(H)...) ? real(H) : H
end

ext: jl


text:
module fancy

# Functions
export say_hello
say_hello() = println("Hello!")

export calculate_sum
function calculate_sum(x,y)
    x+y
end

# Values
export prices
prices = Dict([("Horse", 1000), ("Dog", 200)])

# Don't export this function
function unexported_function()
    println("I'm not exported")
end

end


text:
# this file is the run script for the covid19ode module.
# all plotting of results and analysis goes here.
# in julia, all have to do is include this file with the correct environment set

using Gnuplot
using covid19ode

function hello()
    println("hello world")
end


text:
function EM_converge!(hmm_jobs::RemoteChannel, output_hmms::RemoteChannel, no_models::Integer; load_config::LoadConfig=LoadConfig(1:typemax(Int64)-1, 0:typemax(Int64)-1),  EM_func::Function=linear_step, delta_thresh=1e-3, max_iterates=5000, verbose=false)
    while isready(hmm_jobs)
        workerid = myid()
        jobid, start_iterate, hmm, job_norm, observations = load_balancer(no_models, hmm_jobs, load_config)
        jobid == 0 && break #no valid job for this worker according to load_table entry

        start_iterate > max_iterates - 1 && throw(ArgumentError("BHMM chain $jobid is already longer ($start_iterate iterates) than specified max_iterates!"))
        curr_iterate = start_iterate

        #mask calculations here rather than ms_mle_step to prevent recalculation every iterate
        #build array of observation lengths

        obs_lengths = [findfirst(iszero,observations[o,:])-1 for o in 1:size(observations,1)] #mask calculations here rather than mle_step to prevent recalculation every iterate
        EM_func==bw_step && (observations=transpose(observations))

        start_iterate == 1 && put!(output_hmms, (workerid, jobid, curr_iterate, hmm, job_norm, 0.0, false, 0.0)); #on the first iterate return the initial BHMM immediately
        verbose && @info "Fitting BHMM on Wk $workerid, start iterate $start_iterate, $jobid with $(size(hmm)[1]) states and $(size(hmm)[2]) symbols..."

        curr_iterate += 1

        start=time()
        new_hmm, last_norm = EM_func(hmm, observations, obs_lengths)
        delta = abs(lps(job_norm, -last_norm))
        put!(output_hmms, (workerid, jobid, curr_iterate, new_hmm, last_norm, delta, false, time()-start))

        for i in curr_iterate:max_iterates
            start=time()
            new_hmm, norm = EM_func(new_hmm, observations, obs_lengths)
            curr_iterate += 1
            delta = abs(lps(norm, -last_norm))
            if delta < delta_thresh
                put!(output_hmms, (workerid, jobid, curr_iterate, new_hmm, norm, delta, true, time()-start))
                verbose && @info "$jobid converged after $(curr_iterate-1) EM steps"
                break
            else
                put!(output_hmms, (workerid, jobid, curr_iterate, new_hmm, norm, delta, false, time()-start))
                verbose && @info "$jobid EM step $(curr_iterate-1) delta $delta"
                last_norm = norm
            end
        end
    end
end


text:
# This definitely needs cleaning
using OrdinaryDiffEq, DiffEqDevTools, Test, Random
using DiffEqProblemLibrary.ODEProblemLibrary: importodeproblems; importodeproblems()
import DiffEqProblemLibrary.ODEProblemLibrary: prob_ode_linear, prob_ode_2Dlinear
probArr = Vector{ODEProblem}(undef, 2)
probArr[1] = prob_ode_linear

probArr[2] = prob_ode_2Dlinear
Random.seed!(100)
## Convergence Testing
dts = 1 .//2 .^(8:-1:4)
dts1 = 1 .//2 .^(9:-1:5)
testTol = 0.2

for i = 1:2
  global dts
  prob = probArr[i]
  sim = test_convergence(dts,prob,Euler())
  @test abs(sim.𝒪est[:final]-1) < testTol
  sim2 = test_convergence(dts,prob,Heun())
  @test abs(sim2.𝒪est[:l∞]-2) < testTol
  sim2 = test_convergence(dts,prob,Ralston())
  @test abs(sim2.𝒪est[:l∞]-2) < testTol
  sim2 = test_convergence(dts,prob,Midpoint())
  @test abs(sim2.𝒪est[:l∞]-2) < testTol
  sim3 = test_convergence(dts,prob,RK4())
  @test abs(sim3.𝒪est[:l∞]-4) < testTol
  sim4 = test_convergence(dts,prob,BS3())
  @test abs(sim4.𝒪est[:l2]-3) < testTol
  alg = CarpenterKennedy2N54()
  @test abs(test_convergence(dts,prob,alg).𝒪est[:l∞]-OrdinaryDiffEq.alg_order(alg)) < testTol
  sim5 = test_convergence(dts, prob, AB3())
  @test abs(sim5.𝒪est[:l2]-3) < testTol
  sim6 = test_convergence(dts,prob,ABM32())
  @test abs(sim6.𝒪est[:l2]-3) < testTol
  sim7 = test_convergence(dts, prob, AB4())
  @test abs(sim7.𝒪est[:l2]-4) < testTol
  sim8 = test_convergence(dts1,prob,ABM43())  #using dts1 due to floating point error in convergence test
  @test abs(sim8.𝒪est[:l2]-4) < testTol
  sim9 = test_convergence(dts,prob,AB5())
  @test abs(sim9.𝒪est[:l2]-5) < testTol
  sim10 = test_convergence(dts,prob,ABM54())
  @test abs(sim10.𝒪est[:l2]-5) < testTol
  sim101 = test_convergence(dts,prob,VCAB3())
  @test abs(sim101.𝒪est[:l2]-3) < testTol
  sim102 = test_convergence(dts,prob,VCAB4())
  @test abs(sim102.𝒪est[:l2]-4) < testTol
  sim103 = test_convergence(dts,prob,VCAB5())
  @test abs(sim103.𝒪est[:l2]-5) < testTol
  sim104 = test_convergence(dts,prob,VCABM3())
  @test abs(sim104.𝒪est[:l2]-3) < testTol
  sim105 = test_convergence(dts,prob,VCABM4())
  @test abs(sim105.𝒪est[:l2]-4) < testTol
  sim106 = test_convergence(dts,prob,VCABM5())
  @test abs(sim106.𝒪est[:l2]-5) < testTol
  sim160 = test_convergence(dts,prob,Anas5(w=2))
  @test abs(sim160.𝒪est[:l2]-4) < 2*testTol

  println("Stiff Solvers")

  dts = 1 .//2 .^(9:-1:5)

  sim11 = test_convergence(dts,prob,ImplicitEuler(extrapolant = :linear))
  @test abs(sim11.𝒪est[:final]-1) < testTol

  sim112 = test_convergence(dts,prob,ImplicitEuler(nonlinsolve = NLFunctional()))
  @test abs(sim112.𝒪est[:final]-1) < testTol

  sim12 = test_convergence(dts,prob,
          GenericImplicitEuler(nlsolve=OrdinaryDiffEq.NLSOLVEJL_SETUP(autodiff=true)))
  @test abs(sim12.𝒪est[:final]-1) < testTol

  sim122 = test_convergence(dts,prob,
           GenericImplicitEuler(nlsolve=OrdinaryDiffEq.NLSOLVEJL_SETUP(autodiff=false)))

  sim13 = test_convergence(dts,prob,ImplicitMidpoint())
  @test abs(sim13.𝒪est[:final]-2) < testTol

  sim132 = test_convergence(dts,prob,ImplicitMidpoint(nonlinsolve = NLFunctional()))
  @test abs(sim132.𝒪est[:final]-2) < testTol

  sim13 = test_convergence(dts,prob,Trapezoid())
  @test abs(sim13.𝒪est[:final]-2) < testTol

  sim133 = test_convergence(dts,prob,Trapezoid(nonlinsolve = NLFunctional()))
  @test abs(sim133.𝒪est[:final]-2) < testTol

  sim14 = test_convergence(dts,prob,
          GenericTrapezoid(nlsolve=OrdinaryDiffEq.NLSOLVEJL_SETUP(autodiff=true)))
  @test abs(sim14.𝒪est[:final]-2) < testTol

  sim142 = test_convergence(dts,prob,
           GenericTrapezoid(nlsolve=OrdinaryDiffEq.NLSOLVEJL_SETUP(autodiff=false)))
  @test abs(sim142.𝒪est[:final]-2) < testTol

  sim14 = test_convergence(dts,prob,TRBDF2())
  @test abs(sim14.𝒪est[:final]-2) < testTol

  sim152 = test_convergence(dts,prob,TRBDF2(autodiff=false))
  @test abs(sim152.𝒪est[:final]-2) < testTol+0.1

  sim15 = test_convergence(dts,prob,SDIRK2())
  @test abs(sim15.𝒪est[:final]-2) < testTol

  sim152 = test_convergence(dts,prob,SSPSDIRK2())
  @test abs(sim152.𝒪est[:final]-2) < testTol

  sim16 = test_convergence(dts,prob,Kvaerno3())
  @test abs(sim16.𝒪est[:final]-3) < testTol

  sim162 = test_convergence(dts,prob,Kvaerno3(nonlinsolve = NLFunctional()))
  @test abs(sim162.𝒪est[:final]-3) < testTol

  sim17 = test_convergence(dts,prob,KenCarp3())
  @test abs(sim17.𝒪est[:final]-3) < testTol

  #####################################
  # BDF
  #####################################

  sim = test_convergence(dts,prob,ABDF2(extrapolant = :linear))
  @test abs(sim.𝒪est[:final]-2) < testTol
  @test abs(sim.𝒪est[:l2]-2) < testTol
  @test abs(sim.𝒪est[:l∞]-2) < testTol

  sim = test_convergence(dts,prob,ABDF2(nonlinsolve = NLFunctional()))
  @test abs(sim.𝒪est[:final]-2) < testTol
  @test abs(sim.𝒪est[:l2]-2) < testTol
  @test abs(sim.𝒪est[:l∞]-2) < testTol

  # QBDF
  sim = test_convergence(dts,prob,QBDF1())
  @test abs(sim.𝒪est[:final]-1) < testTol
  @test abs(sim.𝒪est[:l2]-1) < testTol
  @test abs(sim.𝒪est[:l∞]-1) < testTol

  sim = test_convergence(dts,prob,QBDF2())
  @test abs(sim.𝒪est[:final]-2) < testTol
  @test abs(sim.𝒪est[:l2]-2) < testTol
  @test abs(sim.𝒪est[:l∞]-2) < testTol

  # QNDF
  sim = test_convergence(dts,prob,QNDF1())
  @test abs(sim.𝒪est[:final]-1) < testTol
  @test abs(sim.𝒪est[:l2]-1) < testTol
  @test abs(sim.𝒪est[:l∞]-1) < testTol

  sim = test_convergence(dts,prob,QNDF2())
  @test abs(sim.𝒪est[:final]-2) < testTol
  @test abs(sim.𝒪est[:l2]-2) < testTol
  @test abs(sim.𝒪est[:l∞]-2) < testTol

  sim = test_convergence(dts,prob,QNDF2(nonlinsolve = NLFunctional()))
  @test abs(sim.𝒪est[:final]-2) < testTol
  @test abs(sim.𝒪est[:l2]-2) < testTol
  @test abs(sim.𝒪est[:l∞]-2) < testTol

  dts = 1 .//2 .^(7:-1:4)
  println("Higher Order")

  sim18 = test_convergence(dts,prob,Cash4())
  @test abs(sim18.𝒪est[:final]-4) < testTol

  sim19 = test_convergence(dts,prob,Hairer4())
  @test abs(sim19.𝒪est[:final]-4) < testTol

  sim110 = test_convergence(dts,prob,Hairer42())
  @test abs(sim110.𝒪est[:final]-4) < testTol

  sim111 = test_convergence(dts,prob,Kvaerno4())
  @test abs(sim111.𝒪est[:final]-4) < testTol

  sim112 = test_convergence(dts,prob,KenCarp4())
  @test abs(sim112.𝒪est[:final]-4) < testTol

  sim113 = test_convergence(dts,prob,Kvaerno5())
  @test abs(sim113.𝒪est[:final]-5) < testTol

  sim114 = test_convergence(dts,prob,KenCarp5())
  @test abs(sim114.𝒪est[:final]-5) < testTol

  sim115 = test_convergence(dts,prob,KenCarp5(nonlinsolve = NLFunctional()))
  @test_broken abs(sim115.𝒪est[:final]-5) < testTol
end


text:
module APIUtils

using ..VectorEngine

using Libdl

# helpers that facilitate working with VE APIs
include("call.jl")
include("enum.jl")

end


text:
RecipesBase.@recipe function f(ma::LinearEstimator)
    RecipesBase.@series begin
        label --> pretty_label(ma) #TODO, make more generic..
        seriestype --> :line
        xlabel --> :x
        ma.Q.marginal_grid,  ma.Q.Q .+ ma.Q.Qo
    end
end


text:
"""
This script shows how to assemble full 360 degree radar sweeps from the DFG object
built with RExFeed.jl
"""

using DistributedFactorGraphs
using IncrementalInference, RoME
using JSON2
using LinearAlgebra

##

# Where to fetch data
# dfgDataFolder = ENV["HOME"]*"/data/rex";
dfgDataFolder = "/tmp/caesar/rex"

# Load the graph
fg = loadDFG("$dfgDataFolder/dfg")

# add the datastore locations
ds = FolderStore{Vector{UInt8}}(:radar, "$dfgDataFolder/data/radar")
addBlobStore!(fg, ds)

ds = FolderStore{Vector{UInt8}}(:gps_fix, "$dfgDataFolder/data/gps")
addBlobStore!(fg, ds)

ds = FolderStore{Vector{UInt8}}(:lidar, "$dfgDataFolder/data/lidar")
addBlobStore!(fg, ds)

##

# Check what's in it
ls(fg)
lsf(fg)  # there should be no factors the first time a session is loaded

# How about radar data entries?
count(v -> :RADAR in listDataEntries(v), getVariables(fg))

# may be too intense as it is fetching ALL data too?
allRadarVariables = filter(v -> :RADAR in listDataEntries(v), getVariables(fg));

allRadarVariables = sortDFG(allRadarVariables);

function fetchRadarMsg(var::DFGVariable)
    entry,rawData = getData(fg, var.label, :RADAR)
    # raw data is json-encoded; this decoding should happen inside getBigData?
    return JSON2.read(IOBuffer(rawData))
end

msg = fetchRadarMsg(allRadarVariables[1])

##

# function azimuth(msg::NamedTuple)
function azimuth(msg)
    """
    Assemble azimuth vector

    This function assembles a vector of azimuth angles for the given ping. It
    handles both the conversion from [0,360) to [-pi,pi), and the wrap around
    scenario when msg.angle_start > msg.angle_end (or, equivalently,
    msg.angle_increment != 0.17578125).

    """
    if msg.angle_start > msg.angle_end
        # wrap around; increments become messy
        delta = (360 - msg.angle_start) + (msg.angle_end)
        angle_step = delta/msg.n
        a = Vector(range(msg.angle_start, msg.angle_start+delta,length = msg.n))
        a[a.>360].-=360.
        #   sa = [Vector(range(msg.angle_start, 360.0, step = 0.17578125)); Vector(reverse(range(msg.angle_end, 0.,step=-0.17578125)))]
        # use reverse to ensure we stop at the correct angle
    else
        a = range(msg.angle_start, msg.angle_end, length=msg.n)
    end

    # @show a
    a = deg2rad.(Vector(a))   # to radians
    a[a.>pi]=a[a.>pi].-(2*pi) # normalize to [-pi,pi)

    return (a)
end

function nearest(v::Vector{<:Real}, q::Real)::Int64
    return findmin(abs.(v.-q))[2]
end

function tocartsector(ping_polar::Array{Any,2},r::Vector{Float64}, a::Vector{Float64}, res::Float64 )
    """
    res - resolution, in m/px
    """

    vx = [r[1]*cos(a[1]), r[1]*cos(a[end]), r[end]*cos(a[1]), r[end]*cos(a[end])];
    vy = [r[1]*sin(a[1]), r[1]*sin(a[end]), r[end]*sin(a[1]), r[end]*sin(a[end])];
    xmin =  minimum(vx)
    xmax = maximum(vx)
    nx = Int64(round((xmax-xmin)/res))
    x = Vector(range(xmin, xmax, length=nx))

    ymin = minimum(vy)
    ymax = maximum(vy)
    ny = Int64(round((ymax-ymin)/res))
    y = Vector(range(ymin, ymax, length=ny))

    ping_cart = zeros(nx, ny)

    for i = 1:nx
        for j=1:ny
            rq = norm([x[i],y[j]])
            if rq > r[end] || rq < r[1]
                continue
            end

            aq = atan(y[j],x[i])
            if aq > a[end] || aq < a[1]
                continue
            end

            # idx_r = Int(ceil((rq-r[1])/(r[end]-r[1])))
            # idx_a = Int(ceil((aq-a[1])/(a[end]-a[1]))) # dangerous!
            # pvt: turns out that using nearest, while arguably more readable, results in a significant performance hit!
            #############
            # w/ nearest:
            # memory estimate:  1.28 GiB
            # allocs estimate:  1386472
            # --------------
            # minimum time:     797.162 ms (22.86% GC)
            # median time:      960.964 ms (22.70% GC)
            # mean time:        962.183 ms (22.91% GC)
            # maximum time:     1.208 s (22.71% GC)
            # --------------
            # samples:          6
            # evals/sample:     1
            ##############
            # w/o nearest:
            # memory estimate:  72.06 MiB
            # allocs estimate:  922543
            # --------------
            # minimum time:     76.095 ms (19.39% GC)
            # median time:      78.501 ms (20.36% GC)
            # mean time:        94.679 ms (33.84% GC)
            # maximum time:     227.771 ms (70.20% GC)
            # --------------
            # samples:          53
            # evals/sample:     1
            idx_r =  nearest(r,  rq)
            idx_a = nearest(a, aq)
            ping_cart[i,j] = ping_polar[idx_r, idx_a]
        end
    end
    return ping_cart
end

function tocart(ping_polar::Array{Any,2},r::Vector{<:Real}, a::Vector{<:Real}, res::Real )
    """
    res - resolution, in m/px
    """

    n = Int64(round((2*r[end]/res)))
    x = Vector(range(-r[end], r[end], length=n))
    y = Vector(range(-r[end], r[end], length=n))
    ping_cart = zeros(n, n)

    for i = 1:n
        for j=1:n
            rq = norm([x[i],y[j]])
            if rq > r[end] || rq < r[1]
                continue
            end

            aq = atan(y[j],x[i])
            if aq > a[end] || aq < a[1]
                continue
            end

            idx_r =  nearest(r,  rq)
            idx_a = nearest(a, aq)
            ping_cart[i,j] = ping_polar[idx_r, idx_a]
        end
    end
    return ping_cart
end

# function assembleping(msg::NamedTuple, resolution::Float64)
function assembleping(msg, resolution::Float64)
    """
    Assemble polar and cartesian pings from radar message.
    """
    ping_polar= reshape(msg.radar_data, (msg.m, msg.n))

    # assemble azimuth and range vectors
    a = azimuth(msg)
    r = Vector(range(msg.range_min, msg.range_max, length=msg.m ))

    # option a: tight cartesian around the sector
    # ping_cart = tocartsector(ping_polar, r, a, resolution)

    # option b: full 360
    ping_cart = tocart(ping_polar, r, a, resolution)

    return (ping_polar, ping_cart)
end

# msg = alldata[1]
# (polar, cart) = assembleping(msg, 2.0);
# fullsweep = round.(UInt8, cart./255)
# imi = imshow(fullsweep)

# for i = 1:20
#     @sync begin
#         @async begin
#             println(i)
#             if mod(i, 10) == 0
#                 fullsweep .&= 0x00
#             end
#             # sleep(0.1)
#         end
#         @async begin
#             msg = alldata[i]
#             # set resolution to 1m/px
#             (polar, cart) = assembleping(msg, 2.0);
#             fullsweep .|= round.(UInt8, cart./255)
#             imshow!(imi["gui"]["canvas"], fullsweep)
#         end
#     end
# end

##

# set resolution to 1m/px
(polar, cart) = assembleping(msg, 1.0);

using ImageMagick
using Images, ImageView
using ImageFiltering
# imshow(polar)
# imshow(cart)

##

output_dir = joinpath(dfgDataFolder,"pings")
if (!isdir(output_dir))
    mkdir(output_dir)
end
save(joinpath(output_dir,"ping.png"),UInt8.(cart))

##
# now let's process all radar variables in here

fullsweep = zeros(size(cart))
last_angle = 0.
sweep = 0;

for  i in 1:length(allRadarVariables)
    msg = fetchRadarMsg(allRadarVariables[i])

    delta_angle = msg.angle_increment
    # TODO: check upper bound; may have valid but coarse sectors
    if (msg.angle_increment<0 && msg.angle_increment > 0.2)
        # wrap around issues caused by driver: angle_increment
        println("Angle increment exceeds bounds (",msg.angle_increment,"); using default value")
        delta_angle = 0.17578125;
    end

    println(String(allRadarVariables[i].label)," - start:",msg.angle_start,", stop:",msg.angle_end, " step:",msg.angle_increment, " steps:",msg.n)
    (polar, cart) = assembleping(msg, 1.0)
    global fullsweep # need global due to julia's scope handling
    global last_angle
    global sweep
    fullsweep = fullsweep + cart

    # TODO: de-rotate angles by vehicle heading, so that test checks for full sweep in world frame
    if (msg.angle_end < msg.angle_start)# || msg.angle_start<last_angle)
        println("!")
        # now that we have a full sweep, we can try to match with the previous one!
        # TODO: register w/  previous sweep
        # save to disk
        fname = join(["full_",string(sweep),".png"])
        save(joinpath(output_dir,fname),UInt8.(clamp.(fullsweep,0,255)))

        # add full sweep to last variable
        ade,adb = addData!(fg, :radar, allRadarVariables[i].label, :RADARSWEEP, Vector{UInt8}(JSON2.write(fullsweep)), mimeType="application/json")

        # reset sweep
        fullsweep = zeros(size(cart))
        sweep = sweep+1
    end
    # fname = join(["polar_",String(allRadarVariables[i].label),".png"])
    # save(joinpath(output_dir,fname),UInt8.(polar))
    fname = join(["cart_",String(allRadarVariables[i].label),".png"])
    save(joinpath(output_dir,fname),UInt8.(cart))
    last_angle = msg.angle_end
end

##

# count(v -> :RADARSWEEP in listDataEntries(v), getVariables(fg))

# # may be too intense as it is fetching ALL data too?
# allSweepVariables = filter(v -> :RADARSWEEP in listDataEntries(v), getVariables(fg));

# allSweepVariables = sortDFG(allSweepVariables)

# getLabel.(allSweepVariables)

# de, db = getData(fg, allSweepVariables[1].label, :RADARSWEEP)

# JSON2.read(IOBuffer(db))

#


text:
#  Copyright 2016, Iain Dunning, Joey Huchette, Miles Lubin, and contributors
#  This Source Code Form is subject to the terms of the Mozilla Public
#  License, v. 2.0. If a copy of the MPL was not distributed with this
#  file, You can obtain one at http://mozilla.org/MPL/2.0/.
#############################################################################
# JuMP
# An algebraic modelling langauge for Julia
# See http://github.com/JuliaOpt/JuMP.jl
#############################################################################
# Hock-Schittkowski Nonlinear Test Suite  -  HS116
# This file is JuMP implementation of the model described in
#  W. Hock, K. Schittkowski, Test Examples for Nonlinear Programming
#  Codes, Lecture Notes in Economics and Mathematical Systems,
#  Springer, No, 187, 1981
# More information, including original model description, at
# http://www.ai7.uni-bayreuth.de/downloads.htm
#
# This problem has a linear objective with quadratic constraints.
#############################################################################

using JuMP
using Base.Test

let

N = 13
a = 0.002
b = 1.262626
c = 1.231059
d = 0.03475
e = 0.975
f = 0.00975

lower = [0.1, 0.1, 0.1, 0.0001, 0.1, 0.1, 0.1, 0.1, 500, 0.1, 1.0, 0.0001, 0.0001, 0.0, 0.0, 0.0]
upper = [1.0, 1.0, 1.0, 0.1, 0.9, 0.9, 1000, 1000, 1000, 500, 150, 150, 150, Inf, Inf, Inf]
start = [0.5  2 0.8  3 0.9  4 0.1  5 0.14  6 0.5  7 489  8 80  9 650 0.5  2 0.8  3 0.9  4 0.1  5 0.14  6 0.5  7 489  8 80  9 650]

m = Model()
@variable(m, lower[i] <= x[i=1:N] <= upper[i], start = start[i])
@NLobjective(m, Min, x[11] + x[12] + x[13])

@NLconstraint(m, x[3] - x[2] >= 0)
@NLconstraint(m, x[2] - x[1] >= 0)
@NLconstraint(m, 1 - a * x[7] + a * x[8] >= 0)
@NLconstraint(m, x[11] + x[12] + x[13] >= 50)
@NLconstraint(m, x[13] - b * x[10] + c * x[3] * x[10] >= 0)
@NLconstraint(m, x[5] - d * x[2] - e * x[2] * x[5] + f * x[2]^2 >= 0)
@NLconstraint(m, x[6] - d * x[3] - e * x[3] * x[6] + f * x[3]^2 >= 0)
@NLconstraint(m, x[4] - d * x[1] - e * x[1] * x[4] + f * x[1]^2 >= 0)
@NLconstraint(m, x[12] - b * x[9] + c * x[2] * x[9] >= 0)
@NLconstraint(m, x[11] - b * x[8] + c * x[1] * x[8] >= 0)
@NLconstraint(m, x[5] * x[7] - x[1] * x[8] - x[4] * x[7] + x[4] * x[8] >= 0)
@NLconstraint(m, 1 - a * (x[2] * x[9] + x[5] * x[8] - x[1] * x[8] - x[6] * x[9]) -
          x[5] - x[6] >= 0)
@NLconstraint(m, x[2] * x[9] - x[3] * x[10] - x[6] * x[9] - 500 * x[2] +
          500 * x[6] + x[2] * x[10] >= 0)
@NLconstraint(m, x[2] - 0.9 - a * (x[2] * x[10] - x[3] * x[10]) >= 0)
@NLconstraint(m, x[11] + x[12] + x[13] <= 250)

solve(m)

@test_approx_eq_eps getobjectivevalue(m) 97.588409 1e-3

end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

"""
Tools for working with the Markdown file format. Mainly for documentation.
"""
module Markdown

import Base: show, ==
import Core: @doc_str

include(joinpath("parse", "config.jl"))
include(joinpath("parse", "util.jl"))
include(joinpath("parse", "parse.jl"))

include(joinpath("Common", "Common.jl"))
include(joinpath("GitHub", "GitHub.jl"))
include(joinpath("IPython", "IPython.jl"))
include(joinpath("Julia", "Julia.jl"))

include(joinpath("render", "plain.jl"))
include(joinpath("render", "html.jl"))
include(joinpath("render", "latex.jl"))
include(joinpath("render", "rst.jl"))

include(joinpath("render", "terminal", "render.jl"))

export readme, license, @md_str, @doc_str

parse(markdown::AbstractString; flavor = julia) = parse(IOBuffer(markdown), flavor = flavor)
parse_file(file::AbstractString; flavor = julia) = parse(read(file, String), flavor = flavor)

readme(pkg::AbstractString; flavor = github) = parse_file(Pkg.dir(pkg, "README.md"), flavor = flavor)
readme(pkg::Module; flavor = github) = readme(string(pkg), flavor = flavor)

license(pkg::AbstractString; flavor = github) = parse_file(Pkg.dir(pkg, "LICENSE.md"), flavor = flavor)
license(pkg::Module; flavor = github) = license(string(pkg), flavor = flavor)

function mdexpr(s, flavor = :julia)
    md = parse(s, flavor = Symbol(flavor))
    esc(toexpr(md))
end

function docexpr(source::LineNumberNode, mod::Module, s, flavor = :julia)
    :($doc_str($(mdexpr(s, flavor)), $(QuoteNode(source)), $mod))
end

macro md_str(s, t...)
    mdexpr(s, t...)
end

function doc_str(md, source::LineNumberNode, mod::Module)
    md.meta[:path] = isa(source.file, Symbol) ? String(source.file) : ""
    md.meta[:module] = mod
    md
end
doc_str(md::AbstractString, source::LineNumberNode, mod::Module) = doc_str(parse(md), source, mod)

macro doc_str(s::AbstractString, t...)
    docexpr(__source__, __module__, s, t...)
end

function Base.display(d::Base.REPL.REPLDisplay, mds::Vector{MD})
    for md in mds
        display(d, md)
    end
end

end

max_stars_repo_name: xorJane/julia

max_issues_repo_name: JuliaLangEs/julieta

max_forks_repo_name: JuliaLangEs/julieta


text:
using Lale
using DataFrames
using CSV

strfeat= laleoperator("Project","lale")(columns=Dict("type"=>"string"))
numfeat = laleoperator("Project","lale")(columns=Dict("type"=>"number"))

OneHotEncoder = laleoperator("OneHotEncoder")
ConcatFeatures = laleoperator("ConcatFeatures","lale")
Tree = laleoperator("DecisionTreeClassifier")
KNN = laleoperator("KNeighborsClassifier")
RF =  laleoperator("RandomForestClassifier")
PCA = laleoperator("PCA")
Hyperopt = laleoperator("Hyperopt","lale")

df = CSV.read(joinpath(dirname(pathof(Lale)),"../demo/old/credit.csv"),DataFrame)
y = df[!,"class"] |> collect
X = df[:,2:end]

train_X,train_y,test_X,test_y = Lale.train_test_split(X, y,testprop=0.20)

prep_to_numbers =
 ((strfeat >> OneHotEncoder(handle_unknown = "ignore")) & numfeat)>> ConcatFeatures
planned_orig = prep_to_numbers >> ( Tree | KNN)
lopt = LalePipeOptimizer(planned_orig,max_evals = 10,cv = 3)
laletrained = fit(lopt,train_X,train_y)
lalepred  = Lale.transform(laletrained,test_X)
score(:accuracy,lalepred,test_y) |> println

pca_tree_planned = prep_to_numbers >> PCA >> RF
lopt = LalePipeOptimizer(pca_tree_planned,max_evals = 10,cv = 3)
laletrained = fit(lopt,train_X,train_y)
lalepred  = Lale.transform(laletrained,test_X)
score(:accuracy,lalepred,test_y) |> println


text:
using Pkg: Pkg
# import GitCommand
using HTTP: HTTP
using RegistryTools: RegistryTools
using Test: Test

function gather_stdlib_uuids()
    return Set{Base.UUID}(x for x in keys(RegistryTools.stdlibs()))
end

@static if Base.VERSION >= v"1.7.0-"
    const collect_registries = Pkg.Registry.reachable_registries
else
    const collect_registries = Pkg.Types.collect_registries
end

is_valid_url(str::AbstractString) = !isempty(HTTP.URI(str).scheme) && isvalid(HTTP.URI(str))

function _include_this_registry(
    registry_spec, registry_deps_names::Vector{<:AbstractString}
)
    all_fieldvalues = []
    for fieldname in [:name, :repo, :url]
        if hasproperty(registry_spec, fieldname)
            fieldvalue = getproperty(registry_spec, fieldname)
            if fieldvalue isa AbstractString
                push!(all_fieldvalues, fieldvalue)
            end
        end
    end
    for fieldvalue in all_fieldvalues
        for registry_deps_name in registry_deps_names
            if strip(fieldvalue) == strip(registry_deps_name)
                return true
            end
            if strip(fieldvalue) == strip("$(registry_deps_name).git")
                return true
            end
            if strip("$(fieldvalue).git") == strip(registry_deps_name)
                return true
            end
        end
    end
    return false
end

# For when you have a registry that has packages with dependencies obtained from
# another dependency registry. For example, packages registered at the BioJuliaRegistry
# that have General dependencies. BJW.
function load_registry_dep_uuids(registry_deps_names::Vector{<:AbstractString}=String[])
    return with_temp_depot() do
        # Get the registries!
        for repo_spec in registry_deps_names
            if is_valid_url(repo_spec)
                Pkg.Registry.add(Pkg.RegistrySpec(; url=repo_spec))
            else
                Pkg.Registry.add(repo_spec)
            end
        end
        # Now use the RegistrySpec's to find the Project.toml's. I know
        # .julia/registires/XYZ/ABC is the most likely place, but this way the
        # function never has to assume. BJW.
        extrauuids = Set{Base.UUID}()
        for spec in collect_registries()
            if _include_this_registry(spec, registry_deps_names)
                reg = Pkg.TOML.parsefile(joinpath(spec.path, "Registry.toml"))
                for x in keys(reg["packages"])
                    push!(extrauuids, Base.UUID(x))
                end
            end
        end
        return extrauuids
    end
end

#########################
# Testing of registries #
#########################

function load_package_data(
    ::Type{T}, path::String, versions::Vector{VersionNumber}
) where {T}
    compressed = Pkg.TOML.parsefile(path)
    compressed = convert(Dict{String,Dict{String,Union{String,Vector{String}}}}, compressed)
    uncompressed = Dict{VersionNumber,Dict{String,T}}()
    # Many of the entries are repeated so we keep a cache so there is no need to re-create
    # a bunch of identical objects
    cache = Dict{String,T}()
    vsorted = sort(versions)
    for (vers, data) in compressed
        vs = Pkg.Types.VersionRange(vers)
        first = length(vsorted) + 1
        # We find the first and last version that are in the range
        # and since the versions are sorted, all versions in between are sorted
        for i in eachindex(vsorted)
            v = vsorted[i]
            v in vs && (first = i; break)
        end
        last = 0
        for i in reverse(eachindex(vsorted))
            v = vsorted[i]
            v in vs && (last = i; break)
        end
        for i in first:last
            v = vsorted[i]
            uv = get!(() -> Dict{String,T}(), uncompressed, v)
            for (key, value) in data
                if haskey(uv, key)
                    error("Overlapping ranges for $(key) in $(repr(path)) for version $v.")
                else
                    Tvalue = if value isa String
                        Tvalue = get!(() -> T(value), cache, value)
                    else
                        Tvalue = T(value)
                    end
                    uv[key] = Tvalue
                end
            end
        end
    end
    return uncompressed
end

function load_deps(depsfile, versions)
    r =
        load_package_data(Base.UUID, depsfile, versions) isa
        Dict{VersionNumber,Dict{String,Base.UUID}}
    return r
end

function load_compat(compatfile, versions)
    r =
        load_package_data(Pkg.Types.VersionSpec, compatfile, versions) isa
        Dict{VersionNumber,Dict{String,Pkg.Types.VersionSpec}}
    return r
end

"""
    test(path)

Run various checks on the registry located at `path`.
Checks for example that all files are parsable and
understandable by Pkg and consistency between Registry.toml
and each Package.toml.

If your registry has packages that have dependencies that are registered in other
registries elsewhere, then you may provide the github urls for those registries
using the `registry_deps` parameter.
"""
function test(path=pwd(); registry_deps::Vector{<:AbstractString}=String[])
    Test.@testset "(Registry|Package|Versions|Deps|Compat).toml" begin
        cd(path) do
            reg = Pkg.TOML.parsefile("Registry.toml")
            reguuids = Set{Base.UUID}(Base.UUID(x) for x in keys(reg["packages"]))
            stdlibuuids = gather_stdlib_uuids()
            registry_dep_uuids = load_registry_dep_uuids(registry_deps)
            alluuids = reguuids ∪ stdlibuuids ∪ registry_dep_uuids

            # Test that each entry in Registry.toml has a corresponding Package.toml
            # at the expected path with the correct uuid and name
            Test.@testset "$(get(data, "name", uuid))" for (uuid, data) in reg["packages"]
                # Package.toml testing
                pkg = Pkg.TOML.parsefile(abspath(data["path"], "Package.toml"))
                Test.@test Base.UUID(uuid) == Base.UUID(pkg["uuid"])
                Test.@test data["name"] == pkg["name"]
                Test.@test Base.isidentifier(data["name"])
                Test.@test haskey(pkg, "repo")

                # Versions.toml testing
                vers = Pkg.TOML.parsefile(abspath(data["path"], "Versions.toml"))
                vnums = VersionNumber.(keys(vers))
                for (v, data) in vers
                    Test.@test VersionNumber(v) isa VersionNumber
                    Test.@test haskey(data, "git-tree-sha1")
                end

                # Deps.toml testing
                depsfile = abspath(data["path"], "Deps.toml")
                if isfile(depsfile)
                    deps = Pkg.TOML.parsefile(depsfile)
                    # Require all deps to exist in the General registry or be a stdlib
                    depuuids = Set{Base.UUID}(
                        Base.UUID(x) for (_, d) in deps for (_, x) in d
                    )
                    if !(depuuids ⊆ alluuids)
                        msg = "It is not the case that all dependencies exist in the General registry"
                        @error msg setdiff(depuuids, alluuids)
                    end
                    Test.@test depuuids ⊆ alluuids
                    # Test that the way Pkg loads this data works
                    Test.@test load_deps(depsfile, vnums)
                    # Make sure the content roundtrips through decompression/compression
                    compressed = RegistryTools.Compress.compress(
                        depsfile, RegistryTools.Compress.load(depsfile)
                    )
                    Test.@test compressed == deps
                end

                # Compat.toml testing
                compatfile = abspath(data["path"], "Compat.toml")
                if isfile(compatfile)
                    compat = Pkg.TOML.parsefile(compatfile)
                    # Test that all names with compat is a dependency
                    compatnames = Set{String}(x for (_, d) in compat for (x, _) in d)
                    if !(
                        isempty(compatnames) ||
                        (length(compatnames) == 1 && "julia" in compatnames)
                    )
                        depnames = Set{String}(
                            x for (_, d) in Pkg.TOML.parsefile(depsfile) for (x, _) in d
                        )
                        push!(depnames, "julia") # All packages has an implicit dependency on julia
                        if !(compatnames ⊆ depnames)
                            throw(
                                ErrorException("Assertion failed: compatnames ⊆ depnames")
                            )
                        end
                    end
                    # Test that the way Pkg loads this data works
                    Test.@test load_compat(compatfile, vnums)
                    # Make sure the content roundtrips through decompression/compression
                    compressed = RegistryTools.Compress.compress(
                        compatfile, RegistryTools.Compress.load(compatfile)
                    )
                    Test.@test compressed == compat
                end
            end
            # Make sure all paths are unique
            path_parts = [splitpath(data["path"]) for (_, data) in reg["packages"]]
            for i in 1:maximum(length, path_parts)
                i_parts = Set(
                    joinpath(x[1:i]...) for
                    x in path_parts if get(x, i, nothing) !== nothing
                )
                i_parts′ = Set(
                    joinpath(lowercase.(x[1:i])...) for
                    x in path_parts if get(x, i, nothing) !== nothing
                )
                Test.@test length(i_parts) == length(i_parts′)
            end
        end
    end
    return nothing
end


text:
export
    ccifrm,
    cgv2el,
    chbder,
    cidfrm,
    ckcls,
    ckcov!,
    ckcov,
    ckgp,
    ckgpav,
    cklpf,
    ckobj!,
    ckobj,
    ckopn,
    ckupf,
    ckw01,
    ckw02,
    ckw03,
    ckw05,
    clight,
    clpool,
    cmprss,
    cnmfrm,
    conics,
    convrt,
    cpos,
    cposr,
    cvpool,
    cyllat,
    cylrec,
    cylsph

"""
    ccifrm(frclss, clssid)

Return the frame name, frame ID, and center associated with a given frame class and class ID.

### Arguments ###

- `frclss`: Class of frame
- `clssid`: Class ID of frame

### Output ###

Returns `nothing` if no frame was found or

- `frcode`: ID code of the frame
- `frname`: Name of the frame
- `center`: ID code of the center of the frame

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/ccifrm_c.html)
"""
function ccifrm(frclss, clssid)
    lenout = 33
    frcode = Ref{SpiceInt}()
    frname = Array{SpiceChar}(undef, lenout)
    center = Ref{SpiceInt}()
    found = Ref{SpiceBoolean}()
    ccall((:ccifrm_c, libcspice), Cvoid,
          (SpiceInt, SpiceInt, SpiceInt, Ref{SpiceInt}, Ref{SpiceChar},
           Ref{SpiceInt}, Ref{SpiceBoolean}),
          frclss, clssid, lenout, frcode, frname, center, found)
    handleerror()
    Bool(found[]) || return nothing
    Int(frcode[]), chararray_to_string(frname), Int(center[])
end

"""
    cgv2el(center, vec1, vec2)

Form an ellipse from a center vector and two generating vectors.

### Arguments ###

- `center`: Center vector
- `vec1`: Generating vector
- `vec2`: Generating vector

### Output ###

Returns the ellipse defined by the input vectors.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/cgv2el_c.html)
"""
function cgv2el(center, vec1, vec2)
    @checkdims 3 center vec1 vec2
    ellipse = Ref{Ellipse}()
    ccall((:cgv2el_c, libcspice), Cvoid,
        (Ref{SpiceDouble}, Ref{SpiceDouble}, Ref{SpiceDouble}, Ref{Ellipse}),
        center, vec1, vec2, ellipse)
    ellipse[]
end

"""
    chbder(cp, x2s, x, nderiv)

Given the coefficients for the Chebyshev expansion of a polynomial, this returns the value of the
polynomial and its first `nderiv` derivatives evaluated at the input `x`.

### Arguments ###

- `cp`: Chebyshev polynomial coefficients
- `x2s`: Transformation parameters of polynomial
- `x`: Value for which the polynomial is to be evaluated
- `nderiv`: The number of derivatives to compute

### Output ###

Returns the derivatives of the polynomial.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/chbder_c.html)
"""
function chbder(cp, x2s, x, nderiv)
    @checkdims 2 x2s
    degp = length(cp) - 1
    partdp = Array{SpiceDouble}(undef, 3 * (nderiv + 1))
    dpdxs = Array{SpiceDouble}(undef, nderiv + 1)
    ccall((:chbder_c, libcspice), Cvoid,
          (Ref{SpiceDouble}, SpiceInt, Ref{SpiceDouble}, SpiceDouble, SpiceInt,
           Ref{SpiceDouble}, Ref{SpiceDouble}),
          cp, degp, x2s, x, nderiv, partdp, dpdxs)
    dpdxs
end

"""
    cidfrm(cent)

Retrieve frame ID code and name to associate with a frame center.

### Arguments ###

- `cent`: ID code for an object for which there is a preferred reference frame

### Output ###

Returns `nothing` if no frame was found or

- `frcode`: The ID code of the frame associated with `cent`
- `frname`: The name of the frame with ID `frcode`

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/cidfrm_c.html)
"""
function cidfrm(cent)
    lenout = 33
    frcode = Ref{SpiceInt}()
    frname = Array{SpiceChar}(undef, lenout)
    found = Ref{SpiceBoolean}()
    ccall((:cidfrm_c, libcspice), Cvoid,
          (SpiceInt, SpiceInt, Ref{SpiceInt}, Ref{SpiceChar}, Ref{SpiceBoolean}),
          cent, lenout, frcode, frname, found)
    Bool(found[]) || return nothing
    Int(frcode[]), chararray_to_string(frname)
end

"""
    ckcls(handle)

Close an open CK file.

### Arguments ###

- `handle`: Handle of the CK file to be closed

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/ckcls_c.html)
"""
function ckcls(handle)
    ccall((:ckcls_c, libcspice), Cvoid, (SpiceInt,), handle)
    handleerror()
end

"""
    ckcov(ck, idcode, needav, level, tol, timsys)

Find the coverage window for a specified object in a specified CK file.

### Arguments ###

- `ck`: Name of CK file
- `idcode`: ID code of object
- `needav`: Flag indicating whether angular velocity is needed
- `level`: Coverage level:  "SEGMENT" OR "INTERVAL"
- `tol`: Tolerance in ticks
- `timsys`: Time system used to represent coverage

### Output ###

Window giving coverage for `idcode`

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/ckcov_c.html)
"""
function ckcov(ck, idcode, needav, level, tol, timsys)
    cover = SpiceDoubleCell(20000)
    ckcov!(ck, idcode, needav, level, tol, timsys, cover)
end

"""
    ckcov!(ck, idcode, needav, level, tol, timsys, cover)

Find the coverage window for a specified object in a specified CK file.

### Arguments ###

- `ck`: Name of CK file
- `idcode`: ID code of object
- `needav`: Flag indicating whether angular velocity is needed
- `level`: Coverage level:  "SEGMENT" OR "INTERVAL"
- `tol`: Tolerance in ticks
- `timsys`: Time system used to represent coverage
- `cover`: Window giving coverage for `idcode`. Data already present in `cover`
    will be combined with coverage found for the object designated by `idcode`
    in the file `ck`.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/ckcov_c.html)
"""
function ckcov!(ck, idcode, needav, level, tol, timsys, cover)
    data = cover.data
    GC.@preserve data begin
        ccall((:ckcov_c, libcspice), Cvoid,
              (Cstring, SpiceInt, SpiceBoolean, Cstring, SpiceDouble, Cstring,
               Ref{Cell{SpiceDouble}}),
              ck, idcode, needav, level, tol, timsys, cover.cell)
    end
    handleerror()
    cover
end

"""
    ckgp(inst, sclkdp, tol, ref)

Get pointing (attitude) for a specified spacecraft clock time.

### Arguments ###

- `inst`: NAIF ID of instrument, spacecraft, or structure
- `sclkdp`: Encoded spacecraft clock time
- `tol`: Time tolerance
- `ref`: Reference frame

### Outputs ###

Returns `nothing` if the requested pointing is not available or

- `cmat`: C-matrix pointing data
- `clkout`: Output encoded spacecraft clock time

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/ckgp_c.html)
"""
function ckgp(inst, sclkdp, tol, ref)
    cmat = Array{SpiceDouble}(undef, 3, 3)
    clkout = Ref{SpiceDouble}()
    found = Ref{SpiceBoolean}()
    ccall((:ckgp_c, libcspice), Cvoid,
          (SpiceInt, SpiceDouble, SpiceDouble, Cstring,
           Ref{SpiceDouble}, Ref{SpiceDouble}, Ref{SpiceBoolean}),
          inst, sclkdp, tol, ref, cmat, clkout, found)
    handleerror()
    Bool(found[]) || return nothing
    permutedims(cmat), clkout[]
end

"""
    ckgpav(inst, sclkdp, tol, ref)

Get pointing (attitude) and angular velocity for a specified spacecraft clock time.

### Arguments ###

- `inst`: NAIF ID of instrument, spacecraft, or structure
- `sclkdp`: Encoded spacecraft clock time
- `tol`: Time tolerance
- `ref`: Reference frame

### Outputs ###

Returns `nothing` if the requested pointing is not available or

- `cmat`: C-matrix pointing data
- `av`: Angular velocity vector
- `clkout`: Output encoded spacecraft clock time

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/ckgpav_c.html)
"""
function ckgpav(inst, sclkdp, tol, ref)
    cmat = Array{SpiceDouble}(undef, 3, 3)
    av = Array{SpiceDouble}(undef, 3)
    clkout = Ref{SpiceDouble}()
    found = Ref{SpiceBoolean}()
    ccall((:ckgpav_c, libcspice), Cvoid,
        (SpiceInt, SpiceDouble, SpiceDouble, Cstring,
        Ref{SpiceDouble}, Ref{SpiceDouble}, Ref{SpiceDouble}, Ref{SpiceBoolean}),
        inst, sclkdp, tol, ref, cmat, av, clkout, found)
    handleerror()
    Bool(found[]) || return nothing
    permutedims(cmat), av, clkout[]
end

"""
    cklpf(filename)

Load a CK pointing file for use by the CK readers.  Return that file's handle, to be used by other CK
routines to refer to the file.

### Arguments ###

- `filename`: Name of the CK file to be loaded

### Output ###

Loaded file's handle

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/cklpf_c.html)
"""
function cklpf(filename)
    handle = Ref{SpiceInt}()
    ccall((:cklpf_c, libcspice), Cvoid, (Cstring, Ref{SpiceInt}), filename, handle)
    handleerror()
    handle[]
end

"""
    ckobj(ck)

Find the set of ID codes of all objects in a specified CK file.

### Arguments ##

- `ck`: Name of CK file

### Output ###

Set of ID codes of objects in CK file.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/ckobj_c.html)
"""
function ckobj(ck)
    ids = SpiceIntCell(1000)
    ckobj!(ck, ids)
end

"""
    ckobj!(ck, ids)

Find the set of ID codes of all objects in a specified CK file.

### Arguments ##

- `ck`: Name of CK file
- `ids`: Set of ID codes of objects in CK file. Data already present in
    `ids` will be combined with ID code set found for the file `ck`.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/ckobj_c.html)
"""
function ckobj!(ck, ids)
    data = ids.data
    GC.@preserve data begin
        ccall((:ckobj_c, libcspice), Cvoid,
              (Cstring, Ref{Cell{SpiceInt}}),
              ck, ids.cell)
    end
    handleerror()
    ids
end

"""
    ckopn(fname, ifname="CK_file", ncomch=0)

Open a new CK file, returning the handle of the opened file.

### Arguments ###

- `fname`: The name of the CK file to be opened
- `ifname`: The internal filename for the CK (default: "CK_file")
- `ncomch`: The number of characters to reserve for comments (default: 0)

### Output ###

- `handle`: The handle of the opened CK file

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/ckopn_c.html)
"""
function ckopn(fname, ifname="CK_file", ncomch=0)
    handle = Ref{SpiceInt}()
    ccall((:ckopn_c, libcspice), Cvoid,
          (Cstring, Cstring, SpiceInt, Ref{SpiceInt}),
          fname, ifname, ncomch, handle)
    handleerror()
    handle[]
end

"""
    ckupf(handle)

Unload a CK pointing file so that it will no longer be searched by the readers.

### Arguments ###

- `handle`: Handle of CK file to be unloaded

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/ckupf_c.html)
"""
function ckupf(handle)
    ccall((:ckupf_c, libcspice), Cvoid, (SpiceInt,), handle)
end

"""
    ckw01(handle, begtim, endtim, inst, ref, segid, sclkdp, quats, avvs=[zeros(3)])

Add a type 1 segment to a C-kernel.

### Arguments ###

- `handle`: Handle of an open CK file
- `begtim`: The beginning encoded SCLK of the segment
- `endtim`: The ending encoded SCLK of the segment
- `inst`: The NAIF instrument ID code
- `ref`: The reference frame of the segment
- `segid`: Segment identifier
- `sclkdp`: Encoded SCLK times
- `quats`: Quaternions representing instrument pointing
- `avvs`: Angular velocity vectors (optional)

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/ckw01_c.html)
"""
function ckw01(handle, begtim, endtim, inst, ref, segid, sclkdp, quats, avvs=[zeros(3)])
    quats_ = array_to_cmatrix(quats, n=4)
    avvs_ = array_to_cmatrix(avvs, n=3)
    nrec = length(sclkdp)
    avflag = length(avvs) == length(quats)
    ccall((:ckw01_c, libcspice), Cvoid,
          (SpiceInt, SpiceDouble, SpiceDouble, SpiceInt, Cstring, SpiceInt, Cstring, SpiceInt,
           Ref{SpiceDouble}, Ref{SpiceDouble}, Ref{SpiceDouble}),
          handle, begtim, endtim, inst, ref, avflag, segid, nrec, sclkdp, quats_, avvs_)
    handleerror()
end

"""
    ckw02(handle, begtim, endtim, inst, ref, segid, start, stop, quats, avvs, rates)

Write a type 2 segment to a C-kernel.

### Arguments ###

- `handle`: Handle of an open CK file
- `begtim`: The beginning encoded SCLK of the segment
- `endtim`: The ending encoded SCLK of the segment
- `inst`: The NAIF instrument ID code
- `ref`: The reference frame of the segment
- `segid`: Segment identifier
- `start`: Encoded SCLK interval start times
- `stop`: Encoded SCLK interval stop times
- `quats`: Quaternions representing instrument pointing
- `avvs`: Angular velocity vectors
- `rates`: Number of seconds per tick for each interval

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/ckw02_c.html)
"""
function ckw02(handle, begtim, endtim, inst, ref, segid, start, stop, quats, avvs, rates)
    nrec = length(quats)
    @checkdims nrec start stop avvs rates
    quats_ = array_to_cmatrix(quats, n=4)
    avvs_ = array_to_cmatrix(avvs, n=3)
    ccall((:ckw02_c, libcspice), Cvoid,
          (SpiceInt, SpiceDouble, SpiceDouble, SpiceInt, Cstring, Cstring, SpiceInt, Ref{SpiceDouble},
           Ref{SpiceDouble}, Ref{SpiceDouble}, Ref{SpiceDouble}, Ref{SpiceDouble}),
          handle, begtim, endtim, inst, ref, segid, nrec - 1, start, stop, quats_, avvs_, rates)
    handleerror()
end

"""
    ckw03(handle, begtim, endtim, inst, ref, segid, sclkdp, quats, starts, avvs=[zeros(3)])

Add a type 3 segment to a C-kernel.

### Arguments ###

- `handle`: Handle of an open CK file
- `begtim`: The beginning encoded SCLK of the segment
- `endtim`: The ending encoded SCLK of the segment
- `inst`: The NAIF instrument ID code
- `ref`: The reference frame of the segment
- `segid`: Segment identifier
- `sclkdp`: Encoded SCLK times
- `quats`: Quaternions representing instrument pointing
- `starts`: Encoded SCLK interval start times
- `avvs`: Angular velocity vectors (optional)

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/ckw03_c.html)
"""
function ckw03(handle, begtim, endtim, inst, ref, segid, sclkdp, quats, starts, avvs=[zeros(3)])
    nrec = length(quats)
    avflag = length(avvs) == nrec
    nints = length(starts)
    quats_ = array_to_cmatrix(quats, n=4)
    avvs_ = array_to_cmatrix(avvs, n=3)
    @checkdims nrec sclkdp
    ccall((:ckw03_c, libcspice), Cvoid,
          (SpiceInt, SpiceDouble, SpiceDouble, SpiceInt, Cstring, SpiceBoolean, Cstring,
           SpiceInt, Ref{SpiceDouble}, Ref{SpiceDouble}, Ref{SpiceDouble}, SpiceInt,
           Ref{SpiceDouble}),
          handle, begtim, endtim, inst, ref, avflag, segid, nrec - 1,
          sclkdp, quats_, avvs_, nints, starts)
    handleerror()
end

"""
    ckw05(handle, subtyp, degree, begtim, endtim, inst, ref, avflag, segid, sclkdp, packts,
          rate, nints, starts)

Write a type 5 segment to a CK file.

### Arguments ###

- `handle`: Handle of an open CK file
- `subtyp`: CK type 5 subtype code
- `degree`: Degree of interpolating polynomials
- `begtim`: The beginning encoded SCLK of the segment
- `endtim`: The ending encoded SCLK of the segment
- `inst`: The NAIF instrument ID code
- `ref`: The reference frame of the segment
- `avflag`: True if the segment will contain angular velocity
- `segid`: Segment identifier
- `sclkdp`: Encoded SCLK times
- `packts`: Array of packets
- `rate`: Nominal SCLK rate in seconds per tick
- `nints`: Number of intervals
- `starts`: Encoded SCLK interval start times

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/ckw05_c.html)
"""
function ckw05(handle, subtyp, degree, begtim, endtim, inst, ref, avflag, segid, sclkdp, packts,
               rate, nints, starts)
    n = length(packts)
    @checkdims n sclkdp
    if subtyp == 0
        packts_ = array_to_cmatrix(packts, n=8)
    elseif subtyp == 1
        packts_ = array_to_cmatrix(packts, n=4)
    elseif subtyp == 2
        packts_ = array_to_cmatrix(packts, n=14)
    elseif subtyp == 3
        packts_ = array_to_cmatrix(packts, n=7)
    end
    ccall((:ckw05_c, libcspice), Cvoid,
          (SpiceInt, SpiceInt, SpiceInt, SpiceDouble, SpiceDouble, SpiceInt, Cstring, SpiceBoolean,
           Cstring, SpiceInt, Ref{SpiceDouble}, Ref{SpiceDouble},
           SpiceDouble, SpiceInt, Ref{SpiceDouble}),
          handle, subtyp, degree, begtim, endtim, inst, ref, avflag, segid, n, sclkdp, packts_, rate,
          nints, starts)
    handleerror()
end

@deprecate cleard empty!

"""
    cleard(array)

!!! warning "Deprecated"
    Use `empty!(array)` instead.
"""
cleard

"""
    clight()

Returns the speed of light in vacuo (km/sec).

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/clight_c.html)
"""
function clight()
    ccall((:clight_c, libcspice), Cdouble, ())
end

"""
    clpool()

Remove all variables from the kernel pool. Watches on kernel variables are retained.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/clpool_c.html)
"""
function clpool()
    ccall((:clpool_c, libcspice), Cvoid, ())
end

"""
    cmprss(delim, n, input)

Compress a character string by removing occurrences of more than `n` consecutive occurrences
of a specified character.

### Arguments ###

- `delim`: Delimiter to be compressed
- `n`: Maximum consecutive occurrences of delim
- `input`: Input string

### Output ###

Returns the compressed string.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/cmprss_c.html)
"""
function cmprss(delim, n, input)
    lenout = length(input)
    output = Array{SpiceChar}(undef, lenout)
    ccall((:cmprss_c, libcspice), Cvoid,
          (Cchar, SpiceInt, Cstring, SpiceInt, Ref{SpiceChar}),
          first(delim), n, input, lenout, output)
    handleerror()
    chararray_to_string(output)
end

"""
    cnmfrm(cname)

Retrieve frame ID code and name to associate with an object.

### Arguments ###

- `cname`: Name of the object to find a frame for

### Output ###

Returns a tuple of the ID code and the name of the frame associated with `cname` or `nothing` if no
frame is found.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/cnmfrm_c.html)
"""
function cnmfrm(cname, lenout=256)
    frcode = Ref{SpiceInt}()
    frname = Array{SpiceChar}(undef, lenout)
    found = Ref{SpiceBoolean}()
    ccall((:cnmfrm_c, libcspice), Cvoid,
          (Cstring, SpiceInt, Ref{SpiceInt}, Ref{SpiceChar}, Ref{SpiceBoolean}),
          cname, lenout, frcode, frname, found)
    Bool(found[]) || return nothing
    frcode[], chararray_to_string(frname)
end

"""
    conics(elts, et)

Determine the state (position, velocity) of an orbiting body from a set of elliptic, hyperbolic,
or parabolic orbital elements.

### Arguments ###

- `elts`: Conic elements
- `et`: Input time

### Output ###

Returns the state of orbiting body at `et`.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/conics_c.html)
"""
function conics(elts, et)
    @checkdims 8 elts
    state = Array{Float64}(undef, 6)
    ccall((:conics_c, libcspice), Cvoid, (Ref{SpiceDouble}, SpiceDouble, Ref{SpiceDouble}), elts, et, state)
    handleerror()
    state
end

"""
    convrt(x, in, out)

Take a measurement `x`, the units associated with `x`, and units to which `x` should be converted;
return `y` - the value of the measurement in the output units.

### Arguments ###

- `x`: Number representing a measurement in some units
- `in`: The units in which x is measured
- `out`: Desired units for the measurement

### Output ###

Returns the measurement in the desired units.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/convrt_c.html)
"""
function convrt(x, in, out)
    y = Ref{SpiceDouble}()
    ccall((:convrt_c, libcspice), Cvoid, (SpiceDouble, Cstring, Cstring, Ref{SpiceDouble}),
          x, in, out, y)
    handleerror()
    y[]
end

"""
    cpos(str, chars, start)

Find the first occurrence in a string of a character belonging to a collection of characters,
starting at a specified location, searching forward.

### Arguments ###

- `str`: Any character string
- `chars`: A collection of characters
- `start`: Position to begin looking for one of chars

### Output ###

Returns the index of the first character of `str` that is one of the characters in string
`chars`. Returns -1 if none of the characters was found.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/cpos_c.html)
"""
function cpos(str, chars, start)
    idx = Int(ccall((:cpos_c, libcspice), SpiceInt,
                    (Cstring, Cstring, SpiceInt),
                    str, chars, start - 1))
    handleerror()
    idx != -1 ? idx + 1 : idx
end

"""
    cposr(str, chars, start)

Find the first occurrence in a string of a character belonging to a collection of characters,
starting at a specified location, searching in reverse.

### Arguments ###

- `str`: Any character string
- `chars`: A collection of characters
- `start`: Position to begin looking for one of chars

### Output ###

Returns the index of the last character of `str` that is one of the characters in string
`chars`. Returns -1 if none of the characters was found.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/cposr_c.html)
"""
function cposr(str, chars, start)
    idx = Int(ccall((:cposr_c, libcspice), SpiceInt,
                    (Cstring, Cstring, SpiceInt),
                    str, chars, start - 1))
    handleerror()
    idx != -1 ? idx + 1 : idx
end

"""
    cvpool(agent)

Indicate whether or not any watched kernel variables that have a specified agent on their
notification list have been updated.

### Arguments ###

- `agent`: Name of the agent to check for notices

### Output ###

Returns `true` if variables for `agent` have been updated.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/cvpool_c.html)
"""
function cvpool(agent)
    update = Ref{SpiceBoolean}()
    ccall((:cvpool_c, libcspice), Cvoid, (Cstring, Ref{SpiceBoolean}), agent, update)
    Bool(update[])
end

"""
    cyllat(r, lonc, z)

Convert from cylindrical to latitudinal coordinates.

### Arguments ###

- `r`: Distance of point from z axis
- `lonc`: Cylindrical angle of point from XZ plane (radians)
- `z`: Height of point above XY plane

### Output ###

- `radius`: Radius
- `lon`: Longitude (radians)
- `lat`: Latitude (radians)

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/cyllat_c.html)
"""
function cyllat(r, lonc, z)
    radius = Ref{SpiceDouble}()
    lon = Ref{SpiceDouble}()
    lat = Ref{SpiceDouble}()
    ccall((:cyllat_c, libcspice), Cvoid,
          (SpiceDouble, SpiceDouble, SpiceDouble, Ref{SpiceDouble}, Ref{SpiceDouble}, Ref{SpiceDouble}),
          r, lonc, z, radius, lon, lat)
    radius[], lon[], lat[]
end

"""
    cylrec(r, lon, z)

Convert from cylindrical to rectangular coordinates.

### Arguments ###

- `r`: Distance of the point of interest from z axis
- `lon`: Cylindrical angle (in radians) of the point of interest from XZ plane
- `z`: Height of the point above XY plane

### Output ###

Returns rectangular coordinates of the point of interest.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/cylrec_c.html)
"""
function cylrec(r, lon, z)
    rectan = Array{SpiceDouble}(undef, 3)
    ccall((:cylrec_c, libcspice), Cvoid, (SpiceDouble, SpiceDouble, SpiceDouble, Ref{SpiceDouble}),
          r, lon, z, rectan)
    rectan
end

"""
    cylsph(r, lonc, z)

Convert from cylindrical to spherical coordinates.

### Arguments ###

- `r`: Distance of point from z axis
- `lonc`: Angle (radians) of point from XZ plane
- `z`: Height of point above XY plane

### Output ###

- `radius`: Distance of the point from the origin
- `colat`: Polar angle (co-latitude in radians)
- `lon`: Azimuthal angle (longitude)

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/cylsph_c.html)
"""
function cylsph(r, lonc, z)
    radius = Ref{SpiceDouble}()
    colat = Ref{SpiceDouble}()
    lon = Ref{SpiceDouble}()
    ccall((:cylsph_c, libcspice), Cvoid,
          (SpiceDouble, SpiceDouble, SpiceDouble, Ref{SpiceDouble}, Ref{SpiceDouble}, Ref{SpiceDouble}),
          r, lonc, z, radius, colat, lon)
    radius[], colat[], lon[]
end


text:
#simple multiplot
using Plots
using Random
y = randn(100,4)
p1 = plot(x, y) # Make a line plot
p2 = scatter(x, y) # Make a scatter plot
p3 = scatter(x, y[:,1], ylabel = "y1", lw = 3, title = "Column 1 scatter")
p4 = histogram(x, y) # Four histograms each with 10 points? Why not!
f = plot(p1, p2, p3, p4, layout = (4, 1), legend = false)
display(f)


text:
##################################################################################
# This file is part of ModelBaseEcon.jl
# BSD 3-Clause License
# Copyright (c) 2020-2022, Bank of Canada
# All rights reserved.
##################################################################################

"""
struct EqnNotReadyError <: ModelErrorBase

Concrete error type used to indicate that a given equation has not been prepared
for use in the model yet.
"""
struct EqnNotReadyError <: ModelErrorBase end
msg(::EqnNotReadyError) = "Equation not ready to use."
hint(::EqnNotReadyError) = "Call `@initialize model` or `add_equation!()` first."

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

# Equation expressions typed by the user are of course valid equations, however
# during processing we use recursive algorithms, with the bottom of the recursion
# being a Number or a Symbol.  So we need a type that allows these.
const ExtExpr = Union{Expr,Symbol,Number}

"Placeholder evaluation function to use in Equation costruction while it is being created"
@inline eqnnotready(x...) = throw(EqnNotReadyError())

"""
    mutable struct EqnFlags ... end

Holds information about the equation. Flags can be specified in the model
definition by annotating the equation with `@<flag>` (insert the flag you want
to raise in place of `<flag>`).

Supported flags:
 * `@log lhs = rhs` instructs the model parser to make the residual
   `loc(lhs / rhs)`. Normally the residual is computed as `lhs - rhs`.
 * `@lin lhs = rhs` marks the equation for selective linearization.

"""
mutable struct EqnFlags
    lin::Bool
    log::Bool
    EqnFlags() = new(false, false)
    EqnFlags(lin, log) = new(lin, log)
end

Base.hash(f::EqnFlags, h::UInt) = hash(((f.:($flag) for flag in fieldnames(EqnFlags))...,), h)
Base.:(==)(f1::EqnFlags, f2::EqnFlags) = all(f1.:($flag) == f2.:($flag) for flag in fieldnames(EqnFlags))

export Equation
"""
    struct Equation <: AbstractEquation

Data structure representing a single equation in our state space model.

### Methods (for users)
  * TODO

### Implementation (for developers)
During the phase of definition of the Model, this type simply stores the expression
entered by the user. During @initialize(), the true data structure is constructed.
We need this, because the construction of the equation requires information from
the Model object, which may not be available at the time the equation expression
is first read.
"""
struct Equation <: AbstractEquation
    doc::String
    flags::EqnFlags
    "The original expression entered by the user"
    expr::ExtExpr      # original expression
    """
    The residual expression computed from [`expr`](@ref).
    It is used in the evaluation functions.
    Leads and lags of variables and shocks mentioned in [`expr`](@ref)
    are recorded in [`vinds`](@ref) and translated to symbols in [`vsyms`](@ref).
    """
    resid::Expr     # residual expression (with variables renamed and parameters substituted)
    """
    Indexes of variables/shocks mentioned in [`expr`](@ref). Each index is in the form
    `(t, v)`, where `t` is the lag (if negative), 0, or the lead (if positive) and `v`
    is the index of the variable/shock.
    """
    vinds::Vector{Tuple{Int64,Int64}}   # (t, v) indices of the relevant variables
    "Symbols used in `resid` replacing mentions of variables and shocks"
    vsyms::Vector{Symbol}  # symbols representing the relevant variables in resid
    "maximum lag mentioned in `expr`"
    maxlag::Int64   # for this equation
    "maximum lead mentioned in `expr`"
    maxlead::Int64  # for this equation
    "A callable (function) evaluating the residual. Argument is a vector of Float64 same lenght as `vinds`"
    eval_resid::Function  # function evaluating the residual
    "A callable (function) evaluating the (residual, gradient) pair. Argument is a vector of Float64 same lenght as `vinds`"
    eval_RJ::Function     # Function evaluating the residual and its gradient
end

#
# dummy constructor - just stores the expresstion without any processing
Equation(expr::ExtExpr) = Equation("", EqnFlags(), expr, Expr(:block), [], [], 0, 0, eqnnotready, eqnnotready)

# constructor that computes maxlag and maxlead on the fly
function Equation(doc, flags, expr, resid, vinds, vsyms, eval_resid, eval_RJ)
    # compute `maxlag` and `maxlead`
    maxlag, maxlead = (isempty(vinds) ? (0, 0) : extrema(v[1] for v in vinds) .* (-1, 1))
    # call the default constructor
    return Equation(doc, flags, expr, resid, vinds, vsyms, maxlag, maxlead, eval_resid, eval_RJ)
end

# Allows us to pass a Number of a Symbol or a raw Expr to calls where Equation is expected.
Base.convert(::Type{Equation}, e::ExtExpr) = Equation(e)


text:
# This handles loading all models that are nonsequential (Model type)

# load_layers(load_structure(model)["layers"])
# Create a graph from the above array.
#vcall_cat(a...) = vcall(:cat, 3, )

function graphify(a::Array{Any, 1}, structure_file, weight_file, ip)
    global weight = weights(weight_file)
    res = Dict{Any, Any}()
    for ele in a
        if ele.layer_type == :InputLayer
            res[ele.fields["name"]] = vcall(:identity, ip)
            #println(ele.fields["name"], "-->",res[ele.fields["name"]])
        elseif ele.layer_type == :Add
            inputs = ele.input_nodes[1]
            res[ele.fields["name"]] = vcall(ops[:Add](ele), res[inputs[1][1]], res[inputs[2][1]])
        elseif ele.layer_type == :Concatenate
            inputs = ele.input_nodes[1]
            ips = Array{Any, 1}()
            for ip in inputs
                push!(ips, ip[1])
            end
            if length(ips) == 4
                res[ele.fields["name"]] = vcall(:cat, 3, res[ips[1]],res[ips[2]],res[ips[3]],res[ips[4]])  #ToDo: Find an efficient way
            elseif length(ips) == 3                                                                         # for this.
                res[ele.fields["name"]] = vcall(:cat, 3, res[ips[1]],res[ips[2]],res[ips[3]])
            elseif length(ips) == 2
                res[ele.fields["name"]] = vcall(:cat, 3, res[ips[1]],res[ips[2]])
            end
        elseif ele.layer_type == :Dense
            op_dense = ops[:Dense](ele)[1]
            op_activation = ops[:Dense](ele)[2]
            inputs = ele.input_nodes[1][1][1]
            res[ele.fields["name"]] = vcall(Keras.ops[op_activation](ele), vcall(op_dense, res[inputs]))
        elseif ele.layer_type == :GlobalAveragePooling2D
            inputs = ele.input_nodes[1][1][1]
            res[ele.fields["name"]] = vcall(x -> reshape(x, size(x)[3], size(x)[4]),
                                        vcall(ops[:GlobalAveragePooling2D](ele), res[inputs]))
        else
            inputs = ele.input_nodes[1][1][1]
            res[ele.fields["name"]] = vcall(ops[ele.layer_type](ele), res[inputs])
            #println(ele.fields["name"], "-->",res[ele.fields["name"]])
        end
    end
    return res[get_outputlayer(structure_file)] |> syntax |> eval
end

function graphify_dummy(a::Array{Any, 1}, structure_file, weight_file, ip, num)
    global weight = weights(weight_file)
    res = Dict{Any, Any}()
    for ele in a[1:num]
        if ele.layer_type == :InputLayer
            res[ele.fields["name"]] = vcall(:identity, ip)
            #println(ele.fields["name"], "-->",res[ele.fields["name"]])
        elseif ele.layer_type == :Add
            inputs = ele.input_nodes[1]
            res[ele.fields["name"]] = vcall(ops[:Add](ele), res[inputs[1][1]], res[inputs[2][1]])
        else
            inputs = ele.input_nodes[1][1][1]
            res[ele.fields["name"]] = vcall(ops[ele.layer_type](ele), res[inputs])
            #println(ele.fields["name"], "-->",res[ele.fields["name"]])
        end
    end
    return res[a[num].fields["name"]] |> syntax |> eval
end

function get_outputlayer(structure_file)
    res = JSON.parse(String(read(open(structure_file, "r"))))
    return res["config"]["output_layers"][1][1]
end

function get_op(structure_file,weight_file,  ll)
    dic = graphify(ll,weight_file)
    op = get_outputlayer(structure_file)
    #op2 = vcall(:Chain, dic[op])
    return dic[op] |> syntax |> eval
end

function load_nonsequential_model(structure_file, weight_file)
    global weight = weights(weight_file)
    ll = load_layers(load_structure(structure_file)["layers"])
    return graphify(ll, structure_file, weight_file)
end


text:
struct HEPMASS <: Name end

# mandatory methods
format(::Type{HEPMASS}) = TabularData
task(::Type{HEPMASS}) = TwoClass
source(::Type{HEPMASS}) = "https://archive.ics.uci.edu/ml/datasets/HEPMASS"

function downloadlink(::Type{HEPMASS})
    return "https://archive.ics.uci.edu/ml/machine-learning-databases/00347/" .* [
        "all_train.csv.gz",
        "all_test.csv.gz"
    ]
end

function preprocess(N::Type{HEPMASS})
    kwargs = (istarget = 1, header = true, gzip = true)
    return [
        path -> preprocess_csv(N, path, :train; kwargs...),
        path -> preprocess_csv(N, path, :test; kwargs...),
    ]
end

nattributes(::Type{HEPMASS}) = (28, )
ninstances(::Type{HEPMASS}) = (7000000, 0, 3500000)
positivelabel(::Type{HEPMASS}) = 1

# optional methods
function checksum(::Type{HEPMASS})
    return "becdb27c603c40ba3a746736a224535585b8c19f8d02db7da28d0fc78726c52b"
end

name(::Type{HEPMASS}) = "HEPMASS"
author(::Type{HEPMASS}) = ["Daniel Whiteson"]

function licence(::Type{HEPMASS})
    return """
    This data set is licensed under a Creative Commons Attribution 4.0 International (CC BY 4.0) license. This allows for the sharing and adaptation of the datasets for any purpose, provided that the appropriate credit is given.

    URL: https://creativecommons.org/licenses/by/4.0/legalcode
    """
end

function citation(::Type{HEPMASS})
    return """
    @misc{misc_hepmass_347,
        title        = {{HEPMASS}},
        year         = {2016},
        howpublished = {UCI Machine Learning Repository}
    }
    """
end

max_stars_repo_name: VaclavMacha/DatasetProvider

max_issues_repo_name: VaclavMacha/DatasetProvider

max_forks_repo_name: VaclavMacha/DatasetProvider