text:
# Utilities to deal with fetching/serving actual Pkg resources

const uuid_re = raw"[0-9a-z]{8}-[0-9a-z]{4}-[0-9a-z]{4}-[0-9a-z]{4}-[0-9a-z]{12}(?-i)"
const hash_re = raw"[0-9a-f]{40}"
const registry_re = Regex("^/registry/($uuid_re)/($hash_re)\$")
const resource_re = Regex("""
    ^/registry/$uuid_re/$hash_re\$
  | ^/package/$uuid_re/$hash_re\$
  | ^/artifact/$hash_re\$
""", "x")

"""
    get_registries(server)

Interrogate a storage server for a list of registries, match the response against the
registries we are paying attention to, return dict mapping from registry UUID to its
latest treehash.
"""
function get_registries(server::AbstractString)
    regs = Dict{String,String}()
    response = HTTP.get("$server/registries", status_exception = false)
    if response.status != 200
        @error("Failure to fetch /registries", server, response.status)
        return regs
    end
    for line in eachline(IOBuffer(response.body))
        m = match(registry_re, line)
        if m !== nothing
            uuid, hash = m.captures
            uuid in keys(config.registries) || continue
            regs[uuid] = hash
        else
            @error("invalid response", server, resource="registries", line)
        end
    end
    return regs
end

function prune_empty_parents(child, root)
    if !isdir(child)
        return
    end

    root = rstrip(abspath(root), '/')
    child = rstrip(abspath(child), '/')
    last_child = ""

    while child != root && child != last_child
        if !isempty(readdir(child))
            break
        end

        last_child = child
        child = dirname(child)
        rm(last_child; recursive=true)
    end
    return
end

"""
    write_atomic_lru(f::Function, resource::AbstractString)

Performs an atomic filesystem write by writing out to a file on the same
filesystem as the given `path`, then `move()`'ing the file to its eventual
destination.  Requires write access to the file and the containing folder.
Currently stages changes at "<path>.tmp.<randstring>".  If the return value
of `f()` is `false` or an exception is raised, the write will be aborted.

Also tracks the file with the global LRU cache as configured in `config.cache`.
"""
function write_atomic_lru(f::Function, resource::AbstractString)
    temp_file = temp_resource_filepath(resource)
    try
        # First, write a temp file into the `temp` directory:
        mkpath(dirname(temp_file))
        retval = open(temp_file, "w") do io
            f(temp_file, io)
        end
        if retval !== false
            # Calculate size of the file, notify the cache that we're adding
            # a file of that size, so it may need to shrink the cache:
            new_path = add!(config.cache, resource[2:end], filesize(temp_file))
            mv(temp_file, new_path; force=true)
        end
        return retval
    catch e
        rethrow(e)
    finally
        rm(temp_file; force=true)
        prune_empty_parents(temp_file, joinpath(config.root, "temp"))
    end
end

function resource_filepath(resource::AbstractString)
    # We strip off the leading `/` to pass this into the filecache
    return filepath(config.cache, resource[2:end])
end

function temp_resource_filepath(resource::AbstractString)
    return joinpath(config.root, "temp", string(resource[2:end], ".inprogress"))
end

"""
    url_exists(url)

Send a `HEAD` request to the specified URL, returns `true` if the response is HTTP 200.
"""
function url_exists(url::AbstractString)
    response = HTTP.request("HEAD", url, status_exception = false)
    response.status == 200
end

"""
    verify_registry_hash(uuid, hash)

Verify that the origin git repository knows about the given registry tree hash.
"""
function verify_registry_hash(uuid::AbstractString, hash::AbstractString)
    url = Pkg.Operations.get_archive_url_for_version(config.registries[uuid].upstream_url, hash)
    return url === nothing || url_exists(url)
end

function update_registries()
    # collect current registry hashes from servers
    regs = Dict(uuid => Dict{String,Vector{String}}() for uuid in keys(config.registries))
    servers = Dict(uuid => Vector{String}() for uuid in keys(config.registries))
    for server in config.storage_servers
        for (uuid, hash) in get_registries(server)
            push!(get!(regs[uuid], hash, String[]), server)
            push!(servers[uuid], server)
        end
    end
    # for each hash check what other servers know about it
    changed = false
    for (uuid, hash_info) in regs
        isempty(hash_info) && continue # keep serving what we're serving
        for (hash, hash_servers) in hash_info
            for server in servers[uuid]
                server in hash_servers && continue
                url_exists("$server/registry/$uuid/$hash") || continue
                push!(hash_servers, server)
            end
        end
        # Sort hashes by number of servers that know about it, to
        # serve the "newest" hashes first.
        hashes = collect(keys(hash_info))
        sort!(hashes, by = hash -> length(hash_info[hash]))
        for hash in hashes
            resource = "/registry/$uuid/$hash"

            # If this hash is not already on the filesystem, we might need to fetch it!
            hash_servers = sort!(hash_info[hash])
            if !isfile(resource_filepath(resource))
                # check that the origin repo knows about this hash.  This prevents a
                # rogue storage server from serving malicious registry tarballs.
                if !verify_registry_hash(uuid, hash)
                    @debug("rejecting untrusted registry hash", uuid, hash)
                    continue
                end

                # Try to fetch from our servers, but if something goes wrong, continue
                dl_state = fetch_resource(resource, servers=hash_servers)
                if dl_state === nothing
                    continue
                end
                # If nothing went wrong, pause and wait for the download to finish.
                wait(dl_state.dl_task)
            end

            if config.registries[uuid].latest_hash != hash
                @info("new current registry hash", uuid, hash, hash_servers)
                changed = true
            end

            # we've got a new registry hash to serve
            config.registries[uuid].latest_hash = hash
            break
        end
    end

    # write new registry info to file
    registries_path = joinpath(config.root, "static", "registries")
    if changed || !isfile(registries_path)
        new_registries = joinpath(config.root, "temp", "registries.tmp." * randstring())
        mkpath(dirname(new_registries))
        open(new_registries, "w") do io
            for uuid in sort!(collect(keys(config.registries)))
                println(io, "/registry/$(uuid)/$(config.registries[uuid].latest_hash)")
            end
        end
        mkpath(dirname(registries_path))
        mv(new_registries, registries_path; force=true)
    end
    return changed
end

"""
    NUM_FETCH_STATES

We shard our fetch state across `NUM_FETCH_STATES` different buckets, so that we can
theoretically have `M` different requests going at once for different resources, but the
same resource gets blocked on multiple simultaneous cache misses.  In practice we'll have
many accidental collisions before truly saturating, but that's fine since the happy path
(a cache hit) doesn't even need to enter this block, so unnecessary blocking should
hopefully be quite rare.  In general, we can caluclate the probability of accidental
blocking by taking the average number of simultaneous cache misses `M` and the number of
buckets `N` and plugging them into the birthday problem formula:

    prob_blocked(N, M) = 1 - factorial(big(N))/(big(N)^M * factorial(big(N - M)))

Therefore, making the blind assertion that we, on average, suffer two simultaneous
cache misses, a bucket number of 128 yields an accidental blocking chance of 0.7%.
That sounds pretty good to me, so 128 it is.
"""
const NUM_FETCH_STATES = 128

struct DownloadState
    # The actual resource (from which we can get the path via `resource_filepath()` and `temp_resource_filepath()`
    resource::String
    # Full length of the resource (obtained via the HEAD request to the storage server)
    content_length::Int
    # Task that is doing the downloading
    dl_task::Task
end

struct FetchState
    lock::ReentrantLock
    failed::Set{String}
    inprogress::Dict{String,DownloadState}

    FetchState() = new(ReentrantLock(), Set{String}(), Dict{String,DownloadState}())
end
const FETCH_STATES = [FetchState() for idx in 1:NUM_FETCH_STATES]

function with_fetch_state(f::Function, resource::AbstractString)
    state_idx = mod1(hash(resource), NUM_FETCH_STATES)
    state = FETCH_STATES[state_idx]
    lock(state.lock) do
        return f(state)
    end
end

"""
    select_server(resource, servers)

Given a resource and a list of storage servers, return the storage server that responds
most quickly to a HEAD request for that storage server, as well as the HEAD response
itself so that metadata such as the content-length of the resource can be inspected.
"""
function select_server(resource::AbstractString, servers::Vector{<:AbstractString}; timeout = 5, retries = 2)
    function head_req(server, resource)
        @try_printerror begin
            response = HTTP.head(
                string(server, resource);
                status_exception = false,
                timeout = timeout,
                retries = retries,
            )
            return server, response
        end
    end
    # Launch one task per server, performing a HEAD request
    tasks = [@spawn head_req(server, resource) for server in servers]

    # Wait for the first Task that gives us an HTTP 200 OK, returning that server.
    # If none have it, we return `nothing`. :(
    while !isempty(tasks)
        next_task = wait_first(tasks...)
        deleteat!(tasks, findfirst(t -> t == next_task, tasks))
        @try_printerror begin
            server, http_response = fetch(next_task)
            if http_response.status == 200
                return server, http_response
            end
        end
    end
    return nothing, nothing
end

function content_length(resp::HTTP.Messages.Response)
    for h in resp.headers
        if lowercase(h[1]) == "content-length"
            return parse(Int, h[2])
        end
    end
    return nothing
end

"""
    fetch_resource(resource::AbstractString; servers=config.storage_servers)

Reaches out to the list of storage servers, requesting `resource` from the server that
responds with an HTTP 200 OK the fastest.  Downloads on an asynchronous task stored as a
part of the `DownloadState` structure that this function returns.  Failures to download
are recorded and future downloads of that same resource will be skipped, until
`forget_failures()` is called.  The `DownloadState` object contains within it enough
information to still serve a resource as it is being downloaded in the background task.
"""
function fetch_resource(resource::AbstractString; servers::Vector{String}=config.storage_servers)
    if isempty(servers)
        @error("fetch called with no servers", resource)
        error("fetch called with no servers")
    end

    # with_fetch_state() will wait for a lock
    with_fetch_state(resource) do state
        # check if this has failed to download recently
        if resource in state.failed
            @debug("skipping recently failed download", resource)
            return nothing
        end

        # check if this resource is being downloaded already
        if resource in keys(state.inprogress)
            @debug("detected in-progress download", resource)
            return state.inprogress[resource]
        end

        # If not, let's figure out which storage server we're going to download from:
        server, response = select_server(resource, servers)
        if response === nothing
            @debug("no upstream server", resource, servers)
            return nothing
        end

        # Launch download process in a separate task:
        dl_task = @async begin
            success = download(server, resource, content_length(response))
            lock(state.lock) do
                if success
                    global fetch_hits += 1
                else
                    # If the download failed, wait a bit before retrying
                    push!(state.failed, resource)
                end
                # Once downloading is done, remove this resource from the list of inprogress downloads
                delete!(state.inprogress, resource)
            end
        end

        # Generate a DownloadState, map that to this resource
        state.inprogress[resource] = DownloadState(resource, content_length(response), dl_task)
        return state.inprogress[resource]
    end
end

function forget_failures()
    for idx in 1:NUM_FETCH_STATES
        lock(FETCH_STATES[idx].lock) do
            empty!(FETCH_STATES[idx].failed)
        end
    end
end

"""
    stream_file(io_in::IO, start_byte::Int, length::Int, dl_task::Task, io_out::IO)

Streams from a file `io_in` to an output `IO` object `io_out`.  Reads from `io_in` until
the requested `length` bytes are read, or the download task `dl_task` is finished and it
is certain no more bytes are going to be stored within the file.  Returns the total
number of bytes read from `io_in` and output into `io_out`.
"""
function stream_file(io_in::IO, start_byte::Int, length::Int, dl_task::Task, io_out::IO,
                     buffer::Vector{UInt8} = Vector{UInt8}(undef, 2*1024*1024))
    # Because this file may be an incomplete stream, we need to wait until our start byte
    # is  ready, so we `sleep` in a loop while we attempt to `seek()`
    seek(io_in, start_byte)
    while position(io_in) != start_byte
        sleep(0.01)
        seek(io_in, start_byte)
    end

    transmitted = 0
    while transmitted < length
        # See JuliaLang/julia#36300, can be optimized later to only read r bytes
        # r = min(length(buffer), content_length - t)
        # n = readbytes!(io, buffer, r)
        n = readbytes!(io_in, buffer)

        # If we got nothing, either the file is prematurely truncated, or we're still downloading it somewhere.
        if n == 0
            # If we're still streaming, just sleep for a bit before trying again.
            if dl_task.state != :done
                sleep(0.001)
            else
                # Otherwise, break out, something went wrong. :/
                break
            end
        else
            transmitted += write(io_out, view(buffer, 1:min(n, length - transmitted)))
        end
    end
    return transmitted
end

