text:
url = "https://raw.githubusercontent.com/pyviz/colorcet/master/assets/CETperceptual_csv_0_1/"

filelist = [
"cyclic_grey_15-85_c0_n256.csv",
"cyclic_grey_15-85_c0_n256_s25.csv",
"cyclic_mrybm_35-75_c68_n256.csv",
"cyclic_mrybm_35-75_c68_n256_s25.csv",
"cyclic_mygbm_30-95_c78_n256.csv",
"cyclic_mygbm_30-95_c78_n256_s25.csv",
"cyclic_wrwbw_40-90_c42_n256.csv",
"cyclic_wrwbw_40-90_c42_n256_s25.csv",
"diverging-isoluminant_cjm_75_c23_n256.csv",
"diverging-isoluminant_cjm_75_c24_n256.csv",
"diverging-isoluminant_cjo_70_c25_n256.csv",
"diverging-linear_bjr_30-55_c53_n256.csv",
"diverging-linear_bjy_30-90_c45_n256.csv",
"diverging-rainbow_bgymr_45-85_c67_n256.csv",
"diverging_bkr_55-10_c35_n256.csv",
"diverging_bky_60-10_c30_n256.csv",
"diverging_bwr_40-95_c42_n256.csv",
"diverging_bwr_55-98_c37_n256.csv",
"diverging_cwm_80-100_c22_n256.csv",
"diverging_gkr_60-10_c40_n256.csv",
"diverging_gwr_55-95_c38_n256.csv",
"diverging_gwv_55-95_c39_n256.csv",
"isoluminant_cgo_70_c39_n256.csv",
"isoluminant_cgo_80_c38_n256.csv",
"isoluminant_cm_70_c39_n256.csv",
"linear_bgy_10-95_c74_n256.csv",
"linear_bgyw_15-100_c67_n256.csv",
"linear_bgyw_15-100_c68_n256.csv",
"linear_blue_5-95_c73_n256.csv",
"linear_blue_95-50_c20_n256.csv",
"linear_bmw_5-95_c86_n256.csv",
"linear_bmw_5-95_c89_n256.csv",
"linear_bmy_10-95_c71_n256.csv",
"linear_bmy_10-95_c78_n256.csv",
"linear_gow_60-85_c27_n256.csv",
"linear_gow_65-90_c35_n256.csv",
"linear_green_5-95_c69_n256.csv",
"linear_grey_0-100_c0_n256.csv",
"linear_grey_10-95_c0_n256.csv",
"linear_kry_5-95_c72_n256.csv",
"linear_kry_5-98_c75_n256.csv",
"linear_kryw_0-100_c71_n256.csv",
"linear_kryw_5-100_c64_n256.csv",
"linear_kryw_5-100_c67_n256.csv",
"linear_ternary-blue_0-44_c57_n256.csv",
"linear_ternary-green_0-46_c42_n256.csv",
"linear_ternary-red_0-50_c52_n256.csv",
"rainbow_bgyr_35-85_c72_n256.csv",
"rainbow_bgyr_35-85_c73_n256.csv",
"rainbow_bgyrm_35-85_c69_n256.csv",
"rainbow_bgyrm_35-85_c71_n256.csv"]

function getdatafiles()
	filedata = []
	for fname in filelist
		tempfile = download(url * fname)
		fcontents = open(tempfile) do f
			push!(filedata, (fname, read(f, String)))
		end
	end
	return filedata
end

function build_colorschemes(io, colorcetdata)
	for i in 1:length(colorcetdata)
		schemename = colorcetdata[i][1]
		schemename = replace(schemename, "-" => "_")
		schemename = replace(schemename, ".csv" => "")
		rawcolorvalues = colorcetdata[i][2]
		println(io, "loadcolorscheme(:$(schemename), [")
		for l in split(rawcolorvalues, "\n")
		    !isempty(l) && println(io, "\tRGB($l),")
		end
		println(io, "], \"colorcet\", \"\")\n")
	end
end

# colorcetdata = getdatafiles()

build_colorschemes(stdin, colorcetdata)

open("/tmp/colorcetdata.jl", "w") do f
	build_colorschemes(f, colorcetdata)
end


text:
@testset "Abstract Potential Default Implmentation Unit Tests" begin
    atoms = [
        :Ar => (@SVector [0.0, 0.0, 0.0])u"bohr",
        :Ar => (@SVector [0.5, 0.0, 0.0])u"bohr",
        :H => (@SVector [0.0, 0.5, 0.0])u"bohr",
        :H => (@SVector [0.5, 0.5, 0.0])u"bohr"
    ]
    box = [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]]u"bohr"
    system = periodic_system(atoms, box)

    p = MockAbstractPotential(0.0)

    @test get_rcutoff(p) == Inf
    @test ismissing(get_species(p))

    @test potential_energy(system, p) == 0.0u"hartree"
    @test force(system, p) == fill((@SVector[1.0, 2.0, 3.0])u"hartree/bohr", 4)
    @test virial(system, p) == -6u"hartree"
end

ext: jl


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

module Artifacts

import Base: get, SHA1
using Base.BinaryPlatforms, Base.TOML

export artifact_exists, artifact_path, artifact_meta, artifact_hash,
       select_downloadable_artifacts, find_artifacts_toml, @artifact_str

"""
    parse_toml(path::String)

Uses Base.TOML to parse a TOML file. Do not mutate the returned dictionary.
"""
function parse_toml(path::String)
    # Uses the caching mechanics for toml files in Base
    Base.parsed_toml(path)
end

# keep in sync with Base.project_names and Base.manifest_names
const artifact_names = ("JuliaArtifacts.toml", "Artifacts.toml")

const ARTIFACTS_DIR_OVERRIDE = Ref{Union{String,Nothing}}(nothing)
"""
    with_artifacts_directory(f::Function, artifacts_dir::String)

Helper function to allow temporarily changing the artifact installation and search
directory.  When this is set, no other directory will be searched for artifacts, and new
artifacts will be installed within this directory.  Similarly, removing an artifact will
only effect the given artifact directory.  To layer artifact installation locations, use
the typical Julia depot path mechanism.
"""
function with_artifacts_directory(f::Function, artifacts_dir::String)
    try
        ARTIFACTS_DIR_OVERRIDE[] = artifacts_dir
        f()
    finally
        ARTIFACTS_DIR_OVERRIDE[] = nothing
    end
end

"""
    artifacts_dirs(args...)

Return a list of paths joined into all possible artifacts directories, as dictated by the
current set of depot paths and the current artifact directory override via the method
`with_artifacts_dir()`.
"""
function artifacts_dirs(args...)
    if ARTIFACTS_DIR_OVERRIDE[] === nothing
        return String[abspath(depot, "artifacts", args...) for depot in Base.DEPOT_PATH]
    else
        # If we've been given an override, use _only_ that directory.
        return String[abspath(ARTIFACTS_DIR_OVERRIDE[], args...)]
    end
end

"""
    ARTIFACT_OVERRIDES

Artifact locations can be overridden by writing `Override.toml` files within the artifact
directories of Pkg depots.  For example, in the default depot `~/.julia`, one may create
a `~/.julia/artifacts/Override.toml` file with the following contents:

    78f35e74ff113f02274ce60dab6e92b4546ef806 = "/path/to/replacement"
    c76f8cda85f83a06d17de6c57aabf9e294eb2537 = "fb886e813a4aed4147d5979fcdf27457d20aa35d"

    [d57dbccd-ca19-4d82-b9b8-9d660942965b]
    c_simple = "/path/to/c_simple_dir"
    libfoo = "fb886e813a4aed4147d5979fcdf27457d20aa35d""

This file defines four overrides; two which override specific artifacts identified
through their content hashes, two which override artifacts based on their bound names
within a particular package's UUID.  In both cases, there are two different targets of
the override: overriding to an on-disk location through an absolute path, and
overriding to another artifact by its content-hash.
"""
const ARTIFACT_OVERRIDES = Ref{Union{Dict{Symbol,Any},Nothing}}(nothing)
function load_overrides(;force::Bool = false)::Dict{Symbol, Any}
    if ARTIFACT_OVERRIDES[] !== nothing && !force
        return ARTIFACT_OVERRIDES[]
    end

    # We organize our artifact location overrides into two camps:
    #  - overrides per UUID with artifact names mapped to a new location
    #  - overrides per hash, mapped to a new location.
    #
    # Overrides per UUID/bound name are intercepted upon Artifacts.toml load, and new
    # entries within the "hash" overrides are generated on-the-fly.  Thus, all redirects
    # mechanisticly happen through the "hash" overrides.
    overrides = Dict{Symbol,Any}(
        # Overrides by UUID
        :UUID => Dict{Base.UUID,Dict{String,Union{String,SHA1}}}(),

        # Overrides by hash
        :hash => Dict{SHA1,Union{String,SHA1}}(),
    )

    for override_file in reverse(artifacts_dirs("Overrides.toml"))
        !isfile(override_file) && continue

        # Load the toml file
        depot_override_dict = parse_toml(override_file)

        function parse_mapping(mapping::String, name::String)
            if !isabspath(mapping) && !isempty(mapping)
                mapping = tryparse(Base.SHA1, mapping)
                if mapping === nothing
                    @error("Invalid override in '$(override_file)': entry '$(name)' must map to an absolute path or SHA1 hash!")
                end
            end
            return mapping
        end
        function parse_mapping(mapping::Dict, name::String)
            return Dict(k => parse_mapping(v, name) for (k, v) in mapping)
        end
        # Fallthrough for invalid Overrides.toml files
        parse_mapping(mapping, name::String) = nothing

        for (k, mapping) in depot_override_dict
            # First, parse the mapping. Is it an absolute path, a valid SHA1-hash, or neither?
            mapping = parse_mapping(mapping, k)
            if mapping === nothing
                @error("Invalid override in '$(override_file)': failed to parse entry `$(k)`")
                continue
            end

            # Next, determine if this is a hash override or a UUID/name override
            if isa(mapping, String) || isa(mapping, SHA1)
                # if this mapping is a direct mapping (e.g. a String), store it as a hash override
                local hash_str
                hash = tryparse(Base.SHA1, k)
                if hash === nothing
                    @error("Invalid override in '$(override_file)': Invalid SHA1 hash '$(k)'")
                    continue
                end

                # If this mapping is the empty string, un-override it
                if mapping == ""
                    delete!(overrides[:hash], hash)
                else
                    overrides[:hash][hash] = mapping
                end
            elseif isa(mapping, Dict)
                # Convert `k` into a uuid
                uuid = tryparse(Base.UUID, k)
                if uuid === nothing
                    @error("Invalid override in '$(override_file)': Invalid UUID '$(k)'")
                    continue
                end

                # If this mapping is itself a dict, store it as a set of UUID/artifact name overrides
                ovruuid = overrides[:UUID]::Dict{Base.UUID,Dict{String,Union{String,SHA1}}}
                if !haskey(ovruuid, uuid)
                    ovruuid[uuid] = Dict{String,Union{String,SHA1}}()
                end

                # For each name in the mapping, update appropriately
                for (name, override_value) in mapping
                    # If the mapping for this name is the empty string, un-override it
                    if override_value == ""
                        delete!(ovruuid[uuid], name)
                    else
                        # Otherwise, store it!
                        ovruuid[uuid][name] = override_value
                    end
                end
            else
                @error("Invalid override in '$(override_file)': unknown mapping type for '$(k)': $(typeof(mapping))")
            end
        end
    end

    ARTIFACT_OVERRIDES[] = overrides
    return overrides