function download(server::AbstractString, resource::AbstractString, content_length::Int)
    @info "downloading resource" server=server resource=resource
    hash = basename(resource)

    write_atomic_lru(resource) do temp_file, file_io
        # We've got a moderately complex flow of data here.  We manage it all using
        # asynchronous tasks and processes.  In summary, we need to download the given
        # resource, decompress it, and determine its tree hash.  Complicating this is the
        # fact that the tree hash could be an older tree hash that skipped empty
        # directories, so we need to tree hash it twice.  And we want to do this in a
        # streaming fashion, so as to minimize latency for the user.
        #
        # We therefore:
        #   - Stream HTTP request into `file_io`, as fast as it can go.  This is to
        #     ensure that when we cache miss, we can `cat` out that file to the user
        #     as quickly as possible.
        #   - `cat` `file_io` into `gzip_proc`, to decompress it for live tree hashing,
        #     while `serve_file()` can stream `file_io` out to the client.
        #   - tee `gzip_proc` into two separate `tar_task_*` tasks, to compute the skip-
        #     empty and non-skip-empty hashes simultaneously.
        #   - If either of the hashes match, we keep the file.  If it was originally
        #     identified by the `skip_empty` version of the hash, store it under both.

        # Backing buffers for `tee` nodes
        tar_skip_buffio = BufferStream(16*1024*1024)
        tar_noskip_buffio = BufferStream(16*1024*1024)

        # Create gzip process to decompress for us, using `gzip()` from `Gzip_jll`
        gzip_proc = gzip(gz -> open(`$gz -d`, read=true, write=true))

        # Create tee node splitting gzip -> (skip, noskip)
        tar_tee_task = tee_task(gzip_proc.out, tar_skip_buffio, tar_noskip_buffio)

        # Create two tasks to read in the gzip output and do in-tar tree hashing.
        tar_skip_task = @async Tar.tree_hash(tar_skip_buffio; skip_empty=true)
        tar_noskip_task = @async Tar.tree_hash(tar_noskip_buffio; skip_empty=false)

        # Initiate the HTTP request, and start the data flowing into the file
        http_task = @async HTTP.get(server * resource,
            status_exception = false,
            response_stream = file_io,
        )

        # Read data back out from that file into the decompressor
        file_read_task = @async begin
            open(temp_file, read=true) do read_file_io
                written = stream_file(read_file_io, 0, content_length, http_task, gzip_proc.in)
                global payload_bytes_received += written
            end
            close(gzip_proc.in)
        end

        # Raise warnings about bad HTTP response codes
        response = fetch(http_task)
        if response.status != 200
            @warn "response status $(response.status)"
            return false
        end

        # Wait for the file read and tee tasks to finish
        wait(file_read_task)
        wait(tar_tee_task)

        # Fetch the result of the tarball hash check
        calc_skip_hash = fetch(tar_skip_task)
        calc_noskip_hash = fetch(tar_noskip_task)

        # If nothing matches, freak out.
        if hash != calc_skip_hash && hash != calc_noskip_hash
            @warn "resource hash mismatch" server resource hash calc_skip_hash calc_noskip_hash
            return false
        end

        # If calc_skip_hash matches, then store the file under the true hash as well.
        if hash != calc_noskip_hash && hash == calc_skip_hash
            @warn "archaic skip hash detected" resource hash calc_noskip_hash

            noskip_resource = joinpath(basename(resource), calc_noskip_hash)
            write_atomic_lru(noskip_resource) do noskip_temp_file, noskip_file_io
                close(noskip_file_io)
                rm(noskip_temp_file)
                cp(temp_file, noskip_temp_file)
            end
        end

        return true
    end
end

function serve_file(http::HTTP.Stream,
                    io::IO,
                    content_type::AbstractString,
                    content_encoding::AbstractString;
                    buffer::Vector{UInt8} = Vector{UInt8}(undef, 2*1024*1024),
                    content_length = filesize(io),
                    dl_task::Task = @async(nothing))
    # Initialize range parameters
    startbyte, stopbyte = 0, content_length-1

    # Support single byte ranges
    range_string = HTTP.header(http, "Range")
    if !isempty(range_string)
        # Look for the following patterns: "bytes=a-b", "bytes=a-", and "bytes=-b".
        # Empty captures will fail tryparse and then be replaced with the correct endpoints.
        m = match(r"^\s*bytes\s*=\s*(\d*)\s*-\s*(\d*)\s*$", range_string)
        if m === nothing
            @error "unsupported HTTP Range request, ignoring" range_string
        else
            requested_startbyte = max(something(tryparse(Int, m[1]), startbyte), startbyte)
            requested_stopbyte = min(something(tryparse(Int, m[2]), stopbyte), stopbyte)
            if requested_stopbyte - requested_startbyte >= 0
                startbyte = requested_startbyte
                stopbyte = requested_stopbyte
                HTTP.setstatus(http, 206) # Partial Content
                HTTP.setheader(http, "Content-Range" => "bytes $(startbyte)-$(stopbyte)/$(content_length)")
                content_length = stopbyte - startbyte + 1
            end
        end
    end

    HTTP.setheader(http, "Content-Length" => string(content_length))
    HTTP.setheader(http, "Accept-Ranges" => "bytes")
    HTTP.setheader(http, "Content-Type" => content_type)
    if content_encoding != "identity"
        HTTP.setheader(http, "Content-Encoding" => content_encoding)
    end
    startwrite(http)

    # Account this hit
    global total_hits += 1

    # Only write content if the method is `GET` (e.g. `HEAD` requests get the above headers)
    if http.message.method == "GET"
        transmitted = stream_file(io, startbyte, content_length, dl_task, http, buffer)
        global payload_bytes_transmitted += transmitted
        if transmitted != content_length
            @error("file size mismatch", path, content_length, actual=transmitted)
        end
    end
end


text:
module Truss

using LinearAlgebra
using Plots
using GraphRecipes

export define, solve, show

"""
        truss_object

This object holds all the definitions (material properties, section properties, geometry) and the solution (displacements, reactions, internal forces) for a system of truss elements.

num_dof_per_node::Int64
members::Array{NTuple{4,Int64},1}
node_geometry::Array{Float64,2}
section_properties::Array{Float64,1}
material_properties::Array{Float64,1}
supports::Array{Int64,1}
external_forces::Array{Float64,1}
L::Array{Float64,1}
E::Array{Float64,1}
A::Array{Float64,1}
θ::Array{Float64,1}
K::Array{Float64,2}
Kff::Array{Float64,2}
Ksf::Array{Float64,2}
Kfs::Array{Float64,2}
Kss::Array{Float64,2}
F::Array{Float64,1}
Fs::Array{Float64,1}
Ff::Array{Float64,1}
s::Array{Int64,1}
f::Array{Int64,1}
T::Array{Array{Float64,2},1}
k_element_local::Array{Array{Float64,2},1}
k_element_global::Array{Array{Float64,2},1}
u::Array{Float64,1}
uf::Array{Float64,1}
f_element::Array{Array{Float64,1},2}

"""

mutable struct truss_object

    num_dof_per_node::Int64
    members::Array{NTuple{4,Int64},1}
    node_geometry::Array{Float64,2}
    section_properties::Array{Float64,1}
    material_properties::Array{Float64,1}
    supports::Array{Int64,1}
    external_forces::Array{Float64,1}
    L::Array{Float64,1}
    E::Array{Float64,1}
    A::Array{Float64,1}
    θ::Array{Float64,1}
    K::Array{Float64,2}
    Kff::Array{Float64,2}
    Ksf::Array{Float64,2}
    Kfs::Array{Float64,2}
    Kss::Array{Float64,2}
    F::Array{Float64,1}
    Fs::Array{Float64,1}
    Ff::Array{Float64,1}
    s::Array{Int64,1}
    f::Array{Int64,1}
    T::Array{Array{Float64,2},1}
    k_element_local::Array{Array{Float64,2},1}
    k_element_global::Array{Array{Float64,2},1}
    u::Array{Float64,1}
    uf::Array{Float64,1}
    f_element::Array{Array{Float64,1},2}

    #This allows for the construction of a new object.
    truss_object() = new()

end

"""
        calculate_truss_element_lengths(members, node_geometry)

Calculate the truss member lengths in a structural system.

`members` is a NTuple{4, Int64} containing node i, node j, A, E assignments
`node_geometry` is array with the global x-y coordinates

"""

function calculate_truss_element_lengths(members, node_geometry)

    #Initialize array to hold truss lengths.
    truss_lengths=zeros(Float64, length(members))

    #Loop over each truss element.
    for i=1:length(members)

        #Define truss node numbers.
        node_i=members[i][1]
        node_j=members[i][2]

        #Define truss node coordinates.
        node_i_xy=node_geometry[node_i,:]
        node_j_xy=node_geometry[node_j,:]

        #Calculate the truss element length.
        truss_lengths[i] =norm(node_i_xy - node_j_xy)

    end

    return truss_lengths

end

"""
    define_local_truss_element_stiffness_matrix(E, A, L)

Define a truss element stiffness matrix in its local coordinate system.

`E` is elastic modulus.
`A` is the cross-sectional area.
`L` is truss length.

There are 4 degrees of freedom for the truss element.

      2                 4
    1>^--------------->3^

"""

function define_local_truss_element_stiffness_matrix(E, A, L)

	k_truss_element_local = E * A/ L * [1  0 -1 0
		                  0  0 0 0
		                  -1 0  1 0
	                      0  0  0 0]

	return k_truss_element_local

end

"""
        define_truss_element_orientations(members, node_geometry)

Define the truss element orientation angle.

`members` is a NTuple{4, Int64} containing node i, node j, A, E assignments
`node_geometry` is array with the global x-y coordinates

A positive angle is defined counterclockwise from the horizon, rotating about node i.

"""

function define_truss_element_orientations(members, node_geometry)

    # Initialize the orientation array for all truss elements in a structural system.
    θ = zeros(Float64, length(members))

    # Loop over the truss elements.
    for i=1:length(members)

        point_i = node_geometry[members[i][1], :]
        point_j = node_geometry[members[i][2], :]

        element_vector = point_j - point_i

        #Calculate the orientation angle.
        #https://en.wikipedia.org/wiki/Atan2
        θ[i] = atan(element_vector[2], element_vector[1])

    end

    return θ

end

"""
        vector_rotation_operator(θ)

Define an operator matrix that converts truss element displacements to a new, rotated coordinate system.

`θ` is the orientation angle of a truss element.

A positive angle is defined counterclockwise from the horizon, rotating about node i.

"""

function vector_rotation_operator(θ)

	T = [cos(θ) sin(θ) 0 0
		 -sin(θ) cos(θ) 0 0
		 0 0 cos(θ) sin(θ)
		 0 0  -sin(θ) cos(θ)]

	return T

end

"""
    assign_member_properties(members, member_property, property_order, property_type)

Define an array containing the member property values for each element in a structural system.   For example, define all the cross-sectional areas for all the members in a truss.

`members` is a NTuple{4, Int64} containing node i, node j, A, E assignments.
`member_property` is an array containing all the available properties that a member can take on.
`property order` is the index in the `members` NTuple where the specific property is assigned.

"""

function assign_member_properties(members, member_property, property_order, property_type)

    #Initialize the property array
    property = zeros(Float64, length(members))

    #Loop over all the members.
    for i=1:length(members)

        #Assign a specific property (area, elastic modulus) to each element in a structural system.
        property[i] = member_property[members[i][property_order]][property_type]

    end

    return property

end

"""
    calculate_global_element_stiffness_matrix(k_element_local, θ)

Transform a local element stiffness matrix to global coordinates.

`k_element_local` is the element stiffness matrix in its local coordinate system.
`θ` is the orientation angle of the element.

"""

function calculate_global_element_stiffness_matrix(k_element_local, T)

    #Transform the local element stiffness matrix to global coordinates.
    k_element_global = T' * k_element_local * T

    return k_element_global

end

"""
    assemble_global_stiffness_matrix(node_geometry, members, k_element_global, num_dof_per_node)

Assemble all the element stiffness matrices into the global system stiffness matrix.

`node_geometry` is an array with the global node x-y coordinates.
`members` is a NTuple{4, Int64} containing node i, node j, and property assignments.
`k_element_global` are all the element stiffness matrices in a 3D array.
`num_dof_per_node` defines the number of degrees of freedom per element node.  For a truss, this is 2, for a beam it is 3, ...

"""

function assemble_global_stiffness_matrix(node_geometry, members, k_element_global, num_dof_per_node)

    num_nodes = size(node_geometry)[1]

    k_system_global = zeros(Float64, num_nodes * num_dof_per_node, num_nodes * num_dof_per_node)

    for i=1:length(members)

        node_i = members[i][1]
        node_j = members[i][2]

        node_i_dof = [1;2] .+ (node_i - 1) * num_dof_per_node
        node_j_dof = [1;2] .+ (node_j - 1) * num_dof_per_node

        global_dof = [node_i_dof; node_j_dof]

        k_system_global[global_dof, global_dof] = k_system_global[global_dof, global_dof] + k_element_global[i]

    end

    return k_system_global

end

"""
    define(members, section_properties, material_properties, node_geometry, supports, external_forces)

Define and partition the global system stiffness matrix and external force vector.

`members` is a NTuple{4, Int64} containing node i, node j, and property assignments.
'section_properties' defines a library of section properties.
'material_properties' defines a library of material properties.
`node_geometry` is an array with the global node x-y coordinates.
`supports` defines the degrees of freedom that are free (0) and fixed(1).
`external_forces` is an array of the external forces applied to the nodes, in global coordinates.

Partition the stiffness matrix into Kff (free free), Kfs (free supported), Ksf (supported free), and Kss (supported supported).  Partition the external forces at the free dof into Ff. Also return `f` array of free global degrees of freedom, `s` fixed global degrees of freedom, `T` all the rotation matrices for each truss element, and the `k_element_local` and `k_element_global` matrices.

"""

function define(members, section_properties, material_properties, node_geometry, supports, external_forces)

    truss = truss_object()

    truss.members = members
    truss.section_properties = section_properties
    truss.material_properties = material_properties
    truss.node_geometry = node_geometry
    truss.supports = supports

    num_nodes = size(node_geometry)[1]
    truss.external_forces = zeros(Float64, num_nodes * 2)
    truss.external_forces .= external_forces

    #Define the number of degrees of freedom per node.
    truss.num_dof_per_node = 2

    #Calculate truss member lengths.
    truss.L = calculate_truss_element_lengths(members, node_geometry)

    #Calculate truss member orientations.
    truss.θ = define_truss_element_orientations(members, node_geometry)

    #Define truss properties for stiffness matrix calculations.
    truss.A = assign_member_properties(members, section_properties, 3, 1)
    truss.E = assign_member_properties(members, material_properties, 4, 1)

    #Calculate the local stiffness matrix for each element.
    truss.k_element_local = define_local_truss_element_stiffness_matrix.(truss.A, truss.E, truss.L)

    #Define rotation matrix for each element.
    truss.T = vector_rotation_operator.(truss.θ)

    #Rotate local stiffness matrix into global coordinates.
    truss.k_element_global = calculate_global_element_stiffness_matrix.(truss.k_element_local, truss.T)

    #Assemble the global stiffness matrix.
    truss.K = assemble_global_stiffness_matrix(node_geometry, members, truss.k_element_global, truss.num_dof_per_node)

    #Define free degrees of freedom.
    truss.f = findall(x-> x==0, supports)

    #Define number of free dof.
    num_free_dof = length(truss.f)

    #Define fixed degrees of freedom.
    truss.s = findall(x-> x==1, supports)

    #Define number of fixed dof.
    num_fixed_dof = length(truss.s)

    #Initialize partitioned stiffness matrices.
    truss.Kff = zeros(Float64, (num_free_dof, num_free_dof))
    truss.Ksf = zeros(Float64, (num_fixed_dof, num_free_dof))
    truss.Kfs = zeros(Float64, (num_free_dof, num_fixed_dof))
    truss.Kss = zeros(Float64, (num_fixed_dof, num_fixed_dof))

    #Partition the global stiffness matrix.
    truss.Kff .= truss.K[truss.f, truss.f]
    truss.Ksf .= truss.K[truss.s, truss.f]
    truss.Kfs .= truss.K[truss.f, truss.s]
    truss.Kss .= truss.K[truss.s, truss.s]

    #Partition the external force vectors.
    truss.F = external_forces
    truss.Ff = zeros(Float64, length(truss.f))
    truss.Ff .= truss.F[truss.f]

    return truss

end

"""
    solve(K, Kff, Ksf, F, Ff, s, f, T, k_element_local)

Solve for system displacements, reactions, and internal forces.

`K` is the global stiffness matrix.
'Kff' is the partitioned portion of the global stiffness matrix, free by free.
'Ksf' is the partitioned portion of the global stiffness matrix, support by free.
`F` is an array of the external forces.
`s` is an array of the fixed degrees of freedom.
`f` is an array of the support degrees of freedom.
`T` are all the rotation matrices for each element.
`k_element_local` are all the element local stiffness matrices.

Calculate the global displacement vector `uf` for the free degrees of freedom, provide the full global displacement vector `u`, and calculate all the element internal forces as an array of arrays `f_element`.

"""

function solve(truss)

    #Define the total number of degrees of freedom in the system.
    num_dof = size(truss.K)[1]

    #Initialize the global displacement vector.
    truss.u = zeros(Float64, num_dof)

    #Calculate the global displacements at the free degrees of freedom.
    truss.uf = truss.Kff^-1 * truss.Ff

    #Insert uf into u.
    truss.u[truss.f] .= truss.uf

    #Calculate the global system reactions.
    truss.Fs = truss.Ksf * truss.uf

    #Insert Fs into F.
    truss.F[truss.s] .= truss.Fs

    #Define the number of elements.
    num_elem = length(truss.A)

    #Initialize an array of element force vectors.
    truss.f_element = fill(Float64[], num_elem, 1)

    #Calculate the internal forces in each element.
    for i = 1:num_elem

        #Define node i and j for an element.
        node_i = truss.members[i][1]
        node_j = truss.members[i][2]

        #Calculate the global dof numbers at each element nodes i and j.
        node_i_dof = [1;2] .+ (node_i - 1) * truss.num_dof_per_node
        node_j_dof = [1;2] .+ (node_j - 1) * truss.num_dof_per_node

        #Assemble all the global dof numbers for the element.
        global_dof = [node_i_dof; node_j_dof]

        #Find the element displacements in global coordinates.
        u_element_global = zeros(Float64, length(global_dof))
        u_element_global .= truss.u[global_dof]

        #Transform the element displacements into local coordinates.
        u_element_local = truss.T[i] * u_element_global

        #Calculate the element local forces.
        truss.f_element[i] = truss.k_element_local[i] * u_element_local

    end

    return truss

end

function show(truss_object, scale_x, scale_y)

	δx = truss_object.u[1:2:end]
	δy = truss_object.u[2:2:end]

	num_nodes = floor(Int, length(truss_object.u)/2)

	g = zeros(Int64, (num_nodes, num_nodes))

	start_nodes = [x[1] for x in truss_object.members]
	end_nodes = [x[2] for x in truss_object.members]

	for i = 1:num_nodes

		index = findall(x->x==i, start_nodes)

		g[i, end_nodes[index]] .= 1

	end

	graphplot(g,
          x=truss_object.node_geometry[:,1], y=truss_object.node_geometry[:,2],
          nodeshape=:circle, nodesize=2,
          axis_buffer=0.01,
          curves=false,
          color=:red,
          linewidth=0.5,
	      arrowlengthfrac=0)

	graphplot!(g,
          x=truss_object.node_geometry[:,1] .+ δx .* scale_x, y=truss_object.node_geometry[:,2] .+ δy .* 	scale_y,
          nodeshape=:circle, nodesize=2,
          axis_buffer=0.01,
          curves=false,
          color=:black,
          linewidth=2,
	      arrowlengthfrac=0)

	quiver!(truss_object.node_geometry[:,1],truss_object.node_geometry[:,2],quiver=(truss_object.external_forces[1:2:end]./10, truss_object.external_forces[2:2:end]./10))

	return current()

end

end #module

ext: jl


text:
#!/usr/bin/env julia
################################################################################
#
#  Copyright (c) 2015 Wojciech Migda
#  All rights reserved
#  Distributed under the terms of the Apache 2.0 license
#
################################################################################
#
#  Filename: jl_gen_mfcc_feat.jl
#
#  Decription:
#       Generate MFCC feature vectors
#
#  Authors:
#       Wojciech Migda
#
################################################################################
#
#  History:
#  --------
#  Date         Who  Ticket     Description
#  ----------   ---  ---------  ------------------------------------------------
#  2015-09-16   wm              Initial version
#
################################################################################

addprocs(div(CPU_CORES, 2))

using ArgParse

using HDF5
@everywhere using Worker

function generate(
    MP3_DIR::String,
    MP3_H5_FNAME::String,
    MP3_H5_DSET::String,
    MFCC_H5_FNAME::String,
    MFCC_H5_DSET::String,
    NJOBS::Int,
    NMFCC::Int,
    NFILT::Int,
    ENERGY::Bool
    )

    const TRAIN_MP3 = MP3_DIR

    const TRAIN_X = h5read(MP3_H5_FNAME, MP3_H5_DSET)

    const TRAIN_SUBSETS = int(linspace(0, size(TRAIN_X, 1), (NJOBS + 1)))
    #const TRAIN_SUBSETS = int(linspace(0, 30, 5))
    println(TRAIN_SUBSETS)

    foo = @parallel reduce for i = 1:length(TRAIN_SUBSETS) - 1
        Worker.gen_mfcc(
            MP3_DIR,
            sub(TRAIN_X, TRAIN_SUBSETS[i] + 1:TRAIN_SUBSETS[i + 1]),
            TRAIN_SUBSETS[i] + 1,
            _nfilt=NFILT,
            _numcep=NMFCC,
            _appendEnergy=ENERGY
            )
    end

    gc()

    println(typeof(foo))
    const NROWS = TRAIN_SUBSETS[end]
    println(size(foo[][1][], 1))

    fid = h5open(MFCC_H5_FNAME, "w")

    Xdset = d_create(fid, MFCC_H5_DSET, datatype(Float64), dataspace(size(foo[][1][], 1), NROWS))

    println("Saving Xdset")

    for vt = 1:size(foo, 1)
        pos = foo[vt][2]
        for vx = 1:size(foo[vt][1], 1)
            Xdset[:, pos] = foo[vt][1][vx]
            pos += 1
        end

        gc()
    end

    println("Xdset saved")

    gc()

    close(Xdset)
    close(fid)

    println("Generating done")
end

function parse_commandline(THIS_DIR::ASCIIString)
    s = ArgParseSettings()

    @add_arg_table s begin
        #"--opt1"
        #    help = "an option with an argument"
        "--njobs", "-j"
            help = "number of jobs to schedule for parallel execution"
            arg_type = Int
            default = 1
        "--nmfcc", "-m"
            help = "number of MFCCs"
            arg_type = Int
            default = 13
        "--nfilt", "-f"
            help = "number of MFCC filterbank filters"
            arg_type = Int
            default = 26
        "--energy", "-e"
            help = "replace first MFCCs with spectral energies"
            arg_type = Bool
            default = false
        "--mp3-dir"
            help = "input MP3 files folder"
            arg_type = String
            required = true
        "--h5-mp3-file"
            help = "input H5 store with MP3 filenames"
            arg_type = String
            required = true
        "--h5-mp3-dset"
            help = "input H5 dataset with MP3 filenames"
            arg_type = String
            required = true
        "--h5-mfcc-file"
            help = "output H5 store for MFCC data"
            arg_type = String
            required = true
        "--h5-mfcc-dset"
            help = "output H5 dataset for MFCC data"
            arg_type = String
            required = true

        #"arg1"
        #    help = "a positional argument"
        #    required = true
    end

    return parse_args(s)
end

function main()
    const DATA_DIR = "$(THIS_DIR)/../../data"

    parsed_args = parse_commandline(THIS_DIR)
    println("Parsed args:")
    for (arg,val) in parsed_args
        println("  $arg  =>  $val")
    end

    generate(parsed_args["mp3-dir"], parsed_args["h5-mp3-file"], parsed_args["h5-mp3-dset"], parsed_args["h5-mfcc-file"], parsed_args["h5-mfcc-dset"], parsed_args["njobs"], parsed_args["nmfcc"], parsed_args["nfilt"], parsed_args["energy"])
    #generate("$(DATA_DIR)/test", "testXmp3.h5", "lid/test/X/mp3", "testXmfcc.h5", "lid/test/X/mfcc")

end

const THIS_DIR = dirname(Base.source_path())

if ~isinteractive()
    main()
end


text:
using SimpleSDMLayers
using Test

global anyerrors = false

tests = [
   "construction" => "construction.jl",
   "basics" => "basics.jl",
   "matching" => "matching.jl",
   "overloads" => "overloads.jl",
   "ascii" => "ascii.jl",
   "iteration" => "iteration.jl",
   "rescale" => "rescale.jl",
   "mosaic" => "mosaic.jl",
   "generated" => "generated.jl",
   "import" => "dataread.jl",
   "worldclim" => "worldclim.jl",
   "landcover" => "landcover.jl",
   "chelsa" => "chelsa.jl",
   "coarsen" => "coarsen.jl",
   "plotting" => "plots.jl",
   "GBIF" => "gbif.jl",
   "DataFrames" => "dataframes.jl"
]

for test in tests
   try
      include(test.second)
      println("\033[1m\033[32m✓\033[0m\t$(test.first)")
   catch e
      global anyerrors = true
      println("\033[1m\033[31m×\033[0m\t$(test.first)")
      println("\033[1m\033[38m→\033[0m\ttest/$(test.second)")
      showerror(stdout, e, backtrace())
      println()
      break
   end
end

if anyerrors
   throw("Tests failed")
end


text:
#
#   WORK IN PROGRESS - CURRENTLY NOT INCLUDED
#
#

module SNAP

"""
`SNAPRnlBasis:` define radial basis for hyperspherical harmonics.
```math
R_{nl}(r) = (-i)^l 2^{l+1/2} l! \cdot
      \bigg( \frac{(n+1) (n-l)!}{\pi (n+l+1)!} \bigg)^{1/2} \cdot
      \sin^l\Big( \frac{\pi r}{2 r_0}\Big) \cdot
      C_{n-l}^{l+1}\Big( \cos\Big(\frac{\pi r}{2 r_0}\Big)\Big).
```
"""
struct SNAPRnlBasis{T, VSYM, LSYM, MSYM}
   maxn::Int
   maxl::Int
   spec
   B_pool::VectorPool{Complex{T}}
   dB_pool::VectorPool{T}
end

_varsym(::SNAPRnlBasis{T, VSYM, LSYM, MSYM}) where {T, VSYM, LSYM, MSYM} = VSYM
_nsym(::SNAPRnlBasis{T, VSYM, LSYM, MSYM}) where {T, VSYM, LSYM, MSYM} = NSYM
_n(b, basis::SNAPRnlBasis) = getproperty(b, _nsym(basis))
_rr(X, Rn::SNAPRnlBasis) = getproperty(X, _varsym(Rn))

function evaluate!(Rnl, basis::SNAPRnlBasis, X::AbstractState)
   rr = _rr(X)
   r = norm(rr)

end

end


text:
using Query
using DataFrames
using Test

@testset "Standalone Syntax" begin

@testset "@take operator" begin
    df = DataFrame(a=[1,2,3], b=[3.,2.,1.], c=["a", "b", "c"])

    df2 = df |> @take(2) |> DataFrame

    @test df2 isa DataFrame
    @test size(df2) == (2,3)
    @test df2[:a] == [1,2]
    @test df2[:b] == [3.,2.]
    @test df2[:c] == ["a", "b"]

    df2 = DataFrame(@take(df, 2))

    @test df2 isa DataFrame
    @test size(df2) == (2,3)
    @test df2[:a] == [1,2]
    @test df2[:b] == [3.,2.]
    @test df2[:c] == ["a", "b"]
end