end

# Helpers to map an override to an actual path
map_override_path(x::AbstractString) = String(x)::String
map_override_path(x::SHA1) = artifact_path(x)
map_override_path(x::Nothing) = nothing

"""
    query_override(hash::SHA1; overrides::Dict = load_overrides())

Query the loaded `<DEPOT>/artifacts/Overrides.toml` settings for artifacts that should be
redirected to a particular path or another content-hash.
"""
function query_override(hash::SHA1; overrides::Dict{Symbol,Any} = load_overrides())
    return map_override_path(get(overrides[:hash], hash, nothing))
end
function query_override(pkg::Base.UUID, artifact_name::String; overrides::Dict{Symbol,Any} = load_overrides())
    if haskey(overrides[:UUID], pkg)
        return map_override_path(get(overrides[:UUID][pkg], artifact_name, nothing))
    end
    return nothing
end

"""
    artifact_paths(hash::SHA1; honor_overrides::Bool=true)

Return all possible paths for an artifact given the current list of depots as returned
by `Pkg.depots()`.  All, some or none of these paths may exist on disk.
"""
function artifact_paths(hash::SHA1; honor_overrides::Bool=true)
    # First, check to see if we've got an override:
    if honor_overrides
        override = query_override(hash)
        if override !== nothing
            return [override]
        end
    end

    return artifacts_dirs(bytes2hex(hash.bytes))
end

"""
    artifact_path(hash::SHA1; honor_overrides::Bool=true)

Given an artifact (identified by SHA1 git tree hash), return its installation path.  If
the artifact does not exist, returns the location it would be installed to.

!!! compat "Julia 1.3"
    This function requires at least Julia 1.3.
"""
function artifact_path(hash::SHA1; honor_overrides::Bool=true)
    # Get all possible paths (rooted in all depots)
    possible_paths = artifact_paths(hash; honor_overrides=honor_overrides)

    # Find the first path that exists and return it
    for p in possible_paths
        if isdir(p)
            return p
        end
    end

    # If none exist, then just return the one that would exist within `depots1()`.
    return first(possible_paths)
end

"""
    artifact_exists(hash::SHA1; honor_overrides::Bool=true)

Returns whether or not the given artifact (identified by its sha1 git tree hash) exists
on-disk.  Note that it is possible that the given artifact exists in multiple locations
(e.g. within multiple depots).

!!! compat "Julia 1.3"
    This function requires at least Julia 1.3.
"""
function artifact_exists(hash::SHA1; honor_overrides::Bool=true)
    return any(isdir, artifact_paths(hash; honor_overrides=honor_overrides))
end

"""
    unpack_platform(entry::Dict, name::String, artifacts_toml::String)

Given an `entry` for the artifact named `name`, located within the file `artifacts_toml`,
returns the `Platform` object that this entry specifies.  Returns `nothing` on error.
"""
function unpack_platform(entry::Dict{String,Any}, name::String,
                         artifacts_toml::String)::Union{Nothing,Platform}
    if !haskey(entry, "os")
        @error("Invalid artifacts file at '$(artifacts_toml)': platform-specific artifact entry '$name' missing 'os' key")
        return nothing
    end

    if !haskey(entry, "arch")
        @error("Invalid artifacts file at '$(artifacts_toml)': platform-specific artifact entrty '$name' missing 'arch' key")
        return nothing
    end

    # Collect all String-valued mappings in `entry` and use them as tags
    tags = Dict{String, String}()
    for (k, v) in entry
        if v isa String
            tags[k] = v
        end
    end
    # Removing some known entries that shouldn't be passed through `tags`
    delete!(tags, "os")
    delete!(tags, "arch")
    delete!(tags, "git-tree-sha1")
    return Platform(entry["arch"], entry["os"], tags)
end

function pack_platform!(meta::Dict, p::AbstractPlatform)
    for (k, v) in tags(p)
        if v !== nothing
            meta[k] = v
        end
    end
    return meta
end

"""
    load_artifacts_toml(artifacts_toml::String;
                        pkg_uuid::Union{UUID,Nothing}=nothing)

Loads an `(Julia)Artifacts.toml` file from disk.  If `pkg_uuid` is set to the `UUID` of the
owning package, UUID/name overrides stored in a depot `Overrides.toml` will be resolved.
"""
function load_artifacts_toml(artifacts_toml::String;
                             pkg_uuid::Union{Base.UUID,Nothing} = nothing)
    artifact_dict = parse_toml(artifacts_toml)

    # Process overrides for this `pkg_uuid`
    process_overrides(artifact_dict, pkg_uuid)
    return artifact_dict
end

"""
    process_overrides(artifact_dict::Dict, pkg_uuid::Base.UUID)

When loading an `Artifacts.toml` file, we must check `Override.toml` files to see if any
of the artifacts within it have been overridden by UUID.  If they have, we honor the
overrides by inspecting the hashes of the targeted artifacts, then overriding them to
point to the given override, punting the actual redirection off to the hash-based
override system.  This does not modify the `artifact_dict` object, it merely dynamically
adds more hash-based overrides as `Artifacts.toml` files that are overridden are loaded.
"""
function process_overrides(artifact_dict::Dict, pkg_uuid::Base.UUID)
    # Insert just-in-time hash overrides by looking up the names of anything we need to
    # override for this UUID, and inserting new overrides for those hashes.
    overrides = load_overrides()
    if haskey(overrides[:UUID], pkg_uuid)
        pkg_overrides = overrides[:UUID][pkg_uuid]

        for name in keys(artifact_dict)
            # Skip names that we're not overriding
            if !haskey(pkg_overrides, name)
                continue
            end

            # If we've got a platform-specific friend, override all hashes:
            if isa(artifact_dict[name], Array)
                for entry in artifact_dict[name]
                    hash = SHA1(entry["git-tree-sha1"])
                    overrides[:hash][hash] = overrides[:UUID][pkg_uuid][name]
                end
            elseif isa(artifact_dict[name], Dict)
                hash = SHA1(artifact_dict[name]["git-tree-sha1"])
                overrides[:hash][hash] = overrides[:UUID][pkg_uuid][name]
            end
        end
    end
    return artifact_dict
end

# If someone tries to call process_overrides() with `nothing`, do exactly that
process_overrides(artifact_dict::Dict, pkg_uuid::Nothing) = nothing

"""
    artifact_meta(name::String, artifacts_toml::String;
                  platform::AbstractPlatform = HostPlatform(),
                  pkg_uuid::Union{Base.UUID,Nothing}=nothing)

Get metadata about a given artifact (identified by name) stored within the given
`(Julia)Artifacts.toml` file.  If the artifact is platform-specific, use `platform` to choose the
most appropriate mapping.  If none is found, return `nothing`.

!!! compat "Julia 1.3"
    This function requires at least Julia 1.3.
"""
function artifact_meta(name::String, artifacts_toml::String;
                       platform::AbstractPlatform = HostPlatform(),
                       pkg_uuid::Union{Base.UUID,Nothing}=nothing)
    if !isfile(artifacts_toml)
        return nothing
    end

    # Parse the toml of the artifacts_toml file
    artifact_dict = load_artifacts_toml(artifacts_toml; pkg_uuid=pkg_uuid)
    return artifact_meta(name, artifact_dict, artifacts_toml; platform=platform)
end

function artifact_meta(name::String, artifact_dict::Dict, artifacts_toml::String;
                       platform::AbstractPlatform = HostPlatform())
    if !haskey(artifact_dict, name)
        return nothing
    end
    meta = artifact_dict[name]

    # If it's an array, find the entry that best matches our current platform
    if isa(meta, Vector)
        dl_dict = Dict{AbstractPlatform,Dict{String,Any}}()
        for x in meta
            x::Dict{String}
            dl_dict[unpack_platform(x, name, artifacts_toml)] = x
        end
        meta = select_platform(dl_dict, platform)
    # If it's NOT a dict, complain
    elseif !isa(meta, Dict)
        @error("Invalid artifacts file at $(artifacts_toml): artifact '$name' malformed, must be array or dict!")
        return nothing
    end

    # This is such a no-no, we are going to call it out right here, right now.
    if meta !== nothing && !haskey(meta, "git-tree-sha1")
        @error("Invalid artifacts file at $(artifacts_toml): artifact '$name' contains no `git-tree-sha1`!")
        return nothing
    end

    # Return the full meta-dict.
    return meta
end

"""
    artifact_hash(name::String, artifacts_toml::String;
                  platform::AbstractPlatform = HostPlatform())

Thin wrapper around `artifact_meta()` to return the hash of the specified, platform-
collapsed artifact.  Returns `nothing` if no mapping can be found.

!!! compat "Julia 1.3"
    This function requires at least Julia 1.3.
"""
function artifact_hash(name::String, artifacts_toml::String;
                       platform::AbstractPlatform = HostPlatform(),
                       pkg_uuid::Union{Base.UUID,Nothing}=nothing)::Union{Nothing, SHA1}
    meta = artifact_meta(name, artifacts_toml; platform=platform)
    if meta === nothing
        return nothing
    end

    return SHA1(meta["git-tree-sha1"])
end

function select_downloadable_artifacts(artifact_dict::Dict, artifacts_toml::String;
                                       platform::AbstractPlatform = HostPlatform(),
                                       pkg_uuid::Union{Nothing,Base.UUID} = nothing,
                                       include_lazy::Bool = false)
    artifacts = Dict{String,Any}()
    for name in keys(artifact_dict)
        # Get the metadata about this name for the requested platform
        meta = artifact_meta(name, artifact_dict, artifacts_toml; platform=platform)

        # If there are no instances of this name for the desired platform, skip it
        # Also skip if there's no `download` stanza (e.g. it's only a local artifact)
        # or if it's lazy and we're not explicitly looking for lazy artifacts.
        if meta === nothing || !haskey(meta, "download") || (get(meta, "lazy", false) && !include_lazy)
            continue
        end

        # Else, welcome it into the meta-fold
        artifacts[name] = meta
    end
    return artifacts