@testset "@drop operator" begin
    df = DataFrame(a=[1,2,3], b=[3.,2.,1.], c=["a", "b", "c"])

    df2 = df |> @drop(1) |> DataFrame

    @test df2 isa DataFrame
    @test size(df2) == (2,3)
    @test df2[:a] == [2,3]
    @test df2[:b] == [2.,1.]
    @test df2[:c] == ["b","c"]

    df2 = DataFrame(@drop(df, 1))

    @test df2 isa DataFrame
    @test size(df2) == (2,3)
    @test df2[:a] == [2,3]
    @test df2[:b] == [2.,1.]
    @test df2[:c] == ["b","c"]
end

end


text:
global dfg,v1,v2,f1

if typeof(dfg) <: CloudGraphsDFG
    @warn "TEST: Nuking all data for user '$(dfg.userId)'!"
    clearUser!!(dfg)
    createDfgSessionIfNotExist(dfg)
end

# Building simple graph...
@testset "Building a simple Graph" begin
    global dfg,v1,v2,f1
    # Use IIF to add the variables and factors
    v1 = addVariable!(dfg, :a, ContinuousScalar, labels = [:POSE])
    v2 = addVariable!(dfg, :b, ContinuousScalar, labels = [:LANDMARK])
    f1 = addFactor!(dfg, [:a; :b], LinearConditional(Normal(50.0,2.0)) )
    # f1 = addFactor!(fg,[:x0], Prior( pd ) )
    # @test_throws Exception addFactor!(dfg, DFGFactor{Int, :Symbol}("f2"), [v1, DFGVariable("Nope")])
end

#test before anything changes
@testset "Producing Dot Files" begin
    global dfg
    todotstr = toDot(dfg)
    #TODO consider using a regex, but for now test both orders
    todota = todotstr == "graph graphname {\n2 [\"label\"=\"a\",\"shape\"=\"ellipse\",\"fillcolor\"=\"red\",\"color\"=\"red\"]\n2 -- 3\n3 [\"label\"=\"abf1\",\"shape\"=\"box\",\"fillcolor\"=\"blue\",\"color\"=\"blue\"]\n1 [\"label\"=\"b\",\"shape\"=\"ellipse\",\"fillcolor\"=\"red\",\"color\"=\"red\"]\n1 -- 3\n}\n"
    todotb = todotstr == "graph graphname {\n2 [\"label\"=\"b\",\"shape\"=\"ellipse\",\"fillcolor\"=\"red\",\"color\"=\"red\"]\n2 -- 3\n3 [\"label\"=\"abf1\",\"shape\"=\"box\",\"fillcolor\"=\"blue\",\"color\"=\"blue\"]\n1 [\"label\"=\"a\",\"shape\"=\"ellipse\",\"fillcolor\"=\"red\",\"color\"=\"red\"]\n1 -- 3\n}\n"
    @test (todota || todotb)
    @test toDotFile(dfg, "something.dot") == nothing
    Base.rm("something.dot")
end

@testset "Testing CRUD, return and Failures from a GraphsDFG" begin
    global dfg
    # fg to copy to
    # creating a whole new graph with the same labels
    T = typeof(dfg)
    if T <: CloudGraphsDFG
        dfg2 = CloudGraphsDFG{SolverParams}("localhost", 7474, "neo4j", "test",
                                            "testUser", "testRobot", "testSession2",
                                            nothing,
                                            nothing,
                                            IncrementalInference.decodePackedType,
                                            IncrementalInference.rebuildFactorMetadata!,
                                            solverParams=SolverParams())
    else
        dfg2 = T()
    end

    # Build a new in-memory IIF graph to transfer into the new graph.
    iiffg = initfg()
    v1 = deepcopy(addVariable!(iiffg, :a, ContinuousScalar))
    v2 = deepcopy(addVariable!(iiffg, :b, ContinuousScalar))
    v3 = deepcopy(addVariable!(iiffg, :c, ContinuousScalar))
    f1 = deepcopy(addFactor!(iiffg, [:a; :b], LinearConditional(Normal(50.0,2.0)) ))
    f2 = deepcopy(addFactor!(iiffg, [:b; :c], LinearConditional(Normal(10.0,1.0)) ))

    # Add it to the new graph.
    @test addVariable!(dfg2, v1)
    @test addVariable!(dfg2, v2)
    @test_throws ErrorException updateVariable!(dfg2, v3)
    @test addVariable!(dfg2, v3)
    @test_throws ErrorException addVariable!(dfg2, v3)
    @test addFactor!(dfg2, [v1, v2], f1)
    @test_throws ErrorException addFactor!(dfg2, [v1, v2], f1)
    @test_throws ErrorException updateFactor!(dfg2, f2)
    @test addFactor!(dfg2, [:b, :c], f2)

    dv3 = deleteVariable!(dfg2, v3)
    #TODO write compare if we want to compare complete one, for now just label
    # @test dv3 == v3
    @test dv3.label == v3.label
    @test_throws ErrorException deleteVariable!(dfg2, v3)

    @test symdiff(ls(dfg2),[:a,:b]) == []
    df2 = deleteFactor!(dfg2, f2)
    #TODO write compare if we want to compare complete one, for now just label
    # @test df2 == f2
    @test df2.label == f2.label
    @test_throws ErrorException deleteFactor!(dfg2, f2)

    @test lsf(dfg2) == [:abf1]

end

@testset "Listing Nodes" begin
    global dfg,v1,v2,f1
    @test length(ls(dfg)) == 2
    @test length(lsf(dfg)) == 1 # Unless we add the prior!
    @test symdiff([:a, :b], getVariableIds(dfg)) == []
    @test getFactorIds(dfg) == [:abf1] # Unless we add the prior!
    #
    @test lsf(dfg, :a) == [f1.label]
    # Tags
    @test ls(dfg, tags=[:POSE]) == [:a]
    @test symdiff(ls(dfg, tags=[:POSE, :LANDMARK]), ls(dfg, tags=[:VARIABLE])) == []
    # Regexes
    @test ls(dfg, r"a") == [v1.label]
    # TODO: Check that this regular expression works on everything else!
    # it works with the .
    # REF: https://stackoverflow.com/questions/23834692/using-regular-expression-in-neo4j
    @test lsf(dfg, r"abf.*") == [f1.label]

    # Accessors
    @test getAddHistory(dfg) == [:a, :b] #, :abf1
    @test getDescription(dfg) != nothing
    @test getLabelDict(dfg) != nothing
    # Existence
    @test exists(dfg, :a) == true
    @test exists(dfg, v1) == true
    @test exists(dfg, :nope) == false
    # Sorting of results
    # TODO - this function needs to be cleaned up
    unsorted = [:x1_3;:x1_6;:l1;:april1] #this will not work for :x1x2f1
    @test sortDFG(unsorted) == sortVarNested(unsorted)
    @test_skip sortDFG([:x1x2f1, :x1l1f1]) == [:x1l1f1, :x1x2f1]
end

# Gets
@testset "Gets, Sets, and Accessors" begin
    global dfg,v1,v2,f1
    @test getVariable(dfg, v1.label) == v1
    @test getFactor(dfg, f1.label) == f1
    @test_throws Exception getVariable(dfg, :nope)
    @test_throws Exception getVariable(dfg, "nope")
    @test_throws Exception getFactor(dfg, :nope)
    @test_throws Exception getFactor(dfg, "nope")

    # Sets
    v1Prime = deepcopy(v1)
    @test updateVariable!(dfg, v1Prime) == v1 #Maybe move to crud
    @test updateVariable!(dfg, v1Prime) == getVariable(dfg, v1.label)
    f1Prime = deepcopy(f1)
    @test updateFactor!(dfg, f1Prime) == f1 #Maybe move to crud
    @test updateFactor!(dfg, f1Prime) == getFactor(dfg, f1.label)

    # Accessors
    @test label(v1) == v1.label
    @test tags(v1) == v1.tags
    @test timestamp(v1) == v1.timestamp
    @test estimates(v1) == v1.estimateDict
    @test DistributedFactorGraphs.estimate(v1, :notfound) == nothing
    @test solverData(v1) === v1.solverDataDict[:default]
    @test getData(v1) === v1.solverDataDict[:default]
    @test solverData(v1, :default) === v1.solverDataDict[:default]
    @test solverDataDict(v1) == v1.solverDataDict
    @test internalId(v1) == v1._internalId

    @test softtype(v1) == :ContinuousScalar#Symbol(typeof(st1))
    @test softtype(v2) == :ContinuousScalar#Symbol(typeof(st2))
    @test typeof(getSofttype(v1)) == typeof(ContinuousScalar())

    @test label(f1) == f1.label
    @test tags(f1) == f1.tags
    @test solverData(f1) == f1.data
    # Deprecated functions
    @test data(f1) == f1.data
    @test getData(f1) == f1.data
    # Internal function
    @test internalId(f1) == f1._internalId

    @test getSolverParams(dfg) != nothing
    @test setSolverParams(dfg, getSolverParams(dfg)) == getSolverParams(dfg)

    #solver data is initialized
    @test !isInitialized(dfg, :a)
    @test !isInitialized(v2)

    @test !isInitialized(v2, key=:second)

    # Session, robot, and user small data tests
    # NOTE: CloudGraphDFG isnt supporting this yet.
    if !(typeof(dfg) <: CloudGraphsDFG)
        smallUserData = Dict{Symbol, String}(:a => "42", :b => "Hello")
        smallRobotData = Dict{Symbol, String}(:a => "43", :b => "Hello")
        smallSessionData = Dict{Symbol, String}(:a => "44", :b => "Hello")
        setUserData(dfg, deepcopy(smallUserData))
        setRobotData(dfg, deepcopy(smallRobotData))
        setSessionData(dfg, deepcopy(smallSessionData))
        @test getUserData(dfg) == smallUserData
        @test getRobotData(dfg) == smallRobotData
        @test getSessionData(dfg) == smallSessionData
    end
end

@testset "BigData" begin
    # NOTE: CloudGraphDFG isnt supporting this yet.
    if !(typeof(dfg) <: CloudGraphsDFG)
        oid = zeros(UInt8,12); oid[12] = 0x01
        de1 = MongodbBigDataEntry(:key1, NTuple{12,UInt8}(oid))

        oid = zeros(UInt8,12); oid[12] = 0x02
        de2 = MongodbBigDataEntry(:key2, NTuple{12,UInt8}(oid))

        oid = zeros(UInt8,12); oid[12] = 0x03
        de2_update = MongodbBigDataEntry(:key2, NTuple{12,UInt8}(oid))

        #add
        v1 = getVariable(dfg, :a)
        @test addBigDataEntry!(v1, de1)
        @test addBigDataEntry!(dfg, :a, de2)
        @test addBigDataEntry!(v1, de1)

        #get
        @test deepcopy(de1) == getBigDataEntry(v1, :key1)
        @test deepcopy(de2) == getBigDataEntry(dfg, :a, :key2)
        @test_throws Any getBigDataEntry(v2, :key1)
        @test_throws Any getBigDataEntry(dfg, :b, :key1)

        #update
        @test updateBigDataEntry!(dfg, :a, de2_update)
        @test deepcopy(de2_update) == getBigDataEntry(dfg, :a, :key2)
        @test !updateBigDataEntry!(dfg, :b, de2_update)

        #list
        entries = getBigDataEntries(dfg, :a)
        @test length(entries) == 2
        @test symdiff(map(e->e.key, entries), [:key1, :key2]) == Symbol[]
        @test length(getBigDataEntries(dfg, :b)) == 0

        @test symdiff(getBigDataKeys(dfg, :a), [:key1, :key2]) == Symbol[]
        @test getBigDataKeys(dfg, :b) == Symbol[]

        #delete
        @test deepcopy(de1) == deleteBigDataEntry!(v1, :key1)
        @test getBigDataKeys(v1) == Symbol[:key2]
        #delete from dfg
        @test deepcopy(de2_update) == deleteBigDataEntry!(dfg, :a, :key2)
        @test getBigDataKeys(v1) == Symbol[]
    end
end

@testset "Updating Nodes and Estimates" begin
    global dfg
    #get the variable
    var = getVariable(dfg, :a)
    #make a copy and simulate external changes
    newvar = deepcopy(var)
    estimates(newvar)[:default] = MeanMaxPPE(:default, [150.0], [100.0], [50.0])
    #update
    mergeUpdateVariableSolverData!(dfg, newvar)

    #Check if variable is updated
    var = getVariable(dfg, :a)
    @test estimates(newvar) == estimates(var)

    # Add a new estimate.
    estimates(newvar)[:second] = MeanMaxPPE(:second, [15.0], [10.0], [5.0])

    # Confirm they're different
    @test estimates(newvar) != estimates(var)
    # Persist it.
    mergeUpdateVariableSolverData!(dfg, newvar)
    # Get the latest
    var = getVariable(dfg, :a)
    @test symdiff(collect(keys(estimates(var))), [:default, :second]) == Symbol[]

    #Check if variable is updated
    @test estimates(newvar) == estimates(var)

    # Delete :default and replace to see if new ones can be added
    delete!(estimates(newvar), :default)
    #confirm delete
    @test symdiff(collect(keys(estimates(newvar))), [:second]) == Symbol[]
    # Persist it.
    mergeUpdateVariableSolverData!(dfg, newvar)

    # Get the latest and confirm they're the same, :second
    var = getVariable(dfg, :a)

    # TODO issue #166
    @test estimates(newvar) != estimates(var)
    @test collect(keys(estimates(var))) ==  [:default, :second]

    # @test symdiff(collect(keys(estimates(getVariable(dfg, :a)))), [:default, :second]) == Symbol[]
end

# Connectivity test
@testset "Connectivity Test" begin
    global dfg,v1,v2,f1
    @test isFullyConnected(dfg) == true
    @test hasOrphans(dfg) == false
    addVariable!(dfg, :orphan, ContinuousScalar, labels = [:POSE])
    @test isFullyConnected(dfg) == false
    @test hasOrphans(dfg) == true