end

"""
    select_downloadable_artifacts(artifacts_toml::String;
                                  platform = HostPlatform,
                                  include_lazy = false,
                                  pkg_uuid = nothing)

Returns a dictionary where every entry is an artifact from the given `Artifacts.toml`
that should be downloaded for the requested platform.  Lazy artifacts are included if
`include_lazy` is set.
"""
function select_downloadable_artifacts(artifacts_toml::String;
                                       platform::AbstractPlatform = HostPlatform(),
                                       include_lazy::Bool = false,
                                       pkg_uuid::Union{Nothing,Base.UUID} = nothing)
    if !isfile(artifacts_toml)
        return Dict{String,Any}()
    end
    artifact_dict = load_artifacts_toml(artifacts_toml; pkg_uuid=pkg_uuid)
    return select_downloadable_artifacts(artifact_dict, artifacts_toml; platform, pkg_uuid, include_lazy)
end

"""
    find_artifacts_toml(path::String)

Given the path to a `.jl` file, (such as the one returned by `__source__.file` in a macro
context), find the `(Julia)Artifacts.toml` that is contained within the containing project (if it
exists), otherwise return `nothing`.

!!! compat "Julia 1.3"
    This function requires at least Julia 1.3.
"""
function find_artifacts_toml(path::String)
    if !isdir(path)
        path = dirname(path)
    end

    # Run until we hit the root directory.
    while dirname(path) != path
        for f in artifact_names
            artifacts_toml_path = joinpath(path, f)
            if isfile(artifacts_toml_path)
                return abspath(artifacts_toml_path)
            end
        end

        # Does a `(Julia)Project.toml` file exist here, in the absence of an Artifacts.toml?
        # If so, stop the search as we've probably hit the top-level of this package,
        # and we don't want to escape out into the larger filesystem.
        for f in Base.project_names
            if isfile(joinpath(path, f))
                return nothing
            end
        end

        # Move up a directory
        path = dirname(path)
    end

    # We never found anything, just return `nothing`
    return nothing
end

# We do this to avoid doing the `joinpath()` work if we don't have to, and also to
# avoid a trailing slash due to `joinpath()`'s habit of including one when the last
# argument is the empty string.
function jointail(dir, tail)
    if !isempty(tail)
        return joinpath(dir, tail)
    else
        return dir
    end
end

function _artifact_str(__module__, artifacts_toml, name, path_tail, artifact_dict, hash, platform, @nospecialize(lazyartifacts))
    moduleroot = Base.moduleroot(__module__)
    if haskey(Base.module_keys, moduleroot)
        # Process overrides for this UUID, if we know what it is
        process_overrides(artifact_dict, Base.module_keys[moduleroot].uuid)
    end

    # If the artifact exists, we're in the happy path and we can immediately
    # return the path to the artifact:
    dirs = artifact_paths(hash; honor_overrides=true)
    for dir in dirs
        if isdir(dir)
            return jointail(dir, path_tail)
        end
    end

    # If not, try determining what went wrong:
    meta = artifact_meta(name, artifact_dict, artifacts_toml; platform)
    if meta !== nothing && get(meta, "lazy", false)
        if lazyartifacts isa Module && isdefined(lazyartifacts, :ensure_artifact_installed)
            if nameof(lazyartifacts) in (:Pkg, :Artifacts)
                Base.depwarn("using Pkg instead of using LazyArtifacts is deprecated", :var"@artifact_str", force=true)
            end
            return jointail(lazyartifacts.ensure_artifact_installed(string(name), artifacts_toml; platform), path_tail)
        end
        error("Artifact $(repr(name)) is a lazy artifact; package developers must call `using LazyArtifacts` in $(__module__) before using lazy artifacts.")
    end

    path_str = if length(dirs) == 1
        "path \"$(first(dirs))\". "
    else
        string("paths:\n", join("  " .* contractuser.(dirs), '\n'), '\n')
    end

    suggestion_str = if query_override(hash) !== nothing
        "Check that your `Overrides.toml` file is correct (https://pkgdocs.julialang.org/v1/artifacts/#Overriding-artifact-locations)."
    else
        "Try `using Pkg; Pkg.instantiate()` to re-install all missing resources."
    end

    error("Artifact $(repr(name)) was not found by looking in the $(path_str)$suggestion_str")
end

raw"""
    split_artifact_slash(name::String)

Splits an artifact indexing string by path deliminters, isolates the first path element,
returning that and the `joinpath()` of the remaining arguments.  This normalizes all path
separators to the native path separator for the current platform.  Examples:

# Examples
```jldoctest; setup = :(using Artifacts: split_artifact_slash)
julia> split_artifact_slash("Foo")
("Foo", "")

julia> ret = split_artifact_slash("Foo/bar/baz.so");

julia> if Sys.iswindows()
            ret == ("Foo", "bar\\baz.so")
       else
            ret == ("Foo", "bar/baz.so")
       end
true

julia> ret = split_artifact_slash("Foo\\bar\\baz.so");

julia> if Sys.iswindows()
            ret == ("Foo", "bar\\baz.so")
       else
            ret == ("Foo", "bar/baz.so")
       end
true
```
"""
function split_artifact_slash(name::String)
    split_name = split(name, r"(/|\\)")
    if length(split_name) == 1
        return (split_name[1], "")
    else
        return (split_name[1], joinpath(split_name[2:end]...))
    end
end

"""
    artifact_slash_lookup(name::String, atifact_dict::Dict,
                          artifacts_toml::String, platform::Platform)

Returns `artifact_name`, `artifact_path_tail`, and `hash` by looking the results up in
the given `artifacts_toml`, first extracting the name and path tail from the given `name`
to support slash-indexing within the given artifact.
"""
function artifact_slash_lookup(name::String, artifact_dict::Dict,
                               artifacts_toml::String, platform::Platform)
    artifact_name, artifact_path_tail = split_artifact_slash(name)

    meta = artifact_meta(artifact_name, artifact_dict, artifacts_toml; platform)
    if meta === nothing
        error("Cannot locate artifact '$(name)' for $(triplet(platform)) in '$(artifacts_toml)'")
    end
    hash = SHA1(meta["git-tree-sha1"])
    return artifact_name, artifact_path_tail, hash
end

"""
    macro artifact_str(name)

Return the on-disk path to an artifact. Automatically looks the artifact up by
name in the project's `(Julia)Artifacts.toml` file. Throws an error on if the
requested artifact is not present. If run in the REPL, searches for the toml
file starting in the current directory, see `find_artifacts_toml()` for more.

If the artifact is marked "lazy" and the package has `using LazyArtifacts`
defined, the artifact will be downloaded on-demand with `Pkg` the first time
this macro tries to compute the path. The files will then be left installed
locally for later.

If `name` contains a forward or backward slash, all elements after the first slash will
be taken to be path names indexing into the artifact, allowing for an easy one-liner to
access a single file/directory within an artifact.  Example:

    ffmpeg_path = @artifact"FFMPEG/bin/ffmpeg"

!!! compat "Julia 1.3"
    This macro requires at least Julia 1.3.

!!! compat "Julia 1.6"
    Slash-indexing requires at least Julia 1.6.
"""
macro artifact_str(name, platform=nothing)
    # Find Artifacts.toml file we're going to load from
    srcfile = string(__source__.file)
    if ((isinteractive() && startswith(srcfile, "REPL[")) || (!isinteractive() && srcfile == "none")) && !isfile(srcfile)
        srcfile = pwd()
    end
    local artifacts_toml = find_artifacts_toml(srcfile)
    if artifacts_toml === nothing
        error(string(
            "Cannot locate '(Julia)Artifacts.toml' file when attempting to use artifact '",
            name,
            "' in '",
            __module__,
            "'",
        ))
    end

    # Load Artifacts.toml at compile time, so that we don't have to use `__source__.file`
    # at runtime, which gets stale if the `.ji` file is relocated.
    local artifact_dict = load_artifacts_toml(artifacts_toml)

    # Invalidate calling .ji file if Artifacts.toml file changes
    Base.include_dependency(artifacts_toml)

    # Check if the user has provided `LazyArtifacts`, and thus supports lazy artifacts
    # If not, check to see if `Pkg` or `Pkg.Artifacts` has been imported.
    lazyartifacts = nothing
    for module_name in (:LazyArtifacts, :Pkg, :Artifacts)
        if isdefined(__module__, module_name)
            lazyartifacts = GlobalRef(__module__, module_name)
            break
        end
    end

    # If `name` is a constant, (and we're using the default `Platform`) we can actually load
    # and parse the `Artifacts.toml` file now, saving the work from runtime.
    if isa(name, AbstractString) && platform === nothing
        # To support slash-indexing, we need to split the artifact name from the path tail:
        platform = HostPlatform()
        artifact_name, artifact_path_tail, hash = artifact_slash_lookup(name, artifact_dict, artifacts_toml, platform)
        return quote
            Base.invokelatest(_artifact_str, $(__module__), $(artifacts_toml), $(artifact_name), $(artifact_path_tail), $(artifact_dict), $(hash), $(platform), $(lazyartifacts))::String
        end
    else
        if platform === nothing
            platform = :($(HostPlatform)())
        end
        return quote
            local platform = $(esc(platform))
            local artifact_name, artifact_path_tail, hash = artifact_slash_lookup($(esc(name)), $(artifact_dict), $(artifacts_toml), platform)
            Base.invokelatest(_artifact_str, $(__module__), $(artifacts_toml), artifact_name, artifact_path_tail, $(artifact_dict), hash, platform, $(lazyartifacts))::String
        end
    end
end

# Support `AbstractString`s, but avoid compilers needing to track backedges for callers
# of these functions in case a user defines a new type that is `<: AbstractString`
with_artifacts_directory(f::Function, artifacts_dir::AbstractString) =
    with_artifacts_directory(f, String(artifacts_dir)::String)
query_override(pkg::Base.UUID, artifact_name::AbstractString; overrides::Dict=load_overrides()) =
    query_override(pkg, String(artifact_name)::String; overrides=convert(Dict{Symbol, Any}(overrides)))
unpack_platform(entry::Dict, name::AbstractString, artifacts_toml::AbstractString) =
    unpack_platform(convert(Dict{String, Any}, entry), String(name)::String, String(artifacts_toml)::String)
load_artifacts_toml(artifacts_toml::AbstractString; kwargs...) =
    load_artifacts_toml(String(artifacts_toml)::String; kwargs...)
artifact_meta(name::AbstractString, artifacts_toml::AbstractString; kwargs...) =
    artifact_meta(String(name)::String, String(artifacts_toml)::String; kwargs...)
artifact_meta(name::AbstractString, artifact_dict::Dict, artifacts_toml::AbstractString; kwargs...) =
    artifact_meta(String(name)::String, artifact_dict, String(artifacts_toml)::String; kwargs...)
artifact_hash(name::AbstractString, artifacts_toml::AbstractString; kwargs...) =
    artifact_hash(String(name)::String, String(artifacts_toml)::String; kwargs...)
select_downloadable_artifacts(artifact_dict::Dict, artifacts_toml::AbstractString; kwargs...) =
    select_downloadable_artifacts(artifact_dict, String(artifacts_toml)::String, kwargs...)
select_downloadable_artifacts(artifacts_toml::AbstractString; kwargs...) =
    select_downloadable_artifacts(String(artifacts_toml)::String, kwargs...)
find_artifacts_toml(path::AbstractString) =
    find_artifacts_toml(String(path)::String)
split_artifact_slash(name::AbstractString) =
    split_artifact_slash(String(name)::String)
artifact_slash_lookup(name::AbstractString, artifact_dict::Dict, artifacts_toml::AbstractString) =
    artifact_slash_lookup(String(name)::String, artifact_dict, String(artifacts_toml)::String)

# Precompilation to reduce latency
precompile(load_artifacts_toml, (String,))
precompile(NamedTuple{(:pkg_uuid,)}, (Tuple{Base.UUID},))
precompile(Core.kwfunc(load_artifacts_toml), (NamedTuple{(:pkg_uuid,), Tuple{Base.UUID}}, typeof(load_artifacts_toml), String))

end # module Artifacts

max_stars_repo_name: Physolia/julia

max_issues_repo_name: Eurus-Holmes/julia

max_forks_repo_name: Nexuscompute/julia


text:
function get_attitude_provider(obj::StateMapper)
    return jcall(obj, "getAttitudeProvider", AttitudeProvider, ())
end

function get_frame(obj::StateMapper)
    return jcall(obj, "getFrame", Frame, ())
end

function get_mu(obj::StateMapper)
    return jcall(obj, "getMu", jdouble, ())
end

function get_orbit_type(obj::StateMapper)
    return jcall(obj, "getOrbitType", OrbitType, ())
end

function get_position_angle_type(obj::StateMapper)
    return jcall(obj, "getPositionAngleType", PositionAngle, ())
end

function get_reference_date(obj::StateMapper)
    return jcall(obj, "getReferenceDate", AbsoluteDate, ())
end

function map_array_to_state(obj::StateMapper, arg0::AbsoluteDate, arg1::Vector{jdouble}, arg2::Vector{jdouble}, arg3::PropagationType)
    return jcall(obj, "mapArrayToState", SpacecraftState, (AbsoluteDate, Vector{jdouble}, Vector{jdouble}, PropagationType), arg0, arg1, arg2, arg3)
end

function map_array_to_state(obj::StateMapper, arg0::jdouble, arg1::Vector{jdouble}, arg2::Vector{jdouble}, arg3::PropagationType)
    return jcall(obj, "mapArrayToState", SpacecraftState, (jdouble, Vector{jdouble}, Vector{jdouble}, PropagationType), arg0, arg1, arg2, arg3)
end

function map_date_to_double(obj::StateMapper, arg0::AbsoluteDate)
    return jcall(obj, "mapDateToDouble", jdouble, (AbsoluteDate,), arg0)
end

function map_double_to_date(obj::StateMapper, arg0::jdouble)
    return jcall(obj, "mapDoubleToDate", AbsoluteDate, (jdouble,), arg0)
end

function map_double_to_date(obj::StateMapper, arg0::jdouble, arg1::AbsoluteDate)
    return jcall(obj, "mapDoubleToDate", AbsoluteDate, (jdouble, AbsoluteDate), arg0, arg1)
end

function map_state_to_array(obj::StateMapper, arg0::SpacecraftState, arg1::Vector{jdouble}, arg2::Vector{jdouble})
    return jcall(obj, "mapStateToArray", void, (SpacecraftState, Vector{jdouble}, Vector{jdouble}), arg0, arg1, arg2)
end

function set_position_angle_type(obj::StateMapper)
    return jcall(obj, "setPositionAngleType", void, ())
end


text:
######################################################################
# reading the data, matching up people, writing out serialized data
#
# NOTE: You don't need to rerun this for exploratory work.
# It will be rerun when Julia is updated.
#
# NOTE: compressing via GZip gives ≈1/6 file size and +20% read time.
######################################################################
using AMDB
using GZip

data_dir = "/Data/AMDB/RawTables/"

srand(5447225802448749441)      # random seed for reproducible data

# read person spells
records = process_compressed_database(joinpath(data_dir, "hv_pn.csv.gz")) do io
    read_person_spells!(io, 0.01)
end

# remove records with no person spells
records = filter((id,record) -> !isempty(record.person_spells), records)

# read AMP data
process_compressed_database(joinpath(data_dir, "hv_epi_uni_roh.csv.gz")) do io
    read_AMP_spells!(io, records)
end

# read VMZ v2 data
process_compressed_database(joinpath(data_dir, "vmz_v2_int.csv.gz")) do io
    read_VMZ_spells!(io, records)
end

# write out data
GZip.open(normpath(joinpath(data_dir,
                            "../subsample/AMDB_subsample.jls.gz")), "w") do io
    serialize(io, records)
end


text:
@testset "Slider" begin
    s = slider("a", [0,0], [3,0], [0,1.5,3])
    @test s.name == "a"
    @test s.vals == [[0,0], [3,0], [0,1.5,3]]

    @test_throws AssertionError slider("a", [0, 1])
    @test_throws AssertionError slider("a", [[1,2], [2,3,4], [1,2]])

    # midpoint
    @test val(s) == 1.5

    b = board("brd")
    ss = J.str(s, b)
    @test ss == "var a=brd.create('slider', [[0.0,0.0],[3.0,0.0],[0.0,1.5,3.0]], {});"
end


text:
function knuthshuffle!(r::AbstractRNG, v::AbstractVector)
    for i in length(v):-1:2
        j = rand(r, 1:i)
        v[i], v[j] = v[j], v[i]
    end
    return v
end
knuthshuffle!(v::AbstractVector) = knuthshuffle!(Base.Random.GLOBAL_RNG, v)

v = collect(1:20)
println("# v = $v\n   -> ", knuthshuffle!(v))

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using IterTools
using BenchmarkTools

struct Window{T}
    xs::T
    width::Int
    offset::Int
    equal_size::Bool #default = true
end

function Window(xs, width::Int, offset::Int; equal_size = true)
    if width < 0
        error("width must be greater than 0")
    end
    if offset < 0
        error("offset mus be greater than 0")
    end
    return Window(xs, width, offset, equal_size)
end

Base.IteratorSize(it::Window) = Base.SizeUnknown()
Base.eltype(it::Window{T}) where T = Vector{eltype(T)}

function Base.iterate(it::Window{T}, state = 1) where T <: AbstractVector
    if state + it.width - 1 <= length(it.xs)
        return (it.xs[state:(state + it.width - 1)], state + it.offset)
    elseif !it.equal_size && state < length(it.xs)
        return (it.xs[state:end], state + it.offset)
    else
        return nothing
    end
end

event_times = 1:10

pw = Window(event_times, 3, 1)
pw = Window(enumerate(1:100), 3, 1)
pw.xs[1]
for w in pw
    @show w
end
collect(pw)

collect(partition(event_times, 3, 1))

@btime collect(Window($event_times, 2, 1));
@btime collect(partition($event_times, 2, 1));

@time collect(Window(1:1000, 2, 1));
@time collect(partition(1:1000, 2, 1));

@code_warntype collect(Window(event_times, 2, 1));
@code_warntype collect(partition(event_times, 2, 1));


text:
struct FileHolder
    fileProcessor::Dict{String, Function}

    FileHolder(processor::Dict{String, Function}) = new(processor)
end

function FileHolder(filenames::Vector{String})
    processor(line::String) = split(line, ' ')
    FileHolder(Dict{String, Function}(filename => processor) for filename in filenames)
end

struct OutSetting
    outFilename::String
    outFileFormat::Symbol

    OutSetting(filename::String) = new(filename, :json)
end

__idfCaculator(value, documentCount) = log(documentCount / (value != 0 ? value : 1))

mutable struct IDF
    idfDict::Dict{String, Int64}
    idfCalculator::Function
    documentCount::Int64

    fileHolder::FileHolder
    outSetting::OutSetting
    IDF(fileHolder::FileHolder, outSetting::OutSetting=OutSetting("out.txt"); idfCalculator::Function=__idfCaculator) =
        new(Dict{String, Int64}(), idfCalculator, 0, fileHolder, outSetting)
end

function fit!(self::IDF)
    # TODO: parallize it
    for (filename, processor) in self.fileHolder.fileProcessor
        open(filename) do file
            statisticIDF!(self, file, processor)
        end
    end
    return self
end

function statisticIDF!(self::IDF, file, processor::Function)
    for (i, line) in enumerate(eachline(file))
        tokens = processor(line)
        for token in Set(tokens)
            self.idfDict[token] = get(self.idfDict, token, 0) + 1
        end
        self.documentCount += 1

        if i % 1e4 == 0
            println(i)
        end
    end
end

function transform()

end

getIdfValue(self::IDF, key::String) = self.caculator(get(self.idfDict, key, 0), documentCount)

function main()
    # 空格分隔后的，分隔词性，取词
    processor(line::String) = map(t -> t[1], map(s -> split(s, "\x01"), split(line, " ")))
    fh = FileHolder(Dict{String, Function}("filePath" => processor))
    idf = IDF(fh)

    fit!(idf)
end

#main()

max_stars_repo_name: gajanlee/toys

max_issues_repo_name: gajanlee/toys

max_forks_repo_name: gajanlee/toys