end

# Adjacency matrices
@testset "Adjacency Matrices" begin
    global dfg,v1,v2,f1
    # Normal
    adjMat = getAdjacencyMatrix(dfg)
    @test size(adjMat) == (2,4)
    @test symdiff(adjMat[1, :], [nothing, :a, :b, :orphan]) == Symbol[]
    @test symdiff(adjMat[2, :], [:abf1, :abf1, :abf1, nothing]) == Symbol[]
    #sparse
    adjMat, v_ll, f_ll = getAdjacencyMatrixSparse(dfg)
    @test size(adjMat) == (1,3)

    # Checking the elements of adjacency, its not sorted so need indexing function
    indexOf = (arr, el1) -> findfirst(el2->el2==el1, arr)
    @test adjMat[1, indexOf(v_ll, :orphan)] == 0
    @test adjMat[1, indexOf(v_ll, :a)] == 1
    @test adjMat[1, indexOf(v_ll, :b)] == 1
    @test symdiff(v_ll, [:a, :b, :orphan]) == Symbol[]
    @test symdiff(f_ll, [:abf1, :abf1, :abf1]) == Symbol[]
end

# Deletions
@testset "Deletions" begin
    deleteFactor!(dfg, :abf1)
    @test getFactorIds(dfg) == []
    deleteVariable!(dfg, :b)
    @test symdiff([:a, :orphan], getVariableIds(dfg)) == []
    #delete last also for the LightGraphs implementation coverage
    deleteVariable!(dfg, :orphan)
    @test symdiff([:a], getVariableIds(dfg)) == []
    deleteVariable!(dfg, :a)
    @test getVariableIds(dfg) == []
end

# Now make a complex graph for connectivity tests
numNodes = 10
#the deletions in last test should have cleared out the fg
# dfg = DistributedFactorGraphs._getDuplicatedEmptyDFG(dfg)
# if typeof(dfg) <: CloudGraphsDFG
#     clearSession!!(dfg)
# end

#change ready and backendset for x7,x8 for improved tests on x7x8f1
verts = map(n -> addVariable!(dfg, Symbol("x$n"), ContinuousScalar, labels = [:POSE]), 1:numNodes)
#TODO fix this to use accessors
verts[7].ready = 1
# verts[7].backendset = 0
verts[8].ready = 0
verts[8].backendset = 1
#call update to set it on cloud
updateVariable!(dfg, verts[7])
updateVariable!(dfg, verts[8])

facts = map(n -> addFactor!(dfg, [verts[n], verts[n+1]], LinearConditional(Normal(50.0,2.0))), 1:(numNodes-1))

@testset "Getting Neighbors" begin
    global dfg,verts
    # Trivial test to validate that intersect([], []) returns order of first parameter
    @test intersect([:x3, :x2, :x1], [:x1, :x2]) == [:x2, :x1]
    # Get neighbors tests
    @test getNeighbors(dfg, verts[1]) == [:x1x2f1]
    neighbors = getNeighbors(dfg, getFactor(dfg, :x1x2f1))
    @test neighbors == [:x1, :x2]
    # Testing aliases
    @test getNeighbors(dfg, getFactor(dfg, :x1x2f1)) == ls(dfg, getFactor(dfg, :x1x2f1))
    @test getNeighbors(dfg, :x1x2f1) == ls(dfg, :x1x2f1)

    # ready and backendset
    @test getNeighbors(dfg, :x5, ready=1) == Symbol[]
    #TODO Confirm: test failed on GraphsDFG, don't know if the order is important for isa variable.
    @test symdiff(getNeighbors(dfg, :x5, ready=0), [:x4x5f1,:x5x6f1]) == []
    @test getNeighbors(dfg, :x5, backendset=1) == Symbol[]
    @test symdiff(getNeighbors(dfg, :x5, backendset=0),[:x4x5f1,:x5x6f1]) == []
    @test getNeighbors(dfg, :x7x8f1, ready=0) == [:x8]
    @test getNeighbors(dfg, :x7x8f1, backendset=0) == [:x7]
    @test getNeighbors(dfg, :x7x8f1, ready=1) == [:x7]
    @test getNeighbors(dfg, :x7x8f1, backendset=1) == [:x8]
    @test getNeighbors(dfg, verts[1], ready=0) == [:x1x2f1]
    @test getNeighbors(dfg, verts[1], ready=1) == Symbol[]
    @test getNeighbors(dfg, verts[1], backendset=0) == [:x1x2f1]
    @test getNeighbors(dfg, verts[1], backendset=1) == Symbol[]

end

@testset "Getting Subgraphs" begin
    # Subgraphs
    dfgSubgraph = getSubgraphAroundNode(dfg, verts[1], 2)
    # Only returns x1 and x2
    @test symdiff([:x1, :x1x2f1, :x2], [ls(dfgSubgraph)..., lsf(dfgSubgraph)...]) == []
    # Test include orphan factorsVoid
    dfgSubgraph = getSubgraphAroundNode(dfg, verts[1], 1, true)
    @test symdiff([:x1, :x1x2f1], [ls(dfgSubgraph)..., lsf(dfgSubgraph)...]) == []
    # Test adding to the dfg
    dfgSubgraph = getSubgraphAroundNode(dfg, verts[1], 2, true, dfgSubgraph)
    @test symdiff([:x1, :x1x2f1, :x2], [ls(dfgSubgraph)..., lsf(dfgSubgraph)...]) == []
    #
    dfgSubgraph = getSubgraph(dfg,[:x1, :x2, :x1x2f1])
    # Only returns x1 and x2
    @test symdiff([:x1, :x1x2f1, :x2], [ls(dfgSubgraph)..., lsf(dfgSubgraph)...]) == []

    # DFG issue #95 - confirming that getSubgraphAroundNode retains order
    # REF: https://github.com/JuliaRobotics/DistributedFactorGraphs.jl/issues/95
    for fId in getVariableIds(dfg)
        # Get a subgraph of this and it's related factors+variables
        dfgSubgraph = getSubgraphAroundNode(dfg, verts[1], 2)
        # For each factor check that the order the copied graph == original
        for fact in getFactors(dfgSubgraph)
            @test fact._variableOrderSymbols == getFactor(dfg, fact.label)._variableOrderSymbols
        end
    end
end

@testset "Summaries and Summary Graphs" begin
    factorFields = fieldnames(DFGFactorSummary)
    variableFields = fieldnames(DFGVariableSummary)

    summary = getSummary(dfg)
    @test symdiff(collect(keys(summary.variables)), ls(dfg)) == Symbol[]
    @test symdiff(collect(keys(summary.factors)), lsf(dfg)) == Symbol[]

    summaryGraph = getSummaryGraph(dfg)
    @test symdiff(ls(summaryGraph), ls(dfg)) == Symbol[]
    @test symdiff(lsf(summaryGraph), lsf(dfg)) == Symbol[]
    # Check all fields are equal for all variables
    for v in ls(summaryGraph)
        for field in variableFields
            if field != :softtypename
            @test getfield(getVariable(dfg, v), field) == getfield(getVariable(summaryGraph, v), field)
            else
                @test softtype(getVariable(dfg, v)) == softtype(getVariable(summaryGraph, v))
        end
    end
    end
    for f in lsf(summaryGraph)
        for field in factorFields
            @test getfield(getFactor(dfg, f), field) == getfield(getFactor(summaryGraph, f), field)
        end
    end
end


text:
module BatchQueryStrategies

using Memento
using Random

import SVDD
import MLLabelUtils
import MLKernels
import JuMP
import Distances
import LinearAlgebra

import ..QueryStrategies:
    MissingLabelTypeException,
    PoolQs,
    get_query_objects,
    multi_kde
import ..SequentialQueryStrategies:
    SequentialPQs,
    qs_score
import IterTools: subsets
import Clustering: kmedoids
import StatsBase: sample

include("batch_qs_base.jl")

include("TopKBatchQs.jl")
include("GappedTopKBatchQs.jl")
include("RandomBatchQs.jl")
include("IterativeBatchQs.jl")
include("EnumFilterHierarchicalBatchQs.jl")
include("FilterHierarchicalBatchQs.jl")
include("ClusterBatchQs.jl")
include("EnumerativeBatchQs.jl")
include("EnsembleBatchQs.jl")
include("ClusterTopKBatchQs.jl")
include("FilterSimilarBatchQs.jl")

export
    BatchPQs,
    ExtendingBatchQs,
    TopKBatchQs, GappedTopKBatchQs,
    RandomBatchQs,
    ClusterBatchQs,

    MultiObjectiveBatchQs,
    IterativeBatchQs,
    FilterHierarchicalBatchQs,
    EnumFilterHierarchicalBatchQs,
    EnumerativeBatchQs,

    EnsembleBatchQs,
    ClusterTopKBatchQs,

    FilterSimilarBatchQs,

    select_batch,
    angle_batch_diversity,
    euclidean_batch_diversity
end


text:
using Printf, Test

include("../../hybmc/termstructures/YieldCurve.jl")
include("../../hybmc/models/DeterministicModel.jl")
include("../../hybmc/simulations/McSimulation.jl")
include("../../hybmc/simulations/Payoffs.jl")

function deterministicPath()
    yc = YieldCurve(0.0)
    model = DeterministicModel("EUR",yc,["USD"],[1.0],[yc])
    sim = McSimulation(model,[0.0],1)
    return Path(sim,1)
end

@testset "test_payoffSetup" begin
    p = deterministicPath()
    a  = 0.5
    px = 1.0
    py = 2.0
    x = Fixed(px)
    y = Fixed(py)
    #
    @test isequal(at(x+y,p), px+py)
    @test isequal(at(x-y,p), px-py)
    @test isequal(at(x*y,p), px*py)
    @test isequal(at(x/y,p), px/py)
    #
    @test isequal(at(x+a,p), px+a)
    @test isequal(at(x-a,p), px-a)
    @test isequal(at(x*a,p), px*a)
    @test isequal(at(x/a,p), px/a)
    #
    @test isequal(at(a+y,p), a+py)
    @test isequal(at(a-y,p), a-py)
    @test isequal(at(a*y,p), a*py)
    @test isequal(at(a/y,p), a/py)
    #
    @test isequal(at(x<y,p),float(px<py))
    @test isequal(at(x<=y,p),float(px<=py))
    @test isequal(at(x==y,p),float(px==py))
    @test isequal(at(x>=y,p),float(px>=py))
    @test isequal(at(x>y,p),float(px>py))
    #
    @test isequal(at(x<a,p),float(px<a))
    @test isequal(at(x<=a,p),float(px<=a))
    @test isequal(at(x==a,p),float(px==a))
    @test isequal(at(x>=a,p),float(px>=a))
    @test isequal(at(x>a,p),float(px>a))
    #
    @test isequal(at(a<y,p),float(a<py))
    @test isequal(at(a<=y,p),float(a<=py))
    @test isequal(at(a==y,p),float(a==py))
    @test isequal(at(a>=y,p),float(a>=py))
    @test isequal(at(a>y,p),float(a>py))
    #
    @test isequal(obsTime(x + y), 0.0)
    @test isequal(obsTime((x+y)%1.0), 1.0)
end

max_stars_repo_name: sschlenkrich/HybridMonteCarlo

max_issues_repo_name: sschlenkrich/HybridMonteCarlo

max_forks_repo_name: sschlenkrich/HybridMonteCarlo


text:
using MOSLab
using Plots
using LaTeXStrings
NpolyDoping = 5.5e17 # Gate Npoly doping concentration in cm^-3
N_a = 5.5e17 # Bulk doping concentration in cm^-3
E_a = 0.0# The ground state energy of the acceptor in eV
t_ox = 3.8e-7 # oxide thickness in cm
t_si = 1e-6
CList2 = reshape( range(colorant"blue", stop=colorant"red",length=7), 1, 7 );

MOSf(T) = MOSStructure(NPoly(NpolyDoping),SiO2(),SemiconductorData(T,BoltzmanDist(),PSilicon(N_a,E_a)),t_ox,t_si) ## Calculate Parameters of a MOS Structure the given parameters using Boltzman Distribution at temperature T
transistor_model(T) = BSIM6Model(MOSf(T),1.0e-4,1.0e-4)
Idf(Vg,Vd,T) = Id(Vg,Vd,0.0,transistor_model(T)) ## Calculate the Current using the ACM Model of a transistor having the parameters from MOSf(T) and W =1.0 um, L = 1.0 um, similar contructors are available for the other models
plot()
for Vgv in 0.2:0.2:1.8
    plot!(0:0.05:1.8,x->1e6*Idf(Vgv,x,300.0),label="\$V_{GB} = $(Vgv) V\$") # plot Id, Vd characteristics for different VGB
end
PSPG = plot!(xlabel=L"V_{DS} \ [V]",ylabel=L"I_{DS} \ [\mu A]",legend=:topleft)

plot()
for (i,Tv) in enumerate(LinRange(200,500,5))
    plot!(0:0.05:1.0,x->1e6*Idf(x,0.1,Tv),label="\$T = $(Tv) K\$",linecolor=CList2[i]) # plot Id, Vg characteristics for different Temperatures
end
PSPD = plot!(xlabel=L"V_{GS} \ [V]",ylabel=L"I_{DS} \ [\mu A]",legend=:topleft,yaxis=:log10)
plot(PSPG,PSPD) # plot both graphs side by side
savefig("TransistorCurves.png")

gmIdf(Vg,Vd,T) = gm(Vg,Vd,0.0,transistor_model(T))/Idf(Vg,Vd,T)
plot()
for (i,Tv) in enumerate(LinRange(200,500,5))
    plot!(0:0.01:1.0,x->gmIdf(x,0.1,Tv),label="\$T = $(Tv) K\$",linecolor=CList2[i]) # plot Id, Vg characteristics for different Temperatures