text:
using PyPlot, Statistics
# (c) Ludovic Raess (ludovic.rass@gmail.com)
@views function viscous_glacier2D()
	# physics - scales
	α       = 10.0                 # slope angle, degrees
	lx, ly  = 300.0, 100.0         # domain length X, Y, m
	ρg      = 900*9.8              # gravity acceleration, m/s^2
	μ       = 1e12                 # viscosity, Pa.s
	# numerics
	nx, ny  = 61, 21               # grid points X, Y
	niτ     = 14000                # number of iters
	# preprocessing
	dx, dy  = lx/(nx-1), ly/(ny-1) # grid spacing X, Y
	dτv     = min(dx,dy)^2/μ/4.1/2 # timestep
	dτp     = 4.1*μ/max(nx,ny)
	(xc, yc ) = (LinRange(-lx/2, lx/2, nx), LinRange(-ly/2, ly/2, ny))
	(Xc2,Yc2) = ([x for x=xc,y=yc], [y for x=xc,y=yc])
	(X2r,Y2r) = (Xc2*cosd(α) - sind(α)*Yc2, Xc2*sind(α) + cosd(α)*Yc2)
	# initial conditions
	Pt      = zeros(Float64, nx  ,ny  )
	∇V      = zeros(Float64, nx  ,ny  )
	Vx      = zeros(Float64, nx+1,ny  )
	Vy      = zeros(Float64, nx  ,ny+1)
	dVxdτ   = zeros(Float64, nx-1,ny-2)
	dVydτ   = zeros(Float64, nx-2,ny-1)
	σ_xx    = zeros(Float64, nx  ,ny  )
	σ_yy    = zeros(Float64, nx  ,ny  )
	τ_xy    = zeros(Float64, nx-1,ny-1)
	# action - iteration loop
	for iτ = 1:niτ
		# pressure and stress
		∇V    .= diff(Vx,dims=1)./dx .+ diff(Vy,dims=2)./dy
		Pt    .= Pt .- ∇V.*dτp
		σ_xx  .= 2.0*μ.*(diff(Vx,dims=1)./dx .- 1/3*∇V) .- Pt
		σ_yy  .= 2.0*μ.*(diff(Vy,dims=2)./dy .- 1/3*∇V) .- Pt
		τ_xy  .= μ.*(diff(Vx[2:end-1,:],dims=2)./dy .+ diff(Vy[:,2:end-1],dims=1)./dx)
		# free surface top boundary
		Pt[:,end]   .= 0
		σ_xx[:,end] .= 0
		τ_xy[:,end] .= 1/3*τ_xy[:,end-1]
		# flow velocities Vx and Vy
		dVxdτ .= diff(σ_xx[:,2:end-1],dims=1)./dx .+ diff(τ_xy,dims=2)./dy .- sind(α).*ρg
		dVydτ .= diff(σ_yy[2:end-1,:],dims=2)./dy .+ diff(τ_xy,dims=1)./dx .- cosd(α).*ρg
		Vx[2:end-1,2:end-1] .= Vx[2:end-1,2:end-1] .+ dVxdτ.*dτv
		Vy[2:end-1,2:end-1] .= Vy[2:end-1,2:end-1] .+ dVydτ.*dτv
		# no sliding bottom boundary
		Vy[:,1] .= -Vy[:,2]
		# visualisation
		if mod(iτ,1000)==0 err = maximum([mean(abs.(∇V[:,1:end-1])), mean(abs.(dVxdτ)), mean(abs.(dVydτ))]); println("error = $err") end
	end
	# visu
	figure("pyplot_figure",figsize=(5.7,5.5)), clf()
	st = 4; s2d = 60*60*24
	(Xp,  Yp ) = (X2r[1:st:end,1:st:end], Y2r[1:st:end,1:st:end])
	(Vxp, Vyp) = (0.5*(Vx[1:st:end-1,1:st:end  ]+Vx[2:st:end,1:st:end]), 0.5*(Vy[1:st:end  ,1:st:end-1]+Vy[1:st:end,2:st:end]))
	subplot(311), pcolor(X2r,Y2r, 1e-3*Pt), plt.axis("off"), plt.title("pressure [kPa]"), plt.colorbar(fraction=0.02, pad=0.03)
	quiver(Xp, Yp, Vxp, Vyp, pivot="mid", color="white")
	subplot(312), pcolor(X2r,Y2r, s2d*Vx[2:end,:]), plt.axis("off"), plt.title("Vel-x [m/d]"), plt.colorbar(fraction=0.02, pad=0.03)
	subplot(313), pcolor(X2r,Y2r, s2d*Vy[:,2:end]), plt.axis("off"), plt.title("Vel-y [m/d]"), plt.colorbar(fraction=0.02, pad=0.03)
	return
end
viscous_glacier2D()


text:
"""
    function addUpdateMediumInteraction_!(p::Program_,abm::Agent,space::SimulationSpace,platform::String)

Generate the functions related with Agent Medium Interaction.
"""
function addUpdateMediumInteraction_!(p::Program_,abm::Agent,space::SimulationSpace,platform::String)

    if "UpdateMediumInteraction" in keys(p.agent.declaredUpdates)

        code = p.agent.declaredUpdates["UpdateMediumInteraction"]

        for (i,j) in enumerate(p.agent.declaredSymbols["Medium"])
            if p.agent.dims == 1
                code = postwalk(x->@capture(x,y_) && y == j ? :(mediumVCopy[idMediumV[ic1_, 1],$i]) : x, code)
            elseif p.agent.dims == 2
                code = postwalk(x->@capture(x,y_) && y == j ? :(mediumVCopy[idMediumV[ic1_, 1],idMediumV[ic1_, 2],$i]) : x, code)
            elseif p.agent.dims == 3
                code = postwalk(x->@capture(x,y_) && y == j ? :(mediumVCopy[idMediumV[ic1_, 1],idMediumV[ic1_, 2],idMediumV[ic1_, 3],$i]) : x, code)
            end
        end

        #Construct functions
        f = simpleFirstLoopWrapInFunction_(platform,:mediumInteractionStep_!,code)
        f = vectorize_(abm,f,p)

        push!(p.declareF.args,
            f)

    end

    return nothing

end

max_stars_repo_name: gatocor/embryogenesisJulia

max_issues_repo_name: gatocor/embryogenesisJulia

max_forks_repo_name: gatocor/embryogenesisJulia


text:
using TensorFlow

sess = TensorFlow.Session()

x = TensorFlow.constant(Float64[1,2])
y = TensorFlow.Variable(Float64[3,4])
z = TensorFlow.placeholder(Float64)

w = exp(x + z + -y)

run(sess, TensorFlow.global_variables_initializer())
res = run(sess, w, Dict(z=>Float64[1,2]))
Base.Test.@test res[1] ≈ exp(-1)


text:
function part1()
	lines = readlines("inputD10.txt")
	scores = Dict{Char, Int}(')' => 3, ']' => 57, '}' => 1197, '>' => 25137)
	opentoclose = Dict{Char, Char}('(' => ')', '[' => ']', '{' => '}', '<' => '>')

	score = 0

	for line in lines
		brackets = Char[]
		for char in line
			if char in keys(opentoclose)
				push!(brackets, char)
			else
				lastopen = pop!(brackets)
				if opentoclose[lastopen] != char
					score += scores[char]
					break
				end
			end
		end
	end

	return score
end

function part2()
	lines = readlines("inputD10.txt")
	opentoclose = Dict{Char, Char}('(' => ')', '[' => ']', '{' => '}', '<' => '>')

	closescores = Dict{Char, Int}(')' => 1, ']' => 2, '}' => 3, '>' => 4)
	linescores = Int[]
	for line in lines
		brackets = Char[]
		corrupted = false
		for char in line
			if char in keys(opentoclose)
				push!(brackets, char)
			else
				lastopen = pop!(brackets)
				if opentoclose[lastopen] != char
					corrupted = true
					break
				end
			end
		end

		if !corrupted
			linescore = 0
			for bracket in reverse(brackets)
				linescore *= 5
				linescore += closescores[opentoclose[bracket]]
			end
			push!(linescores, linescore)
		end
	end

	return sort(linescores)[div(length(linescores) + 1, 2)]
end

println(part1())
println(part2())


text:
function model_selection(input::ModelSelectionInput)

	# For logging
	log_prefix = "GpABC model selection simulation "
	if isa(input, EmulatedModelSelectionInput)
		log_prefix = "GpABC model selection emulation "
	end

	@info "$(log_prefix) Population 1 - ABC Rejection ϵ = $(input.threshold_schedule[1])"

	tracker = initialise_modelselection(input)

	if all_models_dead(tracker)
		@warn "No particles were accepted in population 1 with threshold $(input.threshold_schedule[1])- terminating model selection algorithm"
		return build_modelselection_output(tracker, false)
	end

	if all_but_one_models_dead(tracker)
		@warn "All but one model is dead after population 1 - terminating model selection algorithm"
		return build_modelselection_output(tracker, true)
	end

	for i in 2:length(input.threshold_schedule)
		@info "$(log_prefix) Population $i - ABC SMC ϵ = $(input.threshold_schedule[i])"

		iterate_modelselection!(tracker, input.threshold_schedule[i])

		# Avoid infinite loop if no particles are accepted
		if all_models_dead(tracker)
			@warn "No particles were accepted in population $i with threshold $(input.threshold_schedule[i])- terminating model selection algorithm"
			return build_modelselection_output(tracker, false)
		end

		if all_but_one_models_dead(tracker)
			@warn "All but one model is dead after population $i - terminating model selection algorithm"
			return build_modelselection_output(tracker, true)
		end
	end

	return build_modelselection_output(tracker, true)
end

# Initialises the model selection run and runs the first (rejection) population
function initialise_modelselection(input::SimulatedModelSelectionInput)

	#
	# Check input sizes
	#
	if span(input.model_prior) != input.M
		throw(ArgumentError("There are $(input.M) models but the span of the model prior support is $(span(input.model_prior))"))
	end

	if length(input.distance_simulation_input) != input.M
		throw(ArgumentError("There are $(input.M) models but $(length(input.simulator_functions)) simulator functions"))
	end

	if length(input.parameter_priors) != input.M
		throw(ArgumentError("There are $(input.M) models but $(length(input.parameter_priors)) sets of parameter priors"))
	end

	total_n_accepted = 0
	n_iterations = 1

	#
	# Initialise arrays that will track rejection ABC run for each model - these
	# will be used to create SimulatedCandidateModelTrackers after the rejection ABC run
	#
	cm_trackers = [CandidateModelTracker(length(input.parameter_priors[m]))
								for m in 1:input.M]

    # Initialise variables to hold ABC results
    parameters = zeros(0)
    distance = 0.0
    weight_value = 0.0

	#
	# Compute first population using rejection-ABC
	#
	while total_n_accepted < input.n_particles && n_iterations <= input.max_iter

		m = rand(input.model_prior)
		parameters, weight_value = generate_parameters(input.parameter_priors[m])
		try
			distance = simulate_distance(parameters, input.distance_simulation_input[m])
        catch e
            if isa(e, DimensionMismatch)
                # This prevents the whole code from failing if there is a problem
                # solving the differential equation(s). The exception is thrown by the
                # distance function
                @warn "The summarised simulated data does not have the same size as the summarised reference data. If this is not happening at every iteration it may be due to the behaviour of OrdinaryDiffEq::solve - please check for related warnings. Continuing to the next iteration."
                n_iterations += 1
                continue
            else
                throw(e)
            end
        end

        cm_trackers[m].n_tries += 1

		# If particle accepted
		if distance[1] < input.threshold_schedule[1]
			total_n_accepted += 1
			update_candidatemodeltracker!(cm_trackers[m], parameters, distance[1], weight_value)
		end

		n_iterations += 1
	end

	if n_iterations == input.max_iter+1
		@warn "Simulated model selection reached maximum number of iterations ($(input.max_iter)) on the first population - consider trying more iterations."
	end

	@info join(["GpABC model selection simulation",
		"Completed $(n_iterations-1) iterations, accepting $total_n_accepted particles in total.",
		"Number of accepted parameters by model:",
		join([string("Model ", m, ": ",  cm_trackers[m].n_accepted) for m in 1:input.M], "\t")
		], "\n")

	smc_trackers = [SimulatedABCSMCTracker(
		length(input.parameter_priors[m]),
		[cm_trackers[m].n_accepted],
		[cm_trackers[m].n_tries],
		[input.threshold_schedule[1]],
        cm_trackers[m].n_accepted > 0,
		[cm_trackers[m].population],
		[cm_trackers[m].distances],
		[StatsBase.Weights(cm_trackers[m].weight_values ./ sum(cm_trackers[m].weight_values), 1.0)], 	# Normalise weights for each model
		input.parameter_priors[m],
		input.distance_simulation_input[m],
		1) for m in 1:input.M]

	return SimulatedModelSelectionTracker(input.M,
		input.n_particles,
		[input.threshold_schedule[1]],
		input.model_prior,
		smc_trackers,
		input.max_iter)