end
PSPD = plot!(xlabel=L"V_{GS} \ [V]",ylabel=L"I_{DS} \ [\mu A]",legend=:topright)
savefig("gmId.png")


text:
function _circular_buffer(A, len=length(A), T=eltype(A))
    buff = CircularBuffer{T}(len)
    for val in A
        push!(buff, val)
    end
    return buff
end

"""
    DiscreteBufferVariable
"""
struct DiscreteBufferVariable{T} <: AbstractVector{T}
    buffer::OffsetArray{T,1,CircularBuffer{T}}
end
function DiscreteBufferVariable(x::AbstractVector{<:T}) where T
    len = length(x)
    buff = _circular_buffer(x, len)
    oa = OffsetArray(buff, -length(x)-1)
    return DiscreteBufferVariable{T}(oa)
end
DiscreteBufferVariable(x, n) = DiscreteBufferVariable(fill(x, n))

@inline Base.getindex(dbv::DiscreteBufferVariable) = getindex(dbv.buffer, -1)

@inline Base.setindex!(dbv::DiscreteBufferVariable, val) = push!(dbv, val)

Base.copy(dbv::DiscreteBufferVariable) = DiscreteBufferVariable(copy(dbv.buffer))

function Base.copyto!(dest::DiscreteBufferVariable, src::AbstractVector{<:T}) where T
    for (i_dest, src_elem) = zip(eachindex(dest), src)
        dest[i_dest] = src_elem
    end
    return dest
end

Base.similar(dbv::DiscreteBufferVariable) = DiscreteBufferVariable(similar(dbv.buffer))
Base.similar(dbv::DiscreteBufferVariable, ::Type{T}) where {T} = DiscreteBufferVariable(similar(dbv.buffer, T))

# Pass through methods to the inner CircularBuffer
for fun in [empty!, pop!, push!, popfirst!, pushfirst!, append!, fill!, size, length, getindex, setindex!, axes]
    m = Base.parentmodule(fun)
    f = nameof(fun)
    @eval $m.$f(dbv::DiscreteBufferVariable, args...; kwargs...) = $m.$f(dbv.buffer, args...; kwargs...)
end


text:
println("===== Universal Driver Initialization BEGIN =====")

if isdir(configs[:caserun])
    cd(configs[:caserun])
else
    throw(ErrorException("Caserun directory [ " * configs[:caserun] * " ] does not exist."))
end

function parseCESMTIME(t_str::AbstractString)
    t = zeros(Integer, 4)
    t[1] = parse(Int, t_str[1:4])
    t[2] = parse(Int, t_str[5:6])
    t[3] = parse(Int, t_str[7:8])
    t[4] = parse(Int, t_str[10:17])
    return t
end

function copy_list_to!(
    a,
    b,
)
    for i = 1:length(a)
        b[i][:] = a[i]
    end
end

output_vars = Dict()

stage = :INIT

PTI = ProgramTunnelInfo(
    reverse_role  = true,
    recv_channels = 2,
)

map = NetCDFIO.MapInfo{Float64}(configs[:domain_file])

t_cnt = 1
output_filename = ""
nullbin  = [zeros(Float64, 1)]

ocn_run_time = 0.0
ocn_run_N    = 0

overall_run_time = 0.0
overall_run_N    = 0

println("===== ", OMMODULE.name, " IS READY =====")

beg_time = Base.time()
end_time = Base.time()
while true

    global OMDATA, stage, t_cnt, output_filename, end_time, timeinfo

    new_end_time = Base.time()
    exe_time = new_end_time - end_time
    end_time = new_end_time

    println(format("# Current real time:      : {:s}", Dates.format(now(), "yyyy/mm/dd HH:MM:SS")))
    println(format("# Time counter for RUN    : {:d}", t_cnt))
    println(format("# Total execution time    : {:.01f} m", (end_time - beg_time)/60.0))
    println(format("# Average execution time  : {:.2f} s", floor(end_time - beg_time) / t_cnt))
    println(format("# Last run execution time : {:.2f} s", exe_time))
    println(format("# Stage                   : {}", String(stage)))

    msg = parseMsg( recvData!(PTI, nullbin, which=1) )

    println("==== MESSAGE RECEIVED ====")
    print(json(msg, 4))
    println("==========================")

    # Update time information. It is better do it separately from
    # other action below to make it clear.
    if msg["MSG"] == "INIT"
        timeinfo = ModelTime.TimeInfo(
            parseCESMTIME(msg["CESMTIME"])
        )
    elseif msg["MSG"] == "RUN"
        ModelTime.update!(
            timeinfo,
            parseCESMTIME(msg["CESMTIME"])
        )
    end

    if stage == :INIT && msg["MSG"] == "INIT"

        println(format("Remove archive list file if it exists: {:s}", configs[:archive_list]))
        rm(configs[:archive_list], force=true)

        println("===== INITIALIZING MODEL: ", OMMODULE.name , " =====")

        OMDATA = OMMODULE.init(
            casename     = configs[:casename],
            map          = map,
            timeinfo     = timeinfo,
            configs      = configs,
            read_restart = (msg["READ_RESTART"] == "TRUE") ? true : false,
        )

        # Must declare SharedArray to avoid implicit copying!!!!
        global send_data_list_shared = SharedArray{Float64}[OMDATA.o2x["SST"], OMDATA.o2x["QFLX2ATM"]]
        global recv_data_list_shared = SharedArray{Float64}[]

        global send_data_list = Array{Float64}[OMDATA.o2x["SST"], OMDATA.o2x["QFLX2ATM"]]
        global recv_data_list = Array{Float64}[]

        global x2o_available_varnames  = split(msg["VAR2D"], ",")
        global x2o_wanted_varnames = keys(OMDATA.x2o)
        global x2o_wanted_flag     = [(x2o_available_varnames[i] in x2o_wanted_varnames) for i = 1:length(x2o_available_varnames)]

        println("List of available x2o variables:")
        for (i, varname) in enumerate(x2o_available_varnames)
            push!(recv_data_list_shared, ( x2o_wanted_flag[i] ) ? OMDATA.x2o[varname] :  SharedArray{Float64}(map.lsize))
            push!(recv_data_list       , ( x2o_wanted_flag[i] ) ? OMDATA.x2o[varname] :  zeros(Float64, map.lsize))
            println(format(" ({:d}) {:s} => {:s}", i, varname, ( x2o_wanted_flag[i] ) ? "Wanted" : "Abandoned" ))
        end

        copy_list_to!(send_data_list_shared, send_data_list)

        sendData(PTI, "OK", send_data_list)
        sendData(PTI, "MASKDATA", [map.mask])

        stage = :RUN

    elseif stage == :RUN && msg["MSG"] == "RUN"

        recvData!(
            PTI,
            recv_data_list,
            which=2
        )

        copy_list_to!(recv_data_list, recv_data_list_shared)

        println("Calling ", OMMODULE.name, " to do MAGICAL calculations")

        Δt = parse(Float64, msg["DT"])

        cost = @elapsed let

            OMMODULE.run!(
                OMDATA;
                Δt            = Δt,
                write_restart = msg["WRITE_RESTART"] == "TRUE",
            )

        end

        global ocn_run_time += cost
        global ocn_run_N += 1

        println(format("*** Ocean takes {:.2f} secs. (Avg: {:.2f} secs) ***", cost, ocn_run_time / ocn_run_N))

        copy_list_to!(send_data_list_shared, send_data_list)
        sendData(PTI, "OK", send_data_list)

        t_cnt += 1

    elseif stage == :RUN && msg["MSG"] == "END"

        # move short_term_archive_files to long term archive directory
        println("Waiting 30 seconds for CESM to create archive directory completely.")
        println("Hacking into CESM code will not be sustainable.")
        for t=1:30
            print(format("\r... {:d} second(s) ...", t))
            sleep(1)
        end
        println()
        println("Waiting complete.")
        println("===== Archiving files BEGIN =====")
        sdir = configs[:caseroot]
        ldir = configs[:archive_root]
        mkpath(ldir)
        for line in eachline(configs[:archive_list])

            args = split(line, ",")

            if length(args) == 0
                continue
            end

            action = args[1]
            args = args[2:end]

            if action in ["mv", "cp"]

                fname, src_dir, dst_dir = args

                if ! isdir(dst_dir)
                    mkpath(dst_dir)
                end

                src_file = joinpath(src_dir, fname)
                dst_file = joinpath(dst_dir, fname)

                if isfile(src_file)

                    if action == "mv"
                        mv(src_file, dst_file, force=true)
                        println(format("Moving file: {:s} ( {:s} => {:s} )", fname, src_dir, dst_dir))
                    elseif action == "cp"
                        cp(src_file, dst_file, force=true)
                        println(format("Copying file: {:s} ( {:s} => {:s} )", fname, src_dir, dst_dir))
                    end

                else
                    println("File does not exist: ", src_file)
                end

            elseif action == "rm"
                fname, fdir = args
                rm(joinpath(fdir, fname), force=true)
                println(format("Removing file: {:s} in {:s}", fname, fdir))
            else
                throw(ErrorException(format("Unknown action in archive list: {:s}", action)))
            end

        end

        println("===== Archiving files END =====")

        OMMODULE.final(OMDATA)

        # Print some report... ?

        println("Simulation ends peacefully.")
        break
    else
        #OMMODULE.crash(OMDATA)
        sendData(PTI, "CRASH", send_data_list)
        throw(ErrorException("Unknown status: stage " * string(stage) * ", MSG: " * string(msg["MSG"])))
    end

    flush(stdout)
end

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
using Dates
using Parameters
using DataStructures

##########################################
## TSOActions
##########################################

@with_kw struct TSOActions
    #gen_id, ts => Float (values are firm)
    # Max allowed P of a limitable for a given ts
    limitations::SortedDict{Tuple{String, Dates.DateTime}, Float64} =
        SortedDict{Tuple{String, Dates.DateTime}, Float64}()
    # Imposed P bounds of an pilotable for a given ts
    impositions::SortedDict{Tuple{String, Dates.DateTime}, UncertainValue{Tuple{Float64,Float64}}} =
        SortedDict{Tuple{String, Dates.DateTime}, UncertainValue{Tuple{Float64,Float64}}}()
end

function reset_tso_actions!(tso_actions::TSOActions)
    empty!(tso_actions.limitations)
    empty!(tso_actions.impositions)

    return tso_actions
end

## Limitations
#--------------

function get_limitations(tso_actions::TSOActions)
    return tso_actions.limitations
end
function get_limitations(limitations_dict::SortedDict{Tuple{String, Dates.DateTime}, Float64})
    return limitations_dict
end

function get_limitation(tso_actions, gen_id::String, ts::Dates.DateTime)::Union{Float64, Missing}
    limitations = get_limitations(tso_actions)
    if !haskey(limitations, (gen_id, ts))
        return missing
    else
        return limitations[gen_id, ts]
    end
end
function safeget_limitation(tso_actions, gen_id::String, ts::Dates.DateTime)::Union{Float64, Missing}
    limitation = get_limitation(tso_actions, gen_id, ts)
    if ismissing(limitation)
        msg = @sprintf("Missing limitation value for (gen_id=%s,ts=%s).", gen_id, ts)
        throw( error(msg) )
    else
        return limitations[gen_id, ts]
    end
end

function set_limitation_value!(tso_actions, gen_id::String, ts::Dates.DateTime, value::Float64)
    limitations = get_limitations(tso_actions)
    limitations[gen_id, ts] = value
end

## Impositions
#--------------

function get_impositions(tso_actions::TSOActions)
    return tso_actions.impositions
end

function get_imposition_uncertain_value(tso_actions::TSOActions, gen_id::String, ts::Dates.DateTime)::Union{UncertainValue{Tuple{Float64,Float64}},Missing}
    impositions = get_impositions(tso_actions)
    if haskey(impositions, (gen_id, ts))
        return impositions[gen_id, ts]
    else
        # msg = @sprintf("impositions is not defined for (%s,%s)", gen_id, ts)
        # @warn msg
        return missing
    end
end

function get_imposition_uncertain_value!(tso_actions::TSOActions, gen_id::String, ts::Dates.DateTime)::Union{UncertainValue{Tuple{Float64,Float64}},Missing}
    impositions = get_impositions(tso_actions)
    if haskey(impositions, (gen_id, ts))
        return impositions[gen_id, ts]
    else
        impositions[gen_id, ts] = UncertainValue{Tuple{Float64,Float64}}()
        return impositions[gen_id, ts]
    end
end

function get_imposition(tso_actions::TSOActions, gen_id::String, ts::Dates.DateTime)::Union{Tuple{Float64,Float64},Missing}
    uncertain_value = get_imposition_uncertain_value(tso_actions, gen_id, ts)
    if ismissing(uncertain_value)
        return missing
    end
    return get_value(uncertain_value)
end

function safeget_imposition(tso_actions::TSOActions, gen_id::String, ts::Dates.DateTime)::Tuple{Float64,Float64}
    imposition = get_imposition(tso_actions, gen_id, ts)
    if ismissing(imposition)
        msg = @sprintf("Missing imposition value for gen_id=%s,ts=%s",
                        gen_id, ts)
        throw( error(msg) )
    else
        return imposition
    end
end

function get_imposition(tso_actions::TSOActions, gen_id::String, ts::Dates.DateTime, scenario::String)::Union{Tuple{Float64,Float64}, Missing}
    uncertain_value = get_imposition_uncertain_value(tso_actions, gen_id, ts)
    if ismissing(uncertain_value)
        return missing
    end
    return get_value(uncertain_value, scenario)
end