end

# Initialises the model selection run and runs the first (rejection) population
function initialise_modelselection(input::EmulatedModelSelectionInput{AF, CUD, ER, EPS, ET}) where {
	AF<:AbstractFloat, CUD<:ContinuousUnivariateDistribution,
	ER<:AbstractEmulatorRetraining, EPS<:AbstractEmulatedParticleSelection, ET<:AbstractEmulatorTraining}
	#
	# Check input sizes
	#
	if span(input.model_prior) != input.M
		throw(ArgumentError("There are $(input.M) models but the span of the model prior support is $(span(input.model_prior))"))
	end

	if length(input.parameter_priors) != input.M
		throw(ArgumentError("There are $(input.M) models but $(length(input.parameter_priors)) sets of parameter priors"))
	end

	#
	# Initialise arrays that will track rejection ABC run for each model - these
	# will be used to create CandidateModelTrackers after the rejection ABC run
	#
	cm_trackers = [CandidateModelTracker(length(input.parameter_priors[m]))
								for m in 1:input.M]

	#
	# Train the emulators
	#
	emulators = []
	for m in 1:input.M
		emulator = abc_train_emulator(input.parameter_priors[m], input.emulator_training_input[m])
		push!(emulators, emulator)
	end

    # Initialise variables to hold ABC results
    total_n_accepted = 0
	n_iterations = 1

	#
	# Emulate first population using rejection-ABC
	#
	while total_n_accepted < input.n_particles && n_iterations <= input.max_iter
		sampled_models = rand(input.model_prior, min(input.max_batch_size, input.n_particles-total_n_accepted))
		tries_this_it = [size(sampled_models[sampled_models.==m],1) for m=1:input.M]

		for m=1:input.M
			if tries_this_it[m] > 0

				parameters, weights = generate_parameters(input.parameter_priors[m], tries_this_it[m])
				distances, accepted_batch_idxs = abc_select_emulated_particles(emulators[m], parameters,
				 	input.threshold_schedule[1], input.emulated_particle_selection)
				cm_trackers[m].n_tries += tries_this_it[m]

				# If at least one particle is accepted
				if length(accepted_batch_idxs) > 0
					total_n_accepted += length(accepted_batch_idxs)
					update_candidatemodeltracker!(cm_trackers[m],
												  accepted_batch_idxs,
												  parameters,
												  distances,
												  weights)

				elseif length(accepted_batch_idxs) > tries_this_it[m]
					error("$(length(accepted_batch_idxs))) particles were accepted when model $m was only sampled $(tries_this_it[m]) times!")
				end
			end
		end

		n_iterations += 1
	end

	if n_iterations == input.max_iter+1
		@warn "Emulated model selection reached maximum number of iterations ($(input.max_iter)) without accepting $(input.n_particles) particles on the first population - consider trying more iterations."
	end

	@info join(["GpABC model selection emulation",
			"Completed $(n_iterations-1) iterations, accepting $total_n_accepted particles in total.",
			"Number of accepted parameters by model: ",
			join([string("Model ", m, ": ",  cm_trackers[m].n_accepted) for m in 1:input.M], "\t")
		], "\n")

	smc_trackers = [EmulatedABCSMCTracker{CUD, typeof(emulators[m]), ER, EPS}(
		length(input.parameter_priors[m]),
		[cm_trackers[m].n_accepted],
		[cm_trackers[m].n_tries],
		[input.threshold_schedule[1]],
		cm_trackers[m].n_accepted > 0,
		[cm_trackers[m].population],
		[cm_trackers[m].distances],
		[StatsBase.Weights(cm_trackers[m].weight_values ./ sum(cm_trackers[m].weight_values), 1.0)],
		input.parameter_priors[m],
		input.emulator_training_input[m],
		input.emulator_retraining,
		input.emulated_particle_selection,
		input.max_batch_size,
		input.max_iter,
		[emulators[m]]) for m = 1:input.M]

	return EmulatedModelSelectionTracker(
				input.M,
				input.n_particles,
				[input.threshold_schedule[1]],
				input.model_prior,
				smc_trackers,
				input.max_batch_size,
				input.max_iter)
end

# Perform a subsequent model selection iteration (based on ABC-SMC)
function iterate_modelselection!(tracker::SimulatedModelSelectionTracker,
	threshold::AbstractFloat)

	# Generate kernels for alive models only
	kernels = Dict{Int,Matrix{ContinuousUnivariateDistribution}}()
	for m in 1:tracker.M
		if tracker.smc_trackers[m].n_accepted[end] > 0
			kernels[m] = generate_kernels(tracker.smc_trackers[m].population[end],
										  tracker.smc_trackers[m].priors)
		end
	end

	# Initialise
	cm_trackers = [CandidateModelTracker(tracker.smc_trackers[m].n_params) for m in 1:tracker.M]
	total_n_accepted = 0
	n_iterations = 1
	parameters = zeros(0)
	distance = 0.0
	weight_value = 0.0

	while total_n_accepted < tracker.n_particles && n_iterations <= tracker.max_iter
		# Sample model from prior
		m = rand(tracker.model_prior)
		smc_tracker = tracker.smc_trackers[m]

		if smc_tracker.n_accepted[end] == 0
			continue
		end

		# Do ABC SMC for a single particle
		parameters, weight_value = generate_parameters(smc_tracker.priors, smc_tracker.weights[end], kernels[m])
		try
            distance = simulate_distance(parameters, smc_tracker.distance_simulation_input)
        catch e
            if isa(e, DimensionMismatch)
                # This prevents the whole code from failing if there is a problem
                # solving the differential equation(s). The exception is thrown by the
                # distance function
                @warn "The summarised simulated data does not have the same size as the summarised reference data. If this is not happening at every iteration it may be due to the behaviour of OrdinaryDiffEq::solve - please check for related warnings. Continuing to the next iteration."
                n_iterations += 1
                continue
            else
                throw(e)
            end
        end

        cm_trackers[m].n_tries += 1

		if distance[1] < threshold
			total_n_accepted += 1
			update_candidatemodeltracker!(cm_trackers[m], parameters, distance[1], weight_value)
		end

		n_iterations += 1

	end

	if n_iterations == tracker.max_iter+1
		@warn "Simulated model selection reached maximum number of iterations ($(tracker.max_iter)) on an SMC run - consider trying more iterations."
	end

	@info join(["GpABC model selection simulation",
			"Completed $(n_iterations-1) iterations, accepting $total_n_accepted particles in total.",
			"Number of accepted parameters by model: ",
			join([string("Model ", m, ": ",  cm_trackers[m].n_accepted) for m in 1:tracker.M], "\t")
		], "\n")

	update_modelselection_tracker!(tracker, cm_trackers, threshold)
end

function iterate_modelselection!(
	tracker::EmulatedModelSelectionTracker{CUD, ET, ER, EPS},
	threshold::AbstractFloat) where {CUD<:ContinuousUnivariateDistribution, ET,
	        ER<:AbstractEmulatorRetraining, EPS<:AbstractEmulatedParticleSelection}

	# Initialise
	cm_trackers = [CandidateModelTracker(length(smc_tracker.priors)) for smc_tracker in tracker.smc_trackers]

	# Generate kernels for alive models only
	kernels = Dict{Int,Matrix{ContinuousUnivariateDistribution}}()
	for m in 1:tracker.M
		if tracker.smc_trackers[m].n_accepted[end] > 0
			kernels[m] = generate_kernels(tracker.smc_trackers[m].population[end],
										  tracker.smc_trackers[m].priors)
		end
	end

	# Retrain emulators for alive models only
	retrained_emulators = Array{ET,1}(undef,tracker.M)
	for m in 1:tracker.M
		smc_tracker = tracker.smc_trackers[m]
		if smc_tracker.n_accepted[end] > 0
			particle_sampling_function(batch_size) = generate_parameters_no_weights(batch_size,
		        smc_tracker.population[end],
		        smc_tracker.weights[end],
		        kernels[m])
			retrained_emulators[m] = abc_retrain_emulator(smc_tracker.emulators[end],
				particle_sampling_function,
				threshold,
				smc_tracker.emulator_training_input,
				smc_tracker.emulator_retraining)
		else
			retrained_emulators[m] = smc_tracker.emulators[end]
		end
	end

    # Initialise
    total_n_accepted = 0
	n_iterations = 1

    # emulate
	while total_n_accepted < tracker.n_particles && n_iterations <= tracker.max_iter

		sampled_models = rand(tracker.model_prior, min(tracker.max_batch_size, tracker.n_particles-total_n_accepted))
		tries_this_it = [size(sampled_models[sampled_models.==m],1) for m=1:tracker.M]

		for m = 1:tracker.M
			# Skip a model that wasn't sampled in this iteration or had no accepted particles in the previous population
			if tries_this_it[m] == 0 || tracker.smc_trackers[m].n_accepted[end] == 0
				continue
			end
			smc_tracker = tracker.smc_trackers[m]
			parameters, weights = generate_parameters(tries_this_it[m],
	                                                 smc_tracker.priors,
	                                                 smc_tracker.weights[end],
	                                                 kernels[m])
			cm_trackers[m].n_tries += tries_this_it[m]

			distances, accepted_indices = abc_select_emulated_particles(retrained_emulators[m], parameters, threshold, smc_tracker.selection)

			# If at least one particle is accepted
			if length(accepted_indices) > 0
				total_n_accepted += length(accepted_indices)
				update_candidatemodeltracker!(cm_trackers[m],
												  accepted_indices,
												  parameters,
												  distances,
												  weights)

			elseif length(accepted_indices) > tries_this_it[m]
				error("$(length(accepted_indices)) particles were accepted when model $m was only sampled $(tries_this_it[m]) times!")
			end
		end

		n_iterations += 1
	end

	if n_iterations == tracker.max_iter+1
		@warn "Emulated model selection reached maximum number of iterations ($(tracker.max_iter)) on an SMC run - consider trying more iterations."
	end

	@info join(["GpABC model selection emulation",
			"Completed $(n_iterations-1) iterations, accepting $total_n_accepted particles in total.",
			"Number of accepted parameters by model: ",
			join([string("Model ", m, ": ",  cm_trackers[m].n_accepted) for m in 1:tracker.M], "\t")],
		"\n")

	update_modelselection_tracker!(tracker, cm_trackers, retrained_emulators, threshold)
end

# not exported
function build_modelselection_output(tracker::ModelSelectionTracker, successful_completion::Bool)
	return ModelSelectionOutput(
		tracker.M,
		[[tracker.smc_trackers[m].n_accepted[i] for m = 1:tracker.M] for i in 1:length(tracker.threshold_schedule)],
		tracker.threshold_schedule,
		[buildAbcSmcOutput(tracker.smc_trackers[m]) for m in 1:tracker.M],
		successful_completion
		)
end

# not exported
function all_models_dead(tracker::ModelSelectionTracker)
	return sum([tracker.smc_trackers[m].n_accepted[end] for m = 1:tracker.M]) == 0
end

# not exported
function all_but_one_models_dead(tracker::ModelSelectionTracker)
	return sum([tracker.smc_trackers[m].n_accepted[end] == 0 for m = 1:tracker.M]) == tracker.M-1
end

# not exported
function update_candidatemodeltracker!(
	cm_tracker::CandidateModelTracker,
	parameters::AbstractArray{Float64,2},
	distance::AbstractFloat,
	weight_value::AbstractFloat)

	cm_tracker.n_accepted += 1
	cm_tracker.population = vcat(cm_tracker.population, parameters)
	push!(cm_tracker.distances, distance)
	push!(cm_tracker.weight_values, weight_value)
end

function update_candidatemodeltracker!(
	cm_tracker::CandidateModelTracker,
	accepted_batch_idxs::AbstractArray{Int64,1},
	parameter_batch::AbstractArray{Float64,2},
	distance_batch::AbstractArray{Float64,1},
	weight_values_batch::AbstractArray{Float64,1})

	cm_tracker.n_accepted += length(accepted_batch_idxs)
	cm_tracker.population = vcat(cm_tracker.population, parameter_batch[accepted_batch_idxs,:])
	cm_tracker.distances = vcat(cm_tracker.distances, distance_batch)
	cm_tracker.weight_values = vcat(cm_tracker.weight_values, weight_values_batch[accepted_batch_idxs])
end

# not exported
function update_modelselection_tracker!(
	tracker::SimulatedModelSelectionTracker,
	cm_trackers::AbstractArray{CandidateModelTracker{AF},1},
	threshold::AF) where {AF<:AbstractFloat}
	push!(tracker.threshold_schedule, threshold)
	for m in 1:tracker.M
		push!(tracker.smc_trackers[m].n_accepted, cm_trackers[m].n_accepted)
		push!(tracker.smc_trackers[m].n_tries, cm_trackers[m].n_tries)
		push!(tracker.smc_trackers[m].threshold_schedule, threshold)
		push!(tracker.smc_trackers[m].population, cm_trackers[m].population)
		push!(tracker.smc_trackers[m].distances, cm_trackers[m].distances)
		push!(tracker.smc_trackers[m].weights, normalise(cm_trackers[m].weight_values))
	end
end

# not exported
function update_modelselection_tracker!(
	tracker::EmulatedModelSelectionTracker,
	cm_trackers::AbstractArray{CandidateModelTracker{AF},1},
	retrained_emulators::AbstractArray{GPModel,1},
	threshold::AF) where {AF<:AbstractFloat}

	push!(tracker.threshold_schedule, threshold)
	for m in 1:tracker.M
		push!(tracker.smc_trackers[m].n_accepted, cm_trackers[m].n_accepted)
		push!(tracker.smc_trackers[m].n_tries, cm_trackers[m].n_tries)
		push!(tracker.smc_trackers[m].threshold_schedule, threshold)
		push!(tracker.smc_trackers[m].population, cm_trackers[m].population)
		push!(tracker.smc_trackers[m].distances, cm_trackers[m].distances)
		push!(tracker.smc_trackers[m].weights, normalise(cm_trackers[m].weight_values))
		push!(tracker.smc_trackers[m].emulators, retrained_emulators[m])
	end
end


text:
vec_simple = @stochastic_model begin
    @stage 1 begin
        @parameters begin
            lb = [40.0, 20.0]
            c  = [100., 150.]
            A  = [1. 1.]
            b  = [120.]
        end
        @decision(model, x[i in 1:2] >= lb[i])
        @objective(model, Min, dot(c, x))
        @constraint(model, A * x - b in MOI.Nonpositives(1))
    end
    @stage 2 begin
        @parameters begin
            lb = [0., 0.]
            T  = [-60. 0.;
                  0.  -80.]
            W  = [6. 10.;
                  8. 5.]
        end
        @uncertain q₁ q₂ d₁ d₂
        ub = [d₁, d₂]
        q  = [q₁, q₂]
        @variable(model, lb[i] <= y[i in 1:2] <= ub[i])
        @objective(model, Max, dot(q, y))
        @constraint(model, T * x + W * y in MOI.Nonpositives(2))
    end
end

ξ₁ = Scenario(q₁ = 24.0, q₂ = 28.0, d₁ = 500.0, d₂ = 100.0, probability = 0.4)
ξ₂ = Scenario(q₁ = 28.0, q₂ = 32.0, d₁ = 300.0, d₂ = 300.0, probability = 0.6)

vec_simple_res = SPResult([46.67, 36.25], -855.83, -1518.75, 662.92, 286.92, -1445.92, -568.92)
push!(problems, (vec_simple, [ξ₁,ξ₂], vec_simple_res, "Simple"))


text:
abstract type AbstractMissile <: AbstractEnv end

#-----------------------------------------
# Single Point Mass Vehicle Kinematics
#-----------------------------------------
struct PointMassMissile <: AbstractMissile
    dim::Int  # @assert dim == 2 || dim == 3
    A_max::Float64
end

function State(env::PointMassMissile)
    @unpack dim = env
    return function (p = zeros(dim), v = zeros(dim))
        ComponentArray(p = p, v = v)
    end
end

function Dynamics!(env::PointMassMissile)
    @unpack A_max = env
    @Loggable function dynamics!(dx, x, params, t; u)
        @unpack p, v = x
        @log p
        @log v
        dx.p = v
        if norm(u) >= A_max
            u = normalize(u) * clamp(norm(u), 0, A_max)
        end
        dx.v = u
        @log u
    end
end

#-----------------------------------------
# 1 Pursuer on 1 Evador Engagement
#-----------------------------------------
struct PursuerEvadorMissile <: AbstractMissile
    pursuer::PointMassMissile
    evador::PointMassMissile
end

function State(env::PursuerEvadorMissile)
    @unpack pursuer, evador = env
    return function (x0_pursuer, x0_evador)
        ComponentArray(pursuer = x0_pursuer, evador = x0_evador)
    end
end

function Dynamics!(env::PursuerEvadorMissile, s_guidance::BPNG)
    @unpack pursuer, evador = env
    @unpack Bias = s_guidance
    @Loggable function dynamics!(dx, x, params, t; u_pursuer, u_evador)
        @onlylog r = norm(x.pursuer.p - x.evador.p)
        @onlylog σ_M = acos(clamp(dot(normalize(x.evador.p - x.pursuer.p), normalize(x.pursuer.v)), -1, 1))
        @onlylog a_M_bias, e_v̂_f, ω_bias, μ, Ω_μ = Bias(s_guidance, x, t)

        @nested_log :pursuer Dynamics!(pursuer)(dx.pursuer, x.pursuer, nothing, t; u = u_pursuer)
        @nested_log :evador Dynamics!(evador)(dx.evador, x.evador, nothing, t; u = u_evador)
    end
end

#-----------------------------------------
# Single Point Mass Vehicle Vertical Plane Dynamics
#-----------------------------------------
struct VerticalPlaneDynamicMissile <: AbstractMissile
    m::Float64
    S::Float64
    C_D_0::Float64
    K::Float64
    ρ₀::Float64
    H::Float64
    g::Float64
end

function State(env::VerticalPlaneDynamicMissile)
    return function (x₀, h₀, γ₀, V₀)
        ComponentArray(x = x₀, h = h₀, γ = γ₀, V = V₀)
    end
end

function Params(env::VerticalPlaneDynamicMissile)
    return function (R₀, c₀)
        ComponentArray(R = R₀, c = c₀)
    end
end

function Dynamics!(env::VerticalPlaneDynamicMissile)
    @unpack m, S, C_D_0, K, ρ₀, H, g = env
    @Loggable function dynamics!(ds, s, params, t; u)
        @unpack x, h, γ, V = s
        @unpack R, c = params
        @log x, h, γ, V, u
        @onlylog R, c

        ρ = ρ₀ * exp(-h / H)
        C_L = m * u / (1 / 2 * ρ * V^2 * S)
        C_D = C_D_0 + K * C_L^2

        ds.x = V * cos(γ)
        ds.h = V * sin(γ)
        ds.γ = u / V - g / V * cos(γ)
        ds.V = -C_D / C_L * u - g * sin(γ)
    end
end

#-----------------------------------------
# Single Rigid Body Vehicle Vertical Plane Dynamics
#-----------------------------------------
# [Disclaimer] Source of model parameters
# @inproceedings{Mracek_1997,
# 	author = {Mracek, Curtis P. and Cloutier, James R.},
# 	title = {{Full Envelope Missile Longitudinal Autopilot Design Using the State-Dependent Riccati Equation Method}},
# 	booktitle = {AIAA Guidance, Navigation, and Control Conference},
# 	year = {1997},
# 	month = {August},
# 	address = {New Orleans, LA, USA},
# 	doi = {10.2514/6.1997-3767}
# }