function safeget_imposition(tso_actions, gen_id::String, ts::Dates.DateTime, scenario::String)::Tuple{Float64,Float64}
    imposition = get_imposition(tso_actions, gen_id, ts, scenario)
    if ismissing(imposition)
        msg = @sprintf("No imposition entry in TSOActions for (gen_id=%s,ts=%s,s=%s).", gen_id, ts, scenario)
        throw(error(msg))
    else
        return imposition
    end
end

function get_imposition_level(tso_actions::TSOActions, gen_id::String, ts::Dates.DateTime, scenario::String)::Union{Float64, Missing}
    imposition = get_imposition(tso_actions, gen_id, ts, scenario)
    if ismissing(imposition)
        return missing
    else
        value_min, value_max = imposition
        if value_min != value_max
            msg = @sprintf("TSOActions for (gen_id=%s,ts=%s,s=%s) imposes interval [%s,%s] and not a single power level.",
                            gen_id, ts, scenario, value_min, value_max)
            throw(error(msg))
        else
            return value_min
        end
    end
end

function safeget_imposition_level(tso_actions::TSOActions, gen_id::String, ts::Dates.DateTime, scenario::String)::Float64
    imposition = get_imposition_level(tso_actions, gen_id, ts, scenario)
    if ismissing(imposition)
        msg = @sprintf("TSOActions has no imposition value for (gen_id=%s,ts=%s,s=%s).", gen_id, ts, scenario)
        throw(error(msg))
    else
        return imposition
    end
end

function set_imposition_definitive_value!(tso_actions::TSOActions, gen_id::String, ts::Dates.DateTime, value_min::Float64, value_max::Float64)
    uncertain_value = get_imposition_uncertain_value!(tso_actions, gen_id, ts)
    set_definitive_value!(uncertain_value, (value_min, value_max))
end

function set_imposition_value!(tso_actions::TSOActions, gen_id::String, ts::Dates.DateTime, scenario::String, value_min::Float64, value_max::Float64)
    uncertain_value = get_imposition_uncertain_value!(tso_actions, gen_id, ts)
    set_value!(uncertain_value, scenario, (value_min, value_max))
end

#######################
# Helpers
########################

"""
Performs a partial shallow copy of the input TSOActions
Modifying the returned tso_actions' kept attributes will modify the original one's.
"""
function filter_tso_actions(tso_actions::TSOActions;
                        keep_limitations::Bool=false,
                        keep_impositions::Bool=false)::TSOActions
    limitations_l = keep_limitations ? tso_actions.limitations : SortedDict{Tuple{String, Dates.DateTime}, Float64}()
    impositions_l = keep_impositions ? tso_actions.impositions : SortedDict{Tuple{String, Dates.DateTime}, Tuple{Float64,Float64}}()

    return TSOActions(limitations_l, impositions_l)
end

max_stars_repo_name: klorel/pscopf

max_issues_repo_name: klorel/pscopf

max_forks_repo_name: klorel/pscopf


text:
"Models used to test Omega"
module OmegaTestModels
using OmegaCore

# function t1()
#   μ = uniform(0.0, 1.0)
#   x = normal(μ, 1.0) + normal(μ, 1.0)
#   y = x ==ₛ 1.0
#   (name = :t1, vars = (μ = μ, x = x, y = y), isfinite = true, isdiff = true, cond = y)
# end

# function t2()
#   a = uniform(0.0, 1.0, 100)
#   b = normal(0.0, 1.0)
#   c = a * b
#   y = c >ₛ 0.0
#   (name = :t2, vars = (a = a, b = b, c = c, y = y), isfinite = true, isdiff = true, cond = y)
# end

# function t3()
#   n = uniform(1:10)
#   function x_(ω)
#     x = normal(ω, 0, 1)
#     for i = 1:n(ω)
#       x += normal(ω, 0, 1)
#     end
#     x
#   end
#   x = ciid(x_)
#   (name = :t3, vars = (x = x, n = n))
# end

# const allmodels = [
#   t1(),
#   t2(),
#   t3()]

# "`m` has key `k` and its true"
# hastrue(m, k) = haskey(m, k) && m[k]

# "Does `m` have a condition"
# hascond(m) = haskey(m, :cond)
# "is `m` differentiable"
# isdiff(m) = hastrue(m, :isdiff)
# "does `m` have a finite number of variables"
# isfinite(m) = hastrue(m, :isfinite)

# """Conjoin predicates
#   p = iseven ∧ (x -> x > 0) ∧ (x -> x < 100)
#   p(50)
# """
# function conjoin(preds::Vararg{Function})
#   function pred_conjunct(x...)
#     for pred in preds
#       if !pred(x...)
#         return false
#       end
#     end
#     true
#   end
#   pred_conjunct
# end
# # zt: add spec

# "`conjoin` ∧"
# const ∧ = conjoin

# export isfinite,
#        isdiff,
#        hascond,
#        allmodels,
#        ∧

end

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
using Singular
using Singular: mul!, addeq!, zero!

function shreq!(v :: Vector{spoly{A}}, u :: Vector{spoly{A}}; step :: Int = 1) where A
  (length(v) > 0 && length(u) > 0) || error("vector dimension equal to zero")
  R = parent(v[1])
  for ix in 1:length(u)-step
    v[ix] = u[ix+step]
  end
end

function _addeq!(l0 :: Vector{spoly{A}},
                 l1 :: Vector{spoly{A}},
                 l2 :: Vector{spoly{A}},
                 I  :: sideal{spoly{A}}) where A
  length(l0) == length(l1) == length(l2) || error("vector dimensions not equal")
  ln = length(l0)
  R = parent(l0[1])
  temp1 = zero(R)
  temp2 = zero(R)
  temp3 = zero(R)
  l0_ix_unred = zero(R)
  carry_unred = zero(R)
  carry = zero(R)
  for ix in 1:ln
    #w/o inlining -- l0[ix] = carry+l1[ix]+l2[ix]
    zero!(l0_ix_unred)
    addeq!(l0_ix_unred, l1[ix])
    addeq!(l0_ix_unred, l2[ix])
    addeq!(l0_ix_unred, carry)
    l0[ix] = reduce(l0_ix_unred, I)
    #note that xy ∨ xz ∨ yz == xy + xz + yz
    #
    #w/o inlining -- carry = carry*l1[ix]+carry*l2[ix]+l1[ix]*l2[ix]
    mul!(temp1, l1[ix], l2[ix])
    mul!(temp2, carry, l2[ix])
    mul!(temp3, carry, l1[ix])
    zero!(carry_unred)
    addeq!(carry_unred, temp1)
    addeq!(carry_unred, temp2)
    addeq!(carry_unred, temp3)
    zero!(carry)
    carry = reduce(carry_unred, I)
  end
end

function _addeq!(l0 :: Vector{spoly{A}},
                 l1 :: Vector{spoly{A}},
                 l2 :: Vector{spoly{A}}) where A
  length(l0) == length(l1) == length(l2) || error("vector dimensions not equal")
  ln = length(l0)
  R = parent(l0[1])
  temp1 = zero(R)
  temp2 = zero(R)
  temp3 = zero(R)
  l0_ix_unred = zero(R)
  carry_unred = zero(R)
  carry = zero(R)
  for ix in 1:ln
    #w/o inlining -- l0[ix] = carry+l1[ix]+l2[ix]
    addeq!(l0[ix], l1[ix])
    addeq!(l0[ix], l2[ix])
    addeq!(l0[ix], carry)
    #note that xy ∨ xz ∨ yz == xy + xz + yz
    #
    #w/o inlining -- carry = carry*l1[ix]+carry*l2[ix]+l1[ix]*l2[ix]
    mul!(temp1, l1[ix], l2[ix])
    mul!(temp2, carry, l2[ix])
    mul!(temp3, carry, l1[ix])
    zero!(carry)
    addeq!(carry, temp1)
    addeq!(carry, temp2)
    addeq!(carry, temp3)
  end
end

function subeq!(l0 :: Vector{spoly{A}},
                l1 :: Vector{spoly{A}},
                l2 :: Vector{spoly{A}},
                I  :: sideal{spoly{A}}) where A
  length(l0) == length(l1) == length(l2) || error("vector dimensions not equal")
  ln = length(l0)
  R = parent(l0[1])
  temp1 = zero(R)
  temp2 = zero(R)
  temp3 = zero(R)
  l0_ix_unred = zero(R)
  carry_unred = one(R)
  carry = one(R)
  for ix in 1:ln
    #w/o inlining -- l0[ix] = carry+l1[ix]+l2[ix]+one(R)
    zero!(l0_ix_unred)
    addeq!(l0_ix_unred, one(R))
    addeq!(l0_ix_unred, l2[ix])
    addeq!(l0_ix_unred, l1[ix])
    addeq!(l0_ix_unred, carry)
    l0[ix] = reduce(l0_ix_unred, I)
    #note that xy ∨ xz ∨ yz == xy + xz + yz
    #
    #w/o inlining -- carry = carry*l1[ix]+carry*l2[ix]+l1[ix]*l2[ix]+l1[ix]+carry
    mul!(temp1, l1[ix], l2[ix])
    mul!(temp2, carry, l2[ix])
    mul!(temp3, carry, l1[ix])
    zero!(carry_unred)
    addeq!(carry_unred, carry)
    addeq!(carry_unred, l1[ix])
    addeq!(carry_unred, temp1)
    addeq!(carry_unred, temp2)
    addeq!(carry_unred, temp3)
    zero!(carry)
    carry = reduce(carry_unred, I)
  end
end

function andeq!(l0 :: Vector{spoly{A}},
                l1 :: Vector{spoly{A}},
                l2 :: Vector{spoly{A}},
                I  :: sideal{spoly{A}}) where A
  length(l0) == length(l1) == length(l2) || error("vector dimensions not equal")
  ln = length(l0)
  R = parent(l0[1])
  l0_ix_unred = zero(R)
  for ix in 1:ln
    zero!(l0_ix_unred)
    mul!(l0_ix_unred,l1[ix],l2[ix])
    l0[ix] = reduce(l0_ix_unred, I)
  end
end

vars = ["b$ix" for ix in 0:31]
R, b = PolynomialRing(Fp(2), vars)
I = Ideal(R, [b[ix]^2-b[ix] for ix in 1:32])
I = std(I) #compute Gröbner basis of I
(b0,b1,b2,b3,b4,b5,b6,b7,
 b8,b9,b10,b11,b12,b13,b14,b15,
 b16,b17,b18,b19,b20,b21,b22,b23,
 b24,b25,b26,b27,b28,b29,b30,b31) = b
reduce(b0^2+b1^3, I)

_55555555 = [one(R),zero(R),one(R),zero(R),one(R),zero(R),one(R),zero(R),
             one(R),zero(R),one(R),zero(R),one(R),zero(R),one(R),zero(R),
             one(R),zero(R),one(R),zero(R),one(R),zero(R),one(R),zero(R),
             one(R),zero(R),one(R),zero(R),one(R),zero(R),one(R),zero(R)]

_33333333 = [one(R),one(R),zero(R),zero(R),one(R),one(R),zero(R),zero(R),
             one(R),one(R),zero(R),zero(R),one(R),one(R),zero(R),zero(R),
             one(R),one(R),zero(R),zero(R),one(R),one(R),zero(R),zero(R),
             one(R),one(R),zero(R),zero(R),one(R),one(R),zero(R),zero(R)]

_0f0f0f0f = [one(R),one(R),one(R),one(R),zero(R),zero(R),zero(R),zero(R),
             one(R),one(R),one(R),one(R),zero(R),zero(R),zero(R),zero(R),
             one(R),one(R),one(R),one(R),zero(R),zero(R),zero(R),zero(R),
             one(R),one(R),one(R),one(R),zero(R),zero(R),zero(R),zero(R)]

_01010101 = [one(R),zero(R),zero(R),zero(R),zero(R),zero(R),zero(R),zero(R),
             one(R),zero(R),zero(R),zero(R),zero(R),zero(R),zero(R),zero(R),
             one(R),zero(R),zero(R),zero(R),zero(R),zero(R),zero(R),zero(R),
             one(R),zero(R),zero(R),zero(R),zero(R),zero(R),zero(R),zero(R)]
#mov    edx,edi
#shr    edx,1
#and    edx,0x55555555
#sub    edi,edx
#mov    eax,edi
#shr    edi,0x2
#and    eax,0x33333333
#and    edi,0x33333333
#add    edi,eax
#mov    eax,edi
#shr    eax,0x4
#add    eax,edi
#and    eax,0x0f0f0f0f

edi1 = [b0,b1,b2,b3,b4,b5,b6,b7,
        b8,b9,b10,b11,b12,b13,b14,b15,
        b16,b17,b18,b19,b20,b21,b22,b23,
        b24,b25,b26,b27,b28,b29,b30,b31]
edx1 = [zero(R) for _ in 1:32]
shreq!(edx1, edi1);
edx2 = [zero(R) for _ in 1:32]
andeq!(edx2, edx1,_55555555, I);
edi2 = [zero(R) for _ in 1:32]
subeq!(edi2, edi1, edx2, I);
eax1 = deepcopy(edi2);
edi3 = [zero(R) for _ in 1:32]
shreq!(edi3, edi2; step=2);
eax2 = [zero(R) for _ in 1:32]
andeq!(eax2, eax1, _33333333, I)
edi4 = [zero(R) for _ in 1:32]
andeq!(edi4, edi3, _33333333, I)
edi5 = [zero(R) for _ in 1:32]
_addeq!(edi5, edi4, eax2, I)
eax3 = deepcopy(edi5)
eax4 = [zero(R) for _ in 1:32]
shreq!(eax4, eax3, step=4);
eax5 = [zero(R) for _ in 1:32]
_addeq!(eax5, eax4, edi5, I);
eax6 = [zero(R) for _ in 1:32]
andeq!(eax6, eax5, _0f0f0f0f, I);

#time for imul

vars1 = ["b0_7_1","b0_7_2","b0_7_4","b0_7_8",
         "b8_15_1","b8_15_2","b8_15_4","b8_15_8",
         "b16_23_1","b16_23_2","b16_23_4","b16_23_8",
         "b24_31_1","b24_31_2","b24_31_4","b24_31_8"]
R1, x = PolynomialRing(Fp(2), vars1)
I1 = Ideal(R1, [x[ix]^2-x[ix] for ix in 1:16])
I1 = std(I1)
(x0,x1,x2,x3,x4,x5,x6,x7,
 x8,x9,x10,x11,x12,x13,x14,x15) = x

(b0_7_1,b0_7_2,b0_7_4,b0_7_8,
 b8_15_1,b8_15_2,b8_15_4,b8_15_8,
 b16_23_1,b16_23_2,b16_23_4,b16_23_8,
 b24_31_1,b24_31_2,b24_31_4,b24_31_8) = x

v1 = [x0,x1,x2,x3,zero(R1),zero(R1),zero(R1),zero(R1),
      x4,x5,x6,x7,zero(R1),zero(R1),zero(R1),zero(R1),
      x8,x9,x10,x11,zero(R1),zero(R1),zero(R1),zero(R1),
      x12,x13,x14,x15,zero(R1),zero(R1),zero(R1),zero(R1)]
v2 = [zero(R1),zero(R1),zero(R1),zero(R1),zero(R1),zero(R1),zero(R1),zero(R1),
      x0,x1,x2,x3,zero(R1),zero(R1),zero(R1),zero(R1),
      x4,x5,x6,x7,zero(R1),zero(R1),zero(R1),zero(R1),
      x8,x9,x10,x11,zero(R1),zero(R1),zero(R1),zero(R1)]
v3 = [zero(R1),zero(R1),zero(R1),zero(R1),zero(R1),zero(R1),zero(R1),zero(R1),
      zero(R1),zero(R1),zero(R1),zero(R1),zero(R1),zero(R1),zero(R1),zero(R1),
      x0,x1,x2,x3,zero(R1),zero(R1),zero(R1),zero(R1),
      x4,x5,x6,x7,zero(R1),zero(R1),zero(R1),zero(R1)]
v4 = [zero(R1),zero(R1),zero(R1),zero(R1),zero(R1),zero(R1),zero(R1),zero(R1),
      zero(R1),zero(R1),zero(R1),zero(R1),zero(R1),zero(R1),zero(R1),zero(R1),
      zero(R1),zero(R1),zero(R1),zero(R1),zero(R1),zero(R1),zero(R1),zero(R1),
      x0,x1,x2,x3,zero(R1),zero(R1),zero(R1),zero(R1)]

#imul   eax,eax,0x01010101
#shr    eax,0x18
#ret

v5 = [zero(R1) for _ in 1:32]
_addeq!(v5,v1,v2,I1)
v6 = [zero(R1) for _ in 1:32]
_addeq!(v6,v5,v3,I1)
v7 = [zero(R1) for _ in 1:32]
_addeq!(v7,v6,v4,I1)
v8 = [zero(R1) for _ in 1:32]
shreq!(v8,v7;step=24)

@assert(all(v8[ix]==zero(R1) for ix in 7:32))

#note that
#
#b0_7_1 = all 1-factor products                 with factors in b0,b1,b2,...,b7
#b0_7_2 = all 2-factor products
#b0_7_1*b0_7_2 = all 3-factor products
#b0_7_4 = all 4-factor products
#b0_7_1*b0_7_4 = all 5-factor products
#b0_7_2*b0_7_4 = all 6-factor products
#b0_7_1*b0_7_2*b0_7_4 = all 7-factor products
#b0_7_8 = all 8-factor products

#hence v8[2] is sum of all 2-factor products    with factors in b0,b1,b2,...,b31
#      v8[3] is sum of all 4-factor products
#      v8[4] is sum of all 8-factor products
#      v8[5] is sum of all 16-factor products
#      v8[6] is sum of all 32-factor products


text:
include("ResultTests.jl")
include("TestFuncTests.jl")


text:
using RungeKutta: name, order, eachstage, nstages

@testset "$(rpad("Partitioned Tableau",80))" begin

    for s in 1:5
        for T ∈ (Float64, BigFloat)

            q = Tableau(:qTest, 2s, rand(T,s,s), rand(T,s), rand(T,s))
            p = Tableau(:pTest, 2s, rand(T,s,s), rand(T,s), rand(T,s))

            tab1 = PartitionedTableau{T}(:Test, 2s, s, q, p)
            tab2 = PartitionedTableau{T}(:Test, q, p)
            tab3 = PartitionedTableau(:Test, q, p)
            tab4 = PartitionedTableau(:Test, q)

            @test tab1 ≈ tab2 ≈ tab3
            @test tab1 == tab2 == tab3 != tab4
            @test hash(tab1) == hash(tab2) == hash(tab3) != hash(tab4)

            @test isequal(tab1, tab2)
            @test isequal(tab1, tab3)
            @test !isequal(tab1, tab4)

            @test eltype(tab1) == eltype(tab2) == eltype(tab3) == eltype(tab4) == T

            @test name(tab1) == :Test
            @test order(tab1) == 2s
            @test eachstage(tab1) == 1:s
            @test nstages(tab1) == s

        end

    end

end


text:
mutable struct WishartStatXMLElement <: MyXMLElement
    el::XMLOrNothing
    traitLikelihood_el::TraitLikelihoodXMLElement
    rm_el::Union{Nothing, RepeatedMeasuresXMLElement}
    mbd_el::Union{Nothing, MBDXMLElement}

    function WishartStatXMLElement(like_el::TraitLikelihoodXMLElement,
                mbd_el::MBDXMLElement)
        return new(nothing, like_el, nothing, mbd_el)
    end

    function WishartStatXMLElement(like_el::TraitLikelihoodXMLElement,
                rm_el::RepeatedMeasuresXMLElement)
        return new(nothing, like_el, rm_el, nothing)
    end
end

mutable struct PrecisionGibbsOperatorXMLElement <: OperatorXMLElement
    el::XMLOrNothing
    wishart_stats::WishartStatXMLElement
    weight::Float64

    function PrecisionGibbsOperatorXMLElement(
            traitLikelihood_el::TraitLikelihoodXMLElement,
            mbd_el::MBDXMLElement)

        wishart_stats = WishartStatXMLElement(traitLikelihood_el, mbd_el)
        return new(nothing, wishart_stats, 1.0)
    end

    function PrecisionGibbsOperatorXMLElement(
        traitLikelihood_el::TraitLikelihoodXMLElement,
        rm_el::RepeatedMeasuresXMLElement)

        wishart_stats = WishartStatXMLElement(traitLikelihood_el, rm_el)
        return new(nothing, wishart_stats, 1.0)
    end

end

mutable struct CompoundPrecisionOperatorXMLElement <: OperatorXMLElement
    el::XMLOrNothing
    diff_el::PrecisionGibbsOperatorXMLElement
    res_el::PrecisionGibbsOperatorXMLElement
    weight::Float64

    function CompoundPrecisionOperatorXMLElement(
            diff_el::PrecisionGibbsOperatorXMLElement,
            res_el::PrecisionGibbsOperatorXMLElement)

        diff_ws = diff_el.wishart_stats
        res_ws = res_el.wishart_stats
        @assert !isnothing(diff_ws.mbd_el) && isnothing(diff_ws.rm_el)
        @assert isnothing(res_ws.mbd_el) && !isnothing(res_ws.rm_el)

        return new(nothing, diff_el, res_el, 1.0)
    end

    function CompoundPrecisionOperatorXMLElement(
                tl_el::TraitLikelihoodXMLElement,
                mbd_el::MBDXMLElement,
                rm_el::RepeatedMeasuresXMLElement)

        diff_el = PrecisionGibbsOperatorXMLElement(tl_el, mbd_el)
        res_el = PrecisionGibbsOperatorXMLElement(tl_el, rm_el)
        return CompoundPrecisionOperatorXMLElement(diff_el, res_el)
    end

end

function make_xml(cpo_el::CompoundPrecisionOperatorXMLElement)
    el = new_element(bn.COMPOUND_PRECISION_OPERATOR)
    set_attribute(el, bn.WEIGHT, cpo_el.weight)

    do_el = new_child(el, bn.DIFFUSION_OPERATOR)
    make_xml(cpo_el.diff_el)
    add_child(do_el, cpo_el.diff_el.el)

    ro_el = new_child(el, bn.RESIDUAL_OPERATOR)
    make_xml(cpo_el.res_el)
    add_child(ro_el, cpo_el.res_el.el)

    cpo_el.el = el
    return el
end

function make_xml(ws_el::WishartStatXMLElement)
    make_xml(ws_el.traitLikelihood_el)

    if isnothing(ws_el.rm_el) && !isnothing(ws_el.mbd_el)
        el = new_element(bn.WISHART_STATISTICS)
        return finish_Wishart_el(ws_el, el)
    elseif !isnothing(ws_el.rm_el) && isnothing(ws_el.mbd_el)
        make_xml(ws_el.rm_el)
        el = new_element(bn.REPEATED_MEASURES_WISHART_STATISTICS)
        add_ref_el(el, ws_el.rm_el.el)
        return finish_Wishart_el(ws_el, el)
    else
        error("Not implemented")
    end
end

function finish_Wishart_el(ws_el::WishartStatXMLElement, el::XMLElement)
    add_ref_el(el, ws_el.traitLikelihood_el.el)
    set_attribute(el, bn.TRAIT_NAME, ws_el.traitLikelihood_el.attrs[bn.TRAIT_NAME])
    ws_el.el = el
    return el
end

function make_xml(pgo_el::PrecisionGibbsOperatorXMLElement)
    el = new_element(bn.PRECISION_GIBBS_OPERATOR)
    set_attribute(el, bn.WEIGHT, pgo_el.weight)
    ws_el = pgo_el.wishart_stats
    make_xml(ws_el)
    add_child(el, ws_el.el)
    if isnothing(ws_el.rm_el) && !isnothing(ws_el.mbd_el)
        make_xml(ws_el.mbd_el)
        add_ref_el(el, ws_el.mbd_el.precision_prior.el)
    elseif !isnothing(ws_el.rm_el) && isnothing(ws_el.mbd_el)
        make_xml(ws_el.rm_el)
        add_ref_el(el, ws_el.rm_el.precision_prior.el)
    else
        error("Not implemented")
    end
    pgo_el.el = el
    return el
end


text:
# These functions are only conditionally loaded with Clustering.jl
# Code adapted from @cormullion's [ColorSchemeTools](https://github.com/JuliaGraphics/ColorSchemeTools.jl).
function get_colorscheme(
    img,
    ncolors;
    maxiter=Clustering._kmeans_default_maxiter,
    tol=Clustering._kmeans_default_tol,
)
    # Cluster in Lab color space
    data = reshape(channelview(Lab.(img)), 3, :)
    R = Clustering.kmeans(data, ncolors; maxiter=maxiter, tol=tol)

    # Make color scheme out of cluster centers
    cs = Lab{Float64}[]
    for i in 1:3:length(R.centers)
        push!(cs, Lab(R.centers[i], R.centers[i + 1], R.centers[i + 2]))
    end
    return cs
end

function _colordither(
    ::Type{T},
    img,
    alg,
    ncolors::Int;
    maxiter=Clustering._kmeans_default_maxiter,
    tol=Clustering._kmeans_default_tol,
    kwargs...,
) where {T}
    cs = get_colorscheme(img, ncolors; maxiter=maxiter, tol=tol)
    return _colordither(T, img, alg, cs; kwargs...)
end

"""
    dither!([out,] img, alg::AbstractDither, ncolors; maxiter, tol, kwargs...)

Dither image `img` using algorithm `alg`.
A color palette with `ncolors` is computed by Clustering.jl's K-means clustering.
The amount of `maxiter` and tolerance `tol` default to those exported by Clustering.jl.
"""
dither!(img, alg::AbstractDither, ncolors::Int; kwargs...)

"""
    dither([T::Type,] img, alg::AbstractDither, ncolors; maxiter, tol, kwargs...)

Dither image `img` using algorithm `alg`.
A color palette with `ncolors` is computed by Clustering.jl's K-means clustering.
The amount of `maxiter` and tolerance `tol` default to those exported by Clustering.jl.
"""
dither(::Type, img, alg::AbstractDither, ncolors::Int; kwargs...)


text:
import RayTracing: AzimuthalQuadrature, nazim2

"""
    Quadrature{A<:AzimuthalQuadrature,P<:PolarQuadrature,T<:Real}

Holds information regarding both the azimuthal and polar quadrature and the total weights.
"""
struct Quadrature{T<:Real,A<:AzimuthalQuadrature,P<:PolarQuadrature}
    azimuthal::A
    polar::P
    ω::Matrix{T}
end

function Quadrature(azimuthal::AzimuthalQuadrature{Na,T}, polar::PolarQuadrature{Np,T}) where {Na,Np,T}
    @unpack δs, ωₐ = azimuthal
    @unpack sinθs, ωₚ = polar
    n_azim_2 = nazim2(azimuthal)
    n_polar_2 = npolar2(polar)

    # IDEA: I think we can use n_azim_4 because the matrix has repeated values
    ω = Matrix{T}(undef, n_azim_2, n_polar_2)

    for i in 1:n_azim_2, j in 1:n_polar_2
        ω[i, j] = 4π * ωₐ[i] * ωₚ[j] * δs[i] * sinθs[j]
    end
    # ω .*= 4π

    return Quadrature(azimuthal, polar, ω)
end