Base.@kwdef struct VerticalPlane_RigidBody_Missile <: AbstractMissile
    # s_aerocoeff::Vector{Float64} = [-0.3, 19.373, -31.023, -9.717, -1.948, 40.44, -64.015, 2.922, -11.803]
    # s_airframe::Vector{Float64} = [204.02, 247.439, 0.0409, 0.2286, 150.0, 0.7]
    # s_environment::Vector{Float64} = [9.8, 1.225, 8435.0, 1.4, 286.0, 288.15, 0.0065]

    s_aerocoeff::Vector{Float64} = ComponentArray(aₐ = -0.3, aₙ = 19.373, bₙ = -31.023, cₙ = -9.717, dₙ = -1.948, aₘ = 40.44, bₘ = -64.015, cₘ = 2.922, dₘ = -11.803)
    s_airframe::Vector{Float64} = ComponentArray(m = 204.02, I_yy = 247.439, S = 0.0409, d = 0.2286, ωₐ = 150.0, ζₐ = 0.7)
    s_environment::Vector{Float64} = ComponentArray(g = 9.8, ρ₀ = 1.225, H = 8435.0, γₐ = 1.4, Rₐ = 286.0, T₀ = 288.15, λ = 0.0065)
end

function State(env::VerticalPlane_RigidBody_Missile)
    return function (h₀, V₀, α₀, q₀, θ₀, δ₀, δ̇₀)
        ComponentArray(h = h₀, V = V₀, α = α₀, q = q₀, θ = θ₀, δ = δ₀, δ̇ = δ̇₀)
    end
end

function Dynamics!(env::VerticalPlane_RigidBody_Missile)
    @unpack s_aerocoeff, s_airframe, s_environment = env
    (aₐ, aₙ, bₙ, cₙ, dₙ, aₘ, bₘ, cₘ, dₘ)     = s_aerocoeff
    (m, I_yy, S, d, ωₐ, ζₐ)                 = s_airframe
    (g, ρ₀, H, γₐ, Rₐ, T₀, λ)               = s_environment

    @Loggable function dynamics!(ds, s, params, t; u)
        (h, V, α, q, θ, δ, δ̇) = s
        δ_c = u
        @log h, V, α, q, θ, δ, δ̇, δ_c
        # @onlylog

        ρ = ρ₀ * exp(-h / H)
        Vₛ = sqrt(γₐ * Rₐ * (T₀ - λ * h))
        M = V / Vₛ
        γ = θ - α
        Q = 0.5f0 * ρ * V^2

        C_A = aₐ
        C_N = aₙ * α^3 + bₙ * α * abs(α) + cₙ * (2.0f0 - M / 3.0f0) * α + dₙ * δ
        C_M = aₘ * α^3 + bₘ * α * abs(α) + cₘ * (-7.0f0 + 8.0f0 * M / 3.0f0) * α + dₘ * δ

        α̇ = Q * S / m / V * (C_N * cos(α) - C_A * sin(α)) + g / V * cos(γ) + q
        a_z = V * (q - α̇)

        ds.h = V * sin(γ)
        ds.V = Q * S / m * (C_N * sin(α) + C_A * cos(α)) - g * sin(γ)
        ds.α = α̇
        ds.q = Q * S * d / I_yy * C_M
        ds.θ = q
        ds.δ = δ̇
        ds.δ̇ = -ωₐ^2 * (δ - δ_c) - 2.0f0 * ζₐ * ωₐ * δ̇
    end
end


text:
const UniFinArr = UnivariateFiniteArray

Base.size(u::UniFinArr, args...) =
    size(first(values(u.prob_given_ref)), args...)

function Base.getindex(u::UniFinArr{<:Any,<:Any,R,P,N},
                       i::Integer...) where {R,P,N}
    prob_given_ref = LittleDict{R,P}()
    for ref in keys(u.prob_given_ref)
        prob_given_ref[ref] = getindex(u.prob_given_ref[ref], i...)
    end
    return UnivariateFinite(u.scitype, u.decoder, prob_given_ref)
end

function Base.getindex(u::UniFinArr{<:Any,<:Any,R,P,N},
                       I...) where {R,P,N}
    prob_given_ref = LittleDict{R,Array{P,N}}()
    for ref in keys(u.prob_given_ref)
        prob_given_ref[ref] = getindex(u.prob_given_ref[ref], I...)
    end
    return UnivariateFiniteArray(u.scitype, u.decoder, prob_given_ref)
end

function Base.setindex!(u::UniFinArr{S,V,R,P,N},
                        v::UnivariateFinite{S,V,R,P},
                        i::Integer...) where {S,V,R,P,N}
    for ref in keys(u.prob_given_ref)
       setindex!(u.prob_given_ref[ref], v.prob_given_ref[ref], i...)
    end
    return u
end

_err_incompatible_levels() = throw(DomainError(
    "Cannot concatenate `UnivariateFiniteArray`s with "*
    "different categorical levels (classes), "*
    "or whose levels, when ordered, are not  "*
    "consistently ordered. "))

# terminology:

# "classes"  - full pool of `CategoricalElement`s, even "unseen" ones (those
#             missing from support)
# "levels"   - same thing but in raw form (eg, `Symbol`s) aka "labels"
# "suppport" - those classes with a corresponding probability (the ones
#              named at time of construction of the `UnivariateFiniteArray`)

function Base.cat(us::UnivariateFiniteArray{S,V,R,P,N}...;
                      dims::Integer) where {S,V,R,P,N}

    isempty(us) && return []

    # build combined raw_support and check compatibility of levels:
    u1 = first(us)
    ordered = isordered(classes(u1))
    support_with_duplicates = Dist.support(u1)
    _classes = classes(u1)
    for i in 2:length(us)
        isordered(us[i]) == ordered || _err_incompatible_levels()
        if ordered
            classes(us[i]) ==
                _classes|| _err_incompatible_levels()
        else
            Set(classes(us[i])) ==
                Set(_classes) || _err_incompatible_levels()
        end
        support_with_duplicates =
            vcat(support_with_duplicates, Dist.support(us[i]))
    end
    _support = unique(support_with_duplicates) # no-longer categorical!

    # build the combined `prob_given_class` dictionary:
    pairs = (class => cat((pdf.(u, class) for u in us)..., dims=dims)
             for class in _support)
    prob_given_class = Dict(pairs)

    return UnivariateFinite(prob_given_class, pool=_classes)
end

Base.vcat(us::UnivariateFiniteArray...) = cat(us..., dims=1)
Base.hcat(us::UnivariateFiniteArray...) = cat(us..., dims=2)

## CONVENIENCE METHOD pdf(array_of_univariate_finite, labels)

# next bit is not specific to `UnivariateFiniteArray` but is for any
# abstract array with eltype `UnivariateFinite`.

# this is type piracy that has been adopted only after much
# agonizing over alternatives. Note that pdf.(u, labels) must
# necessarily have a different meaning (and only makes sense if u and
# labels have the same length or labels is a scalar)

function Distributions.pdf(
    u::AbstractArray{UnivariateFinite{S,V,R,P},N},
    C::AbstractVector{<:Union{V, CategoricalValue{V,R}}}) where {S,V,R,P,N}

    ret = Array{P,N+1}(undef, size(u)..., length(C))
    for i in eachindex(C)
        ret[fill(:,N)...,i] = broadcast(pdf, u, C[i])
    end
    return ret
end

## PERFORMANT BROADCASTING OF pdf

# u - a UnivariateFiniteArray
# cv - a CategoricalValue
# v - a vector of CategoricalArrays

# pdf.(u, cv)
function Base.Broadcast.broadcasted(
    ::typeof(pdf),
    u::UniFinArr{S,V,R,P,N},
    cv::CategoricalValue) where {S,V,R,P,N}

    cv in classes(u) || _err_missing_class(cv)

    return get(u.prob_given_ref, int(cv), zeros(P, size(u)))
end

# pdf.(u, v)
function Base.Broadcast.broadcasted(
    ::typeof(pdf),
    u::UniFinArr{S,V,R,P,N},
    v::AbstractArray{<:CategoricalValue{V,R},N}) where {S,V,R,P,N}

    length(u) == length(v) ||throw(DimensionMismatch(
        "Arrays could not be broadcast to a common size; "*
        "got a dimension with lengths $(length(u)) and $(length(v))"))
    for cv in v
        cv in classes(u) || _err_missing_class(c)
    end

    # will use linear indexing:
    v_flat = ((v[i], i) for i in 1:length(v))

    getter((cv, i), dtype) = _getindex(get(u.prob_given_ref, int(cv), nothing), i, dtype)

    ret_flat = getter.(v_flat, P)
    return reshape(ret_flat, size(u))
end

#dummy function
# returns `x[i]` for `Array` inputs `x`
# For non-Array inputs returns `zero(dtype)`
#This avoids using an if statement
_getindex(x::Array,i, dtype)=x[i]
_getindex(::Nothing, i, dtype) = zero(dtype)

# pdf.(u, raw) where raw is scalar or vec
function Base.Broadcast.broadcasted(
    ::typeof(pdf),
    u::UnivariateFiniteArray{S,V,R,P,N},
    raw::Union{V,AbstractArray{V,N}}) where {S,V,R,P,N}

    cat = transform(classes(u), raw)
    return broadcast(pdf, u, cat)
end

## PERFORMANT BROADCASTING OF mode:

function Base.Broadcast.broadcasted(::typeof(mode),
                                    u::UniFinArr{S,V,R,P,N}) where {S,V,R,P,N}
    dic = u.prob_given_ref

    # using linear indexing:
    mode_flat = map(1:length(u)) do i
        max_prob = maximum(dic[ref][i] for ref in keys(dic))
        m = zero(R)
        for ref in keys(dic)
            if dic[ref][i] == max_prob
                m = ref
                break
            end
        end
        return u.decoder(m)
    end
    return reshape(mode_flat, size(u))
end


text:
module EMLDoublePipeIncr
	using DanaTypes
	using DotPlusInheritance
	using Reexport
	@reexport using ...streams.EMLstreams
	import EMLtypes.length
	include("DoublePipeIncr/Properties_Average.jl")
	include("DoublePipeIncr/Properties_In_Out.jl")
	include("DoublePipeIncr/Properties_Wall.jl")
	include("DoublePipeIncr/Physical_Properties.jl")
	include("DoublePipeIncr/Details_Main.jl")
	include("DoublePipeIncr/DoublePipe_HeatTransfer.jl")
	include("DoublePipeIncr/DoublePipe_PressureDrop.jl")
	include("DoublePipeIncr/Main_DoublePipe.jl")
	include("DoublePipeIncr/DoublePipeIncr.jl")
end


text:
using Documenter
using FINUFFT

makedocs(sitename="FINUFFT.jl")

deploydocs(
    repo = "github.com/ludvigak/FINUFFT.jl.git",
)