text:
"""
    buoyancy_gradients(
        param_set,
        bg_model::EnvBuoyGrad{FT, EBG}
    ) where {FT <: Real, EBG <: AbstractEnvBuoyGradClosure}

Returns the vertical buoyancy gradients in the environment, as well as in its dry and cloudy volume fractions.
The dispatch on EnvBuoyGrad type is performed at the EnvBuoyGrad construction time, and the analytical solutions
used here are consistent for both mean fields and conditional fields obtained from assumed distributions
over the conserved thermodynamic variables.
"""
function buoyancy_gradients(
    param_set,
    bg_model::EnvBuoyGrad{FT, EBG},
) where {FT <: Real, EBG <: AbstractEnvBuoyGradClosure}

    g = CPP.grav(param_set)
    molmass_ratio = CPP.molmass_ratio(param_set)
    R_d = CPP.R_d(param_set)
    R_v = CPP.R_v(param_set)

    phase_part = TD.PhasePartition(0.0, 0.0, 0.0) # assuming R_d = R_m
    Π = TD.exner_given_pressure(param_set, bg_model.p0, phase_part)

    ∂b∂θv = g * (R_d / bg_model.alpha0 / bg_model.p0) * Π

    if bg_model.en_cld_frac > 0.0
        ts_sat = thermo_state_pθq(param_set, bg_model.p0, bg_model.θ_liq_ice_sat, bg_model.qt_sat)
        phase_part = TD.PhasePartition(ts_sat)
        lh = TD.latent_heat_liq_ice(param_set, phase_part)
        cp_m = TD.cp_m(ts_sat)
        ∂b∂θl_sat = (
            ∂b∂θv * (1 + molmass_ratio * (1 + lh / R_v / bg_model.t_sat) * bg_model.qv_sat - bg_model.qt_sat) /
            (1 + lh * lh / cp_m / R_v / bg_model.t_sat / bg_model.t_sat * bg_model.qv_sat)
        )
        ∂b∂qt_sat = (lh / cp_m / bg_model.t_sat * ∂b∂θl_sat - ∂b∂θv) * bg_model.θ_sat
    else
        ∂b∂θl_sat = FT(0)
        ∂b∂qt_sat = FT(0)
    end

    ∂b∂z, ∂b∂z_unsat, ∂b∂z_sat = buoyancy_gradient_chain_rule(bg_model, ∂b∂θv, ∂b∂θl_sat, ∂b∂qt_sat)
    return GradBuoy(∂b∂z, ∂b∂z_unsat, ∂b∂z_sat)
end

"""
    buoyancy_gradient_chain_rule(
        bg_model::EnvBuoyGrad{FT, EBG},
        ∂b∂θv::FT,
        ∂b∂θl_sat::FT,
        ∂b∂qt_sat::FT,
    ) where {FT <: Real, EBG <: AbstractEnvBuoyGradClosure}

Returns the vertical buoyancy gradients in the environment, as well as in its dry and cloudy volume fractions,
from the partial derivatives with respect to thermodynamic variables in dry and cloudy volumes.
"""
function buoyancy_gradient_chain_rule(
    bg_model::EnvBuoyGrad{FT, EBG},
    ∂b∂θv::FT,
    ∂b∂θl_sat::FT,
    ∂b∂qt_sat::FT,
) where {FT <: Real, EBG <: AbstractEnvBuoyGradClosure}
    if bg_model.en_cld_frac > FT(0)
        ∂b∂z_θl_sat = ∂b∂θl_sat * bg_model.∂θl∂z_sat
        ∂b∂z_qt_sat = ∂b∂qt_sat * bg_model.∂qt∂z_sat
    else
        ∂b∂z_θl_sat = FT(0)
        ∂b∂z_qt_sat = FT(0)
    end

    ∂b∂z_unsat = bg_model.en_cld_frac < FT(1) ? ∂b∂θv * bg_model.∂θv∂z_unsat : FT(0)

    ∂b∂z_sat = ∂b∂z_θl_sat + ∂b∂z_qt_sat
    ∂b∂z = (1 - bg_model.en_cld_frac) * ∂b∂z_unsat + bg_model.en_cld_frac * ∂b∂z_sat

    return ∂b∂z, ∂b∂z_unsat, ∂b∂z_sat
end


text:
module WorkspaceManager
import UUIDs: UUID
import ..Pluto: Notebook, Cell, PKG_ROOT_DIR, ExpressionExplorer, pluto_filename, trycatch_expr, Token, withtoken, get_pl_env, tamepath
import ..PlutoRunner
import Distributed

"Contains the Julia process (in the sense of `Distributed.addprocs`) to evaluate code in. Each notebook gets at most one `Workspace` at any time, but it can also have no `Workspace` (it cannot `eval` code in this case)."
mutable struct Workspace
    pid::Integer
    module_name::Symbol
    dowork_token::Token
end

Workspace(pid::Integer, module_name::Symbol) = Workspace(pid, module_name, Token())

"These expressions get evaluated inside every newly create module inside a `Workspace`."
const workspace_preamble = [
    :(using Markdown, InteractiveUtils, Main.PlutoRunner),
    :(ENV["GKSwstype"] = "nul"),
    :(show, showable, showerror, repr, string, print, println), # https://github.com/JuliaLang/julia/issues/18181
]

"These expressions get evaluated whenever a new `Workspace` process is created."
const process_preamble = [
    :(ccall(:jl_exit_on_sigint, Cvoid, (Cint,), 0)),
    :(include($(joinpath(PKG_ROOT_DIR, "src", "runner", "PlutoRunner.jl")))),
]

const moduleworkspace_count = Ref(0)
const workspaces = Dict{UUID,Workspace}()

"""Create a workspace for the notebook, optionally in a separate process.

`new_process`: Should future workspaces be created on a separate process (`true`) or on the same one (`false`)? Only workspaces on a separate process can be stopped during execution. Windows currently supports `true` only partially: you can't stop cells on Windows. _Defaults to `get_pl_env("PLUTO_WORKSPACE_USE_DISTRIBUTED")`_"""
function make_workspace(notebook::Notebook, new_process=(get_pl_env("PLUTO_WORKSPACE_USE_DISTRIBUTED") == "true"))::Workspace
    pid = if new_process
        create_workspaceprocess(;environment_path=resolved_environment_path(notebook))
    else
        pid = Distributed.myid()
        # for some reason the PlutoRunner might not be available in Main unless we include the file
        # (even though this is the main process)
        if !Distributed.remotecall_eval(Main, pid, :(isdefined(Main, :PlutoRunner) && PlutoRunner isa Module))
            for expr in process_preamble
                Distributed.remotecall_eval(Main, [pid], expr)
            end
        end
        pid
    end

    module_name = create_emptyworkspacemodule(pid)
    workspace = Workspace(pid, module_name)

    cd_workspace(workspace, notebook.path)

    return workspace
end

"Call `cd(\$path)` inside the workspace. This is done when creating a workspace, and whenever the notebook changes path."
function cd_workspace(workspace, path::AbstractString)
    eval_in_workspace(workspace, quote
        cd($(path |> dirname))
        Sys.set_process_title("Pluto - " * $(path |> basename))
    end)
end

# TODO: move to PlutoRunner
function create_emptyworkspacemodule(pid::Integer)::Symbol
    id = (moduleworkspace_count[] += 1)

    new_workspace_name = if Distributed.myid() == 1
        Symbol("workspace", id)
    else
        Symbol("workspace", id, "_", Distributed.myid())
    end
    workspace_creation = :(module $(new_workspace_name) $(workspace_preamble...) end)

    Distributed.remotecall_eval(Main, [pid], workspace_creation)
    Distributed.remotecall_eval(Main, [pid], :(PlutoRunner.set_current_module($(new_workspace_name |> QuoteNode))))

    new_workspace_name
end

function default_environment_path(::Nothing)
    get_pl_env("PLUTO_DEFAULT_ENVIRONMENT_PATH")
end
default_environment_path(s) = s

"""
    resolved_environment_path(notebook)

Return the packge environment path of given notebook.
"""
function resolved_environment_path(notebook::Notebook)
    if notebook.environment_path === nothing
        nothing
    else
        if isabspath(notebook.environment_path)
            tamepath(notebook.environment_path)
        else
            # relative project path is always relative to
            # the notebook path.
            tamepath(joinpath(dirname(notebook.path), notebook.environment_path))
        end
    end
end

function create_workspaceprocess(;environment_path::Union{String, Nothing}=nothing)::Integer
    pid = Distributed.addprocs(1; exeflags="--project=$(default_environment_path(environment_path))") |> first

    for expr in process_preamble
        Distributed.remotecall_eval(Main, [pid], expr)
    end

    # so that we NEVER break the workspace with an interrupt 🤕
    @async Distributed.remotecall_eval(Main, [pid],
        :(while true
            try
                wait()
            catch end
        end))

    pid
end

"Return the `Workspace` of `notebook`; will be created if none exists yet."
function get_workspace(notebook::Notebook)::Workspace
    if haskey(workspaces, notebook.notebook_id)
        workspaces[notebook.notebook_id]
    else
        workspaces[notebook.notebook_id] = make_workspace(notebook)
    end
end
get_workspace(workspace::Workspace)::Workspace = workspace

"Try our best to delete the workspace. `ProcessWorkspace` will have its worker process terminated."
function unmake_workspace(notebook::Union{Notebook,Workspace})
    workspace = get_workspace(notebook)

    if workspace.pid != Distributed.myid()
        filter!(p -> p.second.pid != workspace.pid, workspaces)
        interrupt_workspace(workspace; verbose=false)
        if workspace.pid != Distributed.myid()
            Distributed.rmprocs(workspace.pid)
        end
    end
end

"Evaluate expression inside the workspace - output is fetched and formatted, errors are caught and formatted. Returns formatted output and error flags.

`expr` has to satisfy `ExpressionExplorer.is_toplevel_expr`."
function eval_format_fetch_in_workspace(notebook::Union{Notebook,Workspace}, expr::Expr, cell_id::UUID, ends_with_semicolon::Bool=false)::NamedTuple{(:output_formatted, :errored, :interrupted, :runtime),Tuple{PlutoRunner.MimedOutput,Bool,Bool,Union{UInt64,Missing}}}
    workspace = get_workspace(notebook)

    # if multiple notebooks run on the same process, then we need to `cd` between the different notebook paths
    if workspace.pid == Distributed.myid() && notebook isa Notebook
        cd_workspace(workspace, notebook.path)
    end

    # We wrap the expression in a try-catch block, because we want to capture and format the exception on the worker itself.
    wrapped = trycatch_expr(expr, workspace.module_name, cell_id)

    # run the code 🏃♀

    # a try block (on this process) to catch an InterruptException
    take!(workspace.dowork_token)
    try
        # we use [pid] instead of pid to prevent fetching output
        Distributed.remotecall_eval(Main, [workspace.pid], wrapped)
        put!(workspace.dowork_token)
    catch exs
        # We don't use a `finally` because the token needs to be back asap
        put!(workspace.dowork_token)
        try
            @assert exs isa CompositeException
            ex = exs.exceptions |> first
            @assert ex isa Distributed.RemoteException
            @assert ex.pid == workspace.pid
            @assert ex.captured.ex isa InterruptException

            return (output_formatted = PlutoRunner.format_output(CapturedException(InterruptException(), [])), errored = true, interrupted = true, runtime = missing)
        catch assertionerr
            showerror(stderr, exs)
            return (output_formatted = PlutoRunner.format_output(CapturedException(exs, [])), errored = true, interrupted = true, runtime = missing)
        end
    end

    # instead of fetching the output value (which might not make sense in our context, since the user can define structs, types, functions, etc), we format the cell output on the worker, and fetch the formatted output.
    # This also means that very big objects are not duplicated in RAM.
    withtoken(workspace.dowork_token) do
        Distributed.remotecall_eval(Main, workspace.pid, :(PlutoRunner.formatted_result_of($cell_id, $ends_with_semicolon)))
    end
end

"Evaluate expression inside the workspace - output is not fetched, errors are rethrown. For internal use."
function eval_in_workspace(notebook::Union{Notebook,Workspace}, expr)
    workspace = get_workspace(notebook)

    Distributed.remotecall_eval(Main, [workspace.pid], :(Core.eval($(workspace.module_name), $(expr |> QuoteNode))))
    nothing
end

"Evaluate expression inside the workspace - output is returned. For internal use."
function eval_fetch_in_workspace(notebook::Union{Notebook,Workspace}, expr)
    workspace = get_workspace(notebook)

    Distributed.remotecall_eval(Main, workspace.pid, :(Core.eval($(workspace.module_name), $(expr |> QuoteNode))))
end

"Fake deleting variables by moving to a new module without re-importing them."
function delete_vars(notebook::Union{Notebook,Workspace}, to_delete::Set{Symbol}, funcs_to_delete::Set{Vector{Symbol}}, module_imports_to_move::Set{Expr}; kwargs...)
    workspace = get_workspace(notebook)

    old_workspace_name = workspace.module_name
    new_workspace_name = create_emptyworkspacemodule(workspace.pid)

    workspace.module_name = new_workspace_name
    Distributed.remotecall_eval(Main, [workspace.pid], :(PlutoRunner.set_current_module($(new_workspace_name |> QuoteNode))))

    Distributed.remotecall_eval(Main, [workspace.pid], :(PlutoRunner.move_vars($(old_workspace_name |> QuoteNode), $(new_workspace_name |> QuoteNode), $to_delete, $funcs_to_delete, $module_imports_to_move)))
end

"Force interrupt (SIGINT) a workspace, return whether succesful"
function interrupt_workspace(notebook::Union{Notebook,Workspace}; verbose=true)::Bool
    workspace = get_workspace(notebook)

    if Sys.iswindows()
        verbose && @warn "Unfortunately, stopping cells is currently not supported on Windows :(
        Maybe the Windows Subsystem for Linux is right for you:
        https://docs.microsoft.com/en-us/windows/wsl"
        return false
    end
    if workspace.pid == Distributed.myid()
        verbose && @warn """Cells in this workspace can't be stopped, because it is not running in a separate workspace. Use `ENV["PLUTO_WORKSPACE_USE_DISTRIBUTED"]` to control whether future workspaces are generated in a separate process."""
        return false
    end
    if isready(workspace.dowork_token)
        verbose && @info "Tried to stop idle workspace - ignoring."
        return true
    end

    # You can force kill a julia process by pressing Ctrl+C five times 🙃
    # But this is not very consistent, so we will just keep pressing Ctrl+C until the workspace isn't running anymore.
    # TODO: this will also kill "pending" evaluations, and any evaluations started within 100ms of the kill. A global "evaluation count" would fix this.
    # TODO: listen for the final words of the remote process on stdout/stderr: "Force throwing a SIGINT"
    try
        verbose && @info "Sending interrupt to process $(workspace.pid)"
        Distributed.interrupt(workspace.pid)

        delay = 5.0 # seconds
        parts = 100

        for _ in 1:parts
            sleep(delay / parts)
            if isready(workspace.dowork_token)
                verbose && println("Cell interrupted!")
                return true
            end
        end

        verbose && println("Still running... starting sequence")
        while !isready(workspace.dowork_token)
            for _ in 1:5
                verbose && print(" 🔥 ")
                Distributed.interrupt(workspace.pid)
                sleep(0.2)
            end
            sleep(1.5)
        end
        verbose && println()
        verbose && println("Cell interrupted!")
        true
    catch ex
        if !(ex isa KeyError)
            @warn "Interrupt failed for unknown reason"
            showerror(ex, stacktrace(catch_backtrace()))
        end
        false
    end
end

end

ext: jl


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

import Base: eltype

abstract type AbstractRemoteRef end

"""
    client_refs

Tracks whether a particular `AbstractRemoteRef`
(identified by its RRID) exists on this worker.

The `client_refs` lock is also used to synchronize access to `.refs` and associated `clientset` state.
"""
const client_refs = WeakKeyDict{AbstractRemoteRef, Nothing}() # used as a WeakKeySet

"""
    Future(w::Int, rrid::RRID, v::Union{Some, Nothing}=nothing)

A `Future` is a placeholder for a single computation
of unknown termination status and time.
For multiple potential computations, see `RemoteChannel`.
See `remoteref_id` for identifying an `AbstractRemoteRef`.
"""
mutable struct Future <: AbstractRemoteRef
    where::Int
    whence::Int
    id::Int
    lock::ReentrantLock
    @atomic v::Union{Some{Any}, Nothing}

    Future(w::Int, rrid::RRID, v::Union{Some, Nothing}=nothing) =
        (r = new(w,rrid.whence,rrid.id,ReentrantLock(),v); return test_existing_ref(r))

    Future(t::NTuple{4, Any}) = new(t[1],t[2],t[3],ReentrantLock(),t[4])  # Useful for creating dummy, zeroed-out instances
end

"""
    RemoteChannel(pid::Integer=myid())

Make a reference to a `Channel{Any}(1)` on process `pid`.
The default `pid` is the current process.

    RemoteChannel(f::Function, pid::Integer=myid())

Create references to remote channels of a specific size and type. `f` is a function that
when executed on `pid` must return an implementation of an `AbstractChannel`.

For example, `RemoteChannel(()->Channel{Int}(10), pid)`, will return a reference to a
channel of type `Int` and size 10 on `pid`.

The default `pid` is the current process.
"""
mutable struct RemoteChannel{T<:AbstractChannel} <: AbstractRemoteRef
    where::Int
    whence::Int
    id::Int

    function RemoteChannel{T}(w::Int, rrid::RRID) where T<:AbstractChannel
        r = new(w, rrid.whence, rrid.id)
        return test_existing_ref(r)
    end

    function RemoteChannel{T}(t::Tuple) where T<:AbstractChannel
        return new(t[1],t[2],t[3])
    end
end

function test_existing_ref(r::AbstractRemoteRef)
    found = getkey(client_refs, r, nothing)
    if found !== nothing
        @assert r.where > 0
        if isa(r, Future)
            # this is only for copying the reference from Future to RemoteRef (just created)
            fv_cache = @atomic :acquire found.v
            rv_cache = @atomic :monotonic r.v
            if fv_cache === nothing && rv_cache !== nothing
                # we have recd the value from another source, probably a deserialized ref, send a del_client message
                send_del_client(r)
                @lock found.lock begin
                    @atomicreplace found.v nothing => rv_cache
                end
            end
        end
        return found::typeof(r)
    end

    client_refs[r] = nothing
    finalizer(finalize_ref, r)
    return r
end

function finalize_ref(r::AbstractRemoteRef)
    if r.where > 0 # Handle the case of the finalizer having been called manually
        if trylock(client_refs.lock) # trylock doesn't call wait which causes yields
            try
                delete!(client_refs.ht, r) # direct removal avoiding locks
                if isa(r, RemoteChannel)
                    send_del_client_no_lock(r)
                else
                    # send_del_client only if the reference has not been set
                    v_cache = @atomic :monotonic r.v
                    v_cache === nothing && send_del_client_no_lock(r)
                    @atomic :monotonic r.v = nothing
                end
                r.where = 0
            finally
                unlock(client_refs.lock)
            end
        else
            finalizer(finalize_ref, r)
            return nothing
        end
    end
    nothing
end

"""
    Future(pid::Integer=myid())

Create a `Future` on process `pid`.
The default `pid` is the current process.
"""
Future(pid::Integer=myid()) = Future(pid, RRID())
Future(w::LocalProcess) = Future(w.id)
Future(w::Worker) = Future(w.id)

RemoteChannel(pid::Integer=myid()) = RemoteChannel{Channel{Any}}(pid, RRID())

function RemoteChannel(f::Function, pid::Integer=myid())
    remotecall_fetch(pid, f, RRID()) do f, rrid
        rv=lookup_ref(rrid, f)
        RemoteChannel{typeof(rv.c)}(myid(), rrid)
    end
end

Base.eltype(::Type{RemoteChannel{T}}) where {T} = eltype(T)

hash(r::AbstractRemoteRef, h::UInt) = hash(r.whence, hash(r.id, h))
==(r::AbstractRemoteRef, s::AbstractRemoteRef) = (r.whence==s.whence && r.id==s.id)

"""
    remoteref_id(r::AbstractRemoteRef) -> RRID

`Future`s and `RemoteChannel`s are identified by fields:

* `where` - refers to the node where the underlying object/storage
  referred to by the reference actually exists.

* `whence` - refers to the node the remote reference was created from.
  Note that this is different from the node where the underlying object
  referred to actually exists. For example calling `RemoteChannel(2)`
  from the master process would result in a `where` value of 2 and
  a `whence` value of 1.

* `id` is unique across all references created from the worker specified by `whence`.

Taken together,  `whence` and `id` uniquely identify a reference across all workers.

`remoteref_id` is a low-level API which returns a `RRID`
object that wraps `whence` and `id` values of a remote reference.
"""
remoteref_id(r::AbstractRemoteRef) = RRID(r.whence, r.id)

"""
    channel_from_id(id) -> c

A low-level API which returns the backing `AbstractChannel` for an `id` returned by
[`remoteref_id`](@ref).
The call is valid only on the node where the backing channel exists.
"""
function channel_from_id(id)
    rv = lock(client_refs) do
        return get(PGRP.refs, id, false)
    end
    if rv === false
        throw(ErrorException("Local instance of remote reference not found"))
    end
    return rv.c
end

lookup_ref(rrid::RRID, f=def_rv_channel) = lookup_ref(PGRP, rrid, f)
function lookup_ref(pg, rrid, f)
    return lock(client_refs) do
        rv = get(pg.refs, rrid, false)
        if rv === false
            # first we've heard of this ref
            rv = RemoteValue(invokelatest(f))
            pg.refs[rrid] = rv
            push!(rv.clientset, rrid.whence)
        end
        return rv
    end::RemoteValue
end

"""
    isready(rr::Future)

Determine whether a [`Future`](@ref) has a value stored to it.

If the argument `Future` is owned by a different node, this call will block to wait for the answer.
It is recommended to wait for `rr` in a separate task instead
or to use a local [`Channel`](@ref) as a proxy:

```julia
p = 1
f = Future(p)
errormonitor(@async put!(f, remotecall_fetch(long_computation, p)))
isready(f)  # will not block
```
"""
function isready(rr::Future)
    v_cache = @atomic rr.v
    v_cache === nothing || return true

    rid = remoteref_id(rr)
    return if rr.where == myid()
        isready(lookup_ref(rid).c)
    else
        remotecall_fetch(rid->isready(lookup_ref(rid).c), rr.where, rid)
    end
end

"""
    isready(rr::RemoteChannel, args...)

Determine whether a [`RemoteChannel`](@ref) has a value stored to it.
Note that this function can cause race conditions, since by the
time you receive its result it may no longer be true. However,
it can be safely used on a [`Future`](@ref) since they are assigned only once.
"""
function isready(rr::RemoteChannel, args...)
    rid = remoteref_id(rr)
    return if rr.where == myid()
        isready(lookup_ref(rid).c, args...)
    else
        remotecall_fetch(rid->isready(lookup_ref(rid).c, args...), rr.where, rid)
    end
end

del_client(rr::AbstractRemoteRef) = del_client(remoteref_id(rr), myid())

del_client(id, client) = del_client(PGRP, id, client)
function del_client(pg, id, client)
    lock(client_refs) do
        _del_client(pg, id, client)
    end
    nothing
end

function _del_client(pg, id, client)
    rv = get(pg.refs, id, false)
    if rv !== false
        delete!(rv.clientset, client)
        if isempty(rv.clientset)
            delete!(pg.refs, id)
            #print("$(myid()) collected $id\n")
        end
    end
    nothing
end

function del_clients(pairs::Vector)
    for p in pairs
        del_client(p[1], p[2])
    end
end

# The task below is coalescing the `flush_gc_msgs` call
# across multiple producers, see `send_del_client`,
# and `send_add_client`.
# XXX: Is this worth the additional complexity?
#      `flush_gc_msgs` has to iterate over all connected workers.
const any_gc_flag = Threads.Condition()
function start_gc_msgs_task()
    errormonitor(
        Threads.@spawn begin
            while true
                lock(any_gc_flag) do
                    # this might miss events
                    wait(any_gc_flag)
                end
                # Use invokelatest() so that custom message transport streams
                # for workers can be defined in a newer world age than the Task
                # which runs the loop here.
                invokelatest(flush_gc_msgs) # handles throws internally
            end
        end
    )
end

# Function can be called within a finalizer
function send_del_client(rr)
    if rr.where == myid()
        del_client(rr)
    elseif id_in_procs(rr.where) # process only if a valid worker
        process_worker(rr)
    end
end

function send_del_client_no_lock(rr)
    # for gc context to avoid yields
    if rr.where == myid()
        _del_client(PGRP, remoteref_id(rr), myid())
    elseif id_in_procs(rr.where) # process only if a valid worker
        process_worker(rr)
    end
end

function publish_del_msg!(w::Worker, msg)
    lock(w.msg_lock) do
        push!(w.del_msgs, msg)
        @atomic w.gcflag = true
    end
    lock(any_gc_flag) do
        notify(any_gc_flag)
    end
end

function process_worker(rr)
    w = worker_from_id(rr.where)::Worker
    msg = (remoteref_id(rr), myid())

    # Needs to aquire a lock on the del_msg queue
    T = Threads.@spawn begin
        publish_del_msg!($w, $msg)
    end
    Base.errormonitor(T)

    return
end

function add_client(id, client)
    lock(client_refs) do
        rv = lookup_ref(id)
        push!(rv.clientset, client)
    end
    nothing
end

function add_clients(pairs::Vector)
    for p in pairs
        add_client(p[1], p[2]...)
    end
end

function send_add_client(rr::AbstractRemoteRef, i)
    if rr.where == myid()
        add_client(remoteref_id(rr), i)
    elseif (i != rr.where) && id_in_procs(rr.where)
        # don't need to send add_client if the message is already going
        # to the processor that owns the remote ref. it will add_client
        # itself inside deserialize().
        w = worker_from_id(rr.where)
        lock(w.msg_lock) do
            push!(w.add_msgs, (remoteref_id(rr), i))
            @atomic w.gcflag = true
        end
        lock(any_gc_flag) do
            notify(any_gc_flag)
        end
    end
end

channel_type(rr::RemoteChannel{T}) where {T} = T

function serialize(s::ClusterSerializer, f::Future)
    v_cache = @atomic f.v
    if v_cache === nothing
        p = worker_id_from_socket(s.io)
        (p !== f.where) && send_add_client(f, p)
    end
    invoke(serialize, Tuple{ClusterSerializer, Any}, s, f)
end

function serialize(s::ClusterSerializer, rr::RemoteChannel)
    p = worker_id_from_socket(s.io)
    (p !== rr.where) && send_add_client(rr, p)
    invoke(serialize, Tuple{ClusterSerializer, Any}, s, rr)
end

function deserialize(s::ClusterSerializer, t::Type{<:Future})
    fc = invoke(deserialize, Tuple{ClusterSerializer, DataType}, s, t) # deserialized copy
    f2 = Future(fc.where, RRID(fc.whence, fc.id), fc.v) # ctor adds to client_refs table

    # 1) send_add_client() is not executed when the ref is being serialized
    #    to where it exists, hence do it here.
    # 2) If we have received a 'fetch'ed Future or if the Future ctor found an
    #    already 'fetch'ed instance in client_refs (Issue #25847), we should not
    #    track it in the backing RemoteValue store.
    f2v_cache = @atomic f2.v
    if f2.where == myid() && f2v_cache === nothing
        add_client(remoteref_id(f2), myid())
    end
    f2
end

function deserialize(s::ClusterSerializer, t::Type{<:RemoteChannel})
    rr = invoke(deserialize, Tuple{ClusterSerializer, DataType}, s, t)
    if rr.where == myid()
        # send_add_client() is not executed when the ref is being
        # serialized to where it exists
        add_client(remoteref_id(rr), myid())
    end
    # call ctor to make sure this rr gets added to the client_refs table
    RemoteChannel{channel_type(rr)}(rr.where, RRID(rr.whence, rr.id))
end

# Future and RemoteChannel are serializable only in a running cluster.
# Serialize zeroed-out values to non ClusterSerializer objects
function serialize(s::AbstractSerializer, ::Future)
    zero_fut = Future((0,0,0,nothing))
    invoke(serialize, Tuple{AbstractSerializer, Any}, s, zero_fut)
end

function serialize(s::AbstractSerializer, ::RemoteChannel)
    zero_rc = RemoteChannel{Channel{Any}}((0,0,0))
    invoke(serialize, Tuple{AbstractSerializer, Any}, s, zero_rc)
end

# make a thunk to call f on args in a way that simulates what would happen if
# the function were sent elsewhere
function local_remotecall_thunk(f, args, kwargs)
    return ()->invokelatest(f, args...; kwargs...)
end

function remotecall(f, w::LocalProcess, args...; kwargs...)
    rr = Future(w)
    schedule_call(remoteref_id(rr), local_remotecall_thunk(f, args, kwargs))
    return rr
end

function remotecall(f, w::Worker, args...; kwargs...)
    rr = Future(w)
    send_msg(w, MsgHeader(remoteref_id(rr)), CallMsg{:call}(f, args, kwargs))
    return rr
end

"""
    remotecall(f, id::Integer, args...; kwargs...) -> Future

Call a function `f` asynchronously on the given arguments on the specified process.
Return a [`Future`](@ref).
Keyword arguments, if any, are passed through to `f`.
"""
remotecall(f, id::Integer, args...; kwargs...) = remotecall(f, worker_from_id(id), args...; kwargs...)

function remotecall_fetch(f, w::LocalProcess, args...; kwargs...)
    v=run_work_thunk(local_remotecall_thunk(f,args, kwargs), false)
    return isa(v, RemoteException) ? throw(v) : v
end

function remotecall_fetch(f, w::Worker, args...; kwargs...)
    # can be weak, because the program will have no way to refer to the Ref
    # itself, it only gets the result.
    oid = RRID()
    rv = lookup_ref(oid)
    rv.waitingfor = w.id
    send_msg(w, MsgHeader(RRID(0,0), oid), CallMsg{:call_fetch}(f, args, kwargs))
    v = take!(rv)
    lock(client_refs) do
        delete!(PGRP.refs, oid)
    end
    return isa(v, RemoteException) ? throw(v) : v
end

"""
    remotecall_fetch(f, id::Integer, args...; kwargs...)

Perform `fetch(remotecall(...))` in one message.
Keyword arguments, if any, are passed through to `f`.
Any remote exceptions are captured in a
[`RemoteException`](@ref) and thrown.

See also [`fetch`](@ref) and [`remotecall`](@ref).

# Examples
```julia-repl
\$ julia -p 2

julia> remotecall_fetch(sqrt, 2, 4)
2.0

julia> remotecall_fetch(sqrt, 2, -4)
ERROR: On worker 2:
DomainError with -4.0:
sqrt will only return a complex result if called with a complex argument. Try sqrt(Complex(x)).
...
```
"""
remotecall_fetch(f, id::Integer, args...; kwargs...) =
    remotecall_fetch(f, worker_from_id(id), args...; kwargs...)

remotecall_wait(f, w::LocalProcess, args...; kwargs...) = wait(remotecall(f, w, args...; kwargs...))

function remotecall_wait(f, w::Worker, args...; kwargs...)
    prid = RRID()
    rv = lookup_ref(prid)
    rv.waitingfor = w.id
    rr = Future(w)
    send_msg(w, MsgHeader(remoteref_id(rr), prid), CallWaitMsg(f, args, kwargs))
    v = fetch(rv.c)
    lock(client_refs) do
        delete!(PGRP.refs, prid)
    end
    isa(v, RemoteException) && throw(v)
    return rr
end

"""
    remotecall_wait(f, id::Integer, args...; kwargs...)

Perform a faster `wait(remotecall(...))` in one message on the `Worker` specified by worker id `id`.
Keyword arguments, if any, are passed through to `f`.

See also [`wait`](@ref) and [`remotecall`](@ref).
"""
remotecall_wait(f, id::Integer, args...; kwargs...) =
    remotecall_wait(f, worker_from_id(id), args...; kwargs...)

function remote_do(f, w::LocalProcess, args...; kwargs...)
    # the LocalProcess version just performs in local memory what a worker
    # does when it gets a :do message.
    # same for other messages on LocalProcess.
    thk = local_remotecall_thunk(f, args, kwargs)
    schedule(Task(thk))
    nothing
end

function remote_do(f, w::Worker, args...; kwargs...)
    send_msg(w, MsgHeader(), RemoteDoMsg(f, args, kwargs))
    nothing
end

"""
    remote_do(f, id::Integer, args...; kwargs...) -> nothing

Executes `f` on worker `id` asynchronously.
Unlike [`remotecall`](@ref), it does not store the
result of computation, nor is there a way to wait for its completion.

A successful invocation indicates that the request has been accepted for execution on
the remote node.

While consecutive `remotecall`s to the same worker are serialized in the order they are
invoked, the order of executions on the remote worker is undetermined. For example,
`remote_do(f1, 2); remotecall(f2, 2); remote_do(f3, 2)` will serialize the call
to `f1`, followed by `f2` and `f3` in that order. However, it is not guaranteed that `f1`
is executed before `f3` on worker 2.

Any exceptions thrown by `f` are printed to [`stderr`](@ref) on the remote worker.

Keyword arguments, if any, are passed through to `f`.
"""
remote_do(f, id::Integer, args...; kwargs...) = remote_do(f, worker_from_id(id), args...; kwargs...)

# have the owner of rr call f on it
function call_on_owner(f, rr::AbstractRemoteRef, args...)
    rid = remoteref_id(rr)
    if rr.where == myid()
        f(rid, args...)
    else
        remotecall_fetch(f, rr.where, rid, args...)
    end
end

function wait_ref(rid, caller, args...)
    v = fetch_ref(rid, args...)
    if isa(v, RemoteException)
        if myid() == caller
            throw(v)
        else
            return v
        end
    end
    nothing
end

"""
    wait(r::Future)

Wait for a value to become available for the specified [`Future`](@ref).
"""
wait(r::Future) = (v_cache = @atomic r.v; v_cache !== nothing && return r; call_on_owner(wait_ref, r, myid()); r)

"""
    wait(r::RemoteChannel, args...)

Wait for a value to become available on the specified [`RemoteChannel`](@ref).
"""
wait(r::RemoteChannel, args...) = (call_on_owner(wait_ref, r, myid(), args...); r)

"""
    fetch(x::Future)

Wait for and get the value of a [`Future`](@ref). The fetched value is cached locally.
Further calls to `fetch` on the same reference return the cached value. If the remote value
is an exception, throws a [`RemoteException`](@ref) which captures the remote exception and backtrace.
"""
function fetch(r::Future)
    v_cache = @atomic r.v
    v_cache !== nothing && return something(v_cache)

    if r.where == myid()
        rv, v_cache = @lock r.lock begin
            v_cache = @atomic :monotonic r.v
            rv = v_cache === nothing ? lookup_ref(remoteref_id(r)) : nothing
            rv, v_cache
        end

        if v_cache !== nothing
            return something(v_cache)
        else
            v_local = fetch(rv.c)
        end
    else
        v_local = call_on_owner(fetch_ref, r)
    end

    v_cache = @atomic r.v

    if v_cache === nothing # call_on_owner case
        v_old, status = @lock r.lock begin
            @atomicreplace r.v nothing => Some(v_local)
        end
        # status == true - when value obtained through call_on_owner
        # status == false - any other situation: atomicreplace fails, because by the time the lock is obtained cache will be populated
        # why? local put! performs caching and putting into channel under r.lock

        # for local put! use the cached value, for call_on_owner cases just take the v_local as it was just cached in r.v

        # remote calls getting the value from `call_on_owner` used to return the value directly without wrapping it in `Some(x)`
        # so we're doing the same thing here
        if status
            send_del_client(r)
            return v_local
        else # this `v_cache` is returned at the end of the function
            v_cache = v_old
        end
    end

    send_del_client(r)
    something(v_cache)
end

fetch_ref(rid, args...) = fetch(lookup_ref(rid).c, args...)

"""
    fetch(c::RemoteChannel)

Wait for and get a value from a [`RemoteChannel`](@ref). Exceptions raised are the
same as for a [`Future`](@ref). Does not remove the item fetched.
"""
fetch(r::RemoteChannel, args...) = call_on_owner(fetch_ref, r, args...)::eltype(r)

isready(rv::RemoteValue, args...) = isready(rv.c, args...)

"""
    put!(rr::Future, v)

Store a value to a [`Future`](@ref) `rr`.
`Future`s are write-once remote references.
A `put!` on an already set `Future` throws an `Exception`.
All asynchronous remote calls return `Future`s and set the
value to the return value of the call upon completion.
"""
function put!(r::Future, v)
    if r.where == myid()
        rid = remoteref_id(r)
        rv = lookup_ref(rid)
        isready(rv) && error("Future can be set only once")
        @lock r.lock begin
            put!(rv, v) # this notifies the tasks waiting on the channel in fetch
            set_future_cache(r, v) # set the cache before leaving the lock, so that the notified tasks already see it cached
        end
        del_client(rid, myid())
    else
        @lock r.lock begin # same idea as above if there were any local tasks fetching on this Future
            call_on_owner(put_future, r, v, myid())
            set_future_cache(r, v)
        end
    end
    r
end

function set_future_cache(r::Future, v)
    _, ok = @atomicreplace r.v nothing => Some(v)
    ok || error("internal consistency error detected for Future")
end

function put_future(rid, v, caller)
    rv = lookup_ref(rid)
    isready(rv) && error("Future can be set only once")
    put!(rv, v)
    # The caller has the value and hence can be removed from the remote store.
    del_client(rid, caller)
    nothing
end

put!(rv::RemoteValue, args...) = put!(rv.c, args...)
function put_ref(rid, caller, args...)
    rv = lookup_ref(rid)
    put!(rv, args...)
    if myid() == caller && rv.synctake !== nothing
        # Wait till a "taken" value is serialized out - github issue #29932
        lock(rv.synctake)
        unlock(rv.synctake)
    end
    nothing
end

"""
    put!(rr::RemoteChannel, args...)

Store a set of values to the [`RemoteChannel`](@ref).
If the channel is full, blocks until space is available.
Return the first argument.
"""
put!(rr::RemoteChannel, args...) = (call_on_owner(put_ref, rr, myid(), args...); rr)

# take! is not supported on Future

take!(rv::RemoteValue, args...) = take!(rv.c, args...)
function take_ref(rid, caller, args...)
    rv = lookup_ref(rid)
    synctake = false
    if myid() != caller && rv.synctake !== nothing
        # special handling for local put! / remote take! on unbuffered channel
        # github issue #29932
        synctake = true
        lock(rv.synctake)
    end

    v = try
        take!(rv, args...)
    catch e
        # avoid unmatched unlock when exception occurs
        # github issue #33972
        synctake && unlock(rv.synctake)
        rethrow(e)
    end

    isa(v, RemoteException) && (myid() == caller) && throw(v)

    if synctake
        return SyncTake(v, rv)
    else
        return v
    end
end

"""
    take!(rr::RemoteChannel, args...)

Fetch value(s) from a [`RemoteChannel`](@ref) `rr`,
removing the value(s) in the process.
"""
take!(rr::RemoteChannel, args...) = call_on_owner(take_ref, rr, myid(), args...)::eltype(rr)

# close and isopen are not supported on Future

close_ref(rid) = (close(lookup_ref(rid).c); nothing)
close(rr::RemoteChannel) = call_on_owner(close_ref, rr)

isopen_ref(rid) = isopen(lookup_ref(rid).c)
isopen(rr::RemoteChannel) = call_on_owner(isopen_ref, rr)

getindex(r::RemoteChannel) = fetch(r)
getindex(r::Future) = fetch(r)

getindex(r::Future, args...) = getindex(fetch(r), args...)
function getindex(r::RemoteChannel, args...)
    if r.where == myid()
        return getindex(fetch(r), args...)
    end
    return remotecall_fetch(getindex, r.where, r, args...)
end

max_stars_repo_name: BSnelling/julia

max_issues_repo_name: BSnelling/julia

max_forks_repo_name: BSnelling/julia


text:
"""
Thermdynamic data collection

Function names list (order by appearance in code):
    metacyc
    metacycIllParsed
    collect_formation_dg (internal eQuilibrator setup)
    collect_formation_dg (overload without internal eQuilibrator setup)
    collect_reaction_dg
    reaction_dg_bounds
    collect_dGr_bounds
    reaction_bounds
"""

"""
    $(SIGNATURES)
Returns reaction string with a prefixed "metacyc.compound:".
"""
function metacyc(str)
    eQuilibrator._parse_reaction_string(str, "metacyc.compound")
end

"""
    $(SIGNATURES)
Returns modified string with prefixed "metacyc.compound:". If the MetaCyc-based model was created via the Python-based "Moped" module, metabolite and reaction names may have artifacts ("__45__" and "__45__" substrings) that should be replaced by a hyphen or periodt, respectively, before ΔG estimations.
"""
function metacycIllParsed(str)
    nstr = replace(str, "__45__46__45__" => "-46-")
    nstr = replace(nstr, "__46__" => ".")
    nstr = replace(nstr, "__45__" => "-")
    nstr = replace(nstr, "__43__" => "+")
    eQuilibrator._parse_reaction_string(nstr, "metacyc.compound")
end

"""
    collect_formation_dg(
        database,
        db_id;
        dGf_dict = OrderedDict(),
        ph::Float64 = 7.0,
        pmg::Float64 = 2.0,
        i_strengh::Quantity = 100.0u"mM",
        temp::Quantity = 25u"°C",
    )

Collects ΔG of formation for all metabolites in `database` via Component Contribution. The data bank source must be passed by `db_id` (Options: `bigg`, `kegg`, `metanex`, `chebi`). Alternatively a dictionary of ΔGs of formation `dGf_dict` can be passed to extend it by new entries for the database metabolites.
Keyword arguments setup eQuilibrator conditions.

# Example
```
collect_formation_dg(
    ecoli_core,
    bigg
)
```
"""
function collect_formation_dg(
    database,
    db_id;
    dGf_dict = OrderedDict(),
    met_prefix = "",
    ph::Float64 = 7.0,
    pmg::Float64 = 2.0,
    i_strengh::Quantity = 100.0u"mM",
    temp::Quantity = 25u"°C",
)
    println("eQuilibrator is initialized. Please wait...")

    equil = eQuilibrator.Equilibrator(
        pH = ph,
        pMg = pmg,
        ionic_strength = i_strengh,
        temperature = temp,
    )

    n = length(metabolites(database))
    p = Progress(n, 1)
    for met in metabolites(database)
        if met ∉ keys(dGf_dict)
            rxn_str = string(" = ", met[length(met_prefix)+1:end-2])
            try
                @suppress begin
                    dGf_dict[met] =
                        standard_dg_prime(equil, db_id(rxn_str); balance_warn = false)
                end
            catch
                dGf_dict[met] = "NA"
            end
        end
        next!(p)
    end
    return dGf_dict
end

"""
    collect_formation_dg(
        database,
        equil,
        db_id;
        met_prefix = "",
        dGf_dict = OrderedDict(),
    )

Collects ΔG of formation for all metabolites in `database` via Component Contribution. The data bank source must be passed by `db_id` (Options: `bigg`, `kegg`, `metanex`, `chebi`). An eQuilibrator setup has to be passed by `equil` (see also: https://github.com/stelmo/Equilibrator.jl). Alternatively a dictionary of ΔGs of formation `dGf_dict` can be passed to extend it by new entries for the database metabolites.

# Example
```
collect_formation_dg(
    ecoli_core,
    eQuilibrator.Equilibrator(ionic_strength=150.0u"mM"),
    bigg
)
```
"""
function collect_formation_dg(
    database,
    equil,
    db_id;
    met_prefix = "",
    dGf_dict = OrderedDict(),
)
    println("ΔG of formation are collected.")
    n = length(metabolites(database))
    p = Progress(n, 1)
    for met in metabolites(database)
        if met ∉ keys(dGf_dict)
            rxn_str = string(" = ", met[length(met_prefix)+1:end-2])
            try
                @suppress begin
                    dGf_dict[met] =
                        standard_dg_prime(equil, db_id(rxn_str); balance_warn = false)
                end
            catch
                dGf_dict[met] = measurement(1000000u"kJ/mol", 10000000u"kJ/mol") #setting infeasible value
            end
        end
        next!(p)
    end
    return dGf_dict
end

"""
function collect_reaction_dg(
    database,
    equil,
    db_id;
    met_prefix = "",
    dgr_dict = OrderedDict(),
)
Collects standard transformed ΔG of reactions for every reaction in an SBML model. The eQuilibrator setup (see `eQuilibrator.Equilibrator()`) is passed as the `equil` argument. The source database (e.g. BiGG) is passed by `db_id`. If a dictionary with ΔGs of reaction already exist and should be expand, pass the dictionary to the `dgr_dict` keyword argument.
"""
function collect_reaction_dg(
    database,
    equil,
    db_id;
    met_prefix = "",
    dgr_dict = OrderedDict(),
)
    println("Standard ΔG of reactions are collected.")
    n = length(metabolites(database))
    p = Progress(n, 1)
    for rxn in reactions(database)
        if rxn ∉ keys(dgr_dict)
            rxn_string = stoichiometry_string(
                reaction_stoichiometry(database, rxn),
                format_id = x -> db_id(x[length(met_prefix)+1:end-2]),
            )
            try
                @suppress begin
                    dgr_dict[rxn] = standard_dg_prime(equil, rxn_string, balance_warn = false)
                end
            catch
                dgr_dict[rxn] = measurement(1000000u"kJ/mol", 10000000u"kJ/mol") #setting infeasible value
            end
        end
        next!(p)
    end
    return dgr_dict
end

"""
    collect_dGr_bounds(
        database;
        return_opts::Int = 3,
        ph::Float64 = 7.0,
        pmg::Float64 = 2.0,
        i_strength::Quantity = 100.0u"mM",
        temp::Quantity = 25u"°C",
        db_id = bigg,
        met_suffix = "",
    )

DEV NOTE: Necessary for MinFlux/MinRxn procedure, currently work in progress. Function might change severely.
Collects ΔG of reaction for every reaction in `database` via Component Contribution.
Returns dictionary of ΔG of reactions (`return_opts = 1`), the calculated bounds for MinFlux/MinRxn (`return_opts = 2`) or both (`return_opts = 3`).
Kwargs setup eQuilibrator conditions and return.
"""
function collect_dGr_bounds(
    database;
    return_opts::Int = 3,
    ph::Float64 = 7.0,
    pmg::Float64 = 2.0,
    i_strength::Quantity = 100.0u"mM",
    temp::Quantity = 25u"°C",
    db_id = bigg,
    met_prefix = "",
)
    println("eQuilibrator is initialized. Please wait...")
    equil = eQuilibrator.Equilibrator(
        pH = ph,
        pMg = pmg,
        ionic_strength = i_strength,
        temperature = temp,
    )

    R = 8.31446261815324u"J/(K*mol)"
    R = uconvert(u"kJ/(K*mol)", R)

    println("Calculating standard ΔG of reaction:")
    dGr_standard_dict = OrderedDict()
    n = length(reactions(database))
    p = Progress(n, 1)
    for rxn in reactions(database)
        rxn_string = stoichiometry_string(
            reaction_stoichiometry(database, rxn),
            format_id = x -> x[length(met_prefix)+1:end-2],
        )
        try
            @suppress begin
                dGr_standard_dict[rxn] =
                    standard_dg_prime(equil, db_id(rxn_string), balance_warn = false)
            end
        catch
            dGr_standard_dict[rxn] = "NA"
        end
        next!(p)
    end

    println("Calculating flux bounds:")
    dGr_bounds = OrderedDict()
    n = length(reactions(database))
    p2 = Progress(n, 1)

    for (rxn, dG) in dGr_standard_dict
        try
            @suppress begin
                reactants =
                    [abs(v) for (k, v) in reaction_stoichiometry(database, rxn) if v < 0]
                products =
                    [abs(v) for (k, v) in reaction_stoichiometry(database, rxn) if v > 0]
                qmax =
                    prod(fill(0.1, length(products)) .^ products) /
                    prod(fill(1 * 10^-6, length(reactants)) .^ reactants)
                qmin =
                    prod(fill(1 * 10^-6, length(products)) .^ products) /
                    prod(fill(0.1, length(reactants)) .^ reactants)
                u_calc = R * uconvert(u"K", temp) * log(qmax)
                l_calc = R * uconvert(u"K", temp) * log(qmin)
                println(typeof(dG))
                dGu = dG.val + u_calc.val
                dGl = dG.val + l_calc.val
                dGr_bounds[rxn] = (dGl, dGu)
            end
        catch
            dGr_bounds[rxn] = "NA"
        end
        next!(p2)
    end
    if return_opts == 3 # Returns both ΔG of reaction dictionary and the calculated bounds
        return dGr_standard_dict, dGr_bounds
    elseif return_opts == 2
        return dGr_bounds
    elseif return_opts == 1
        return dGr_standard_dict
    else
        println("Return options not specified. Switched to complete return.")
        return dGr_standard_dict, dGr_bounds
    end
end

function reduce_model(database, dg_dictionary)
    for rxn in reactions(database)
        if rxn ∉ keys(dg_dictionary)
            database = remove_reaction(database, rxn)
        end
    end
    return database
end

max_stars_repo_name: VictorVeraFrazao/OptStoic

max_issues_repo_name: VictorVeraFrazao/OptStoic

max_forks_repo_name: VictorVeraFrazao/OptStoic


text:
using JuMP, EAGO

                     m = Model(EAGO.Optimizer)

                     EAGO.register_eago_operators!(m)
                     set_optimizer_attribute(m, "verbosity", 4)
                     set_optimizer_attribute(m, "output_iterations", 1)

                     @variable(m, -1 <= x[i=1:3] <= 1)

                     @variable(m, -1.9515287248442508 <= q <= 2.66866690087751)

                     add_NL_constraint(m, :(swish(0.7910255593621245 + -0.9427161680575309*swish(-0.9881665794300214 + 0.6330293110094045*$(x[1]) + 0.8786704770929474*$(x[2]) + 0.11417841288574326*$(x[3])) + 0.08805479609820388*swish(0.6512390960566696 + -0.0007920608466527312*$(x[1]) + -0.28084831326063275*$(x[2]) + 0.17707194589795927*$(x[3]))) + swish(-0.7742678604770994 + 0.3023856535046492*swish(-0.9881665794300214 + 0.6330293110094045*$(x[1]) + 0.8786704770929474*$(x[2]) + 0.11417841288574326*$(x[3])) + -0.5348059615310992*swish(0.6512390960566696 + -0.0007920608466527312*$(x[1]) + -0.28084831326063275*$(x[2]) + 0.17707194589795927*$(x[3]))) - $q <= 0.0))

                     @objective(m, Min, q)
                     optimize!(m)
                     @show solve_time(m)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
#=
MBP Julia v1.6.1
[ Info: Benchmarks
Compile Native:
  0.000830 seconds (338 allocations: 22.781 KiB, 0.05% compilation time)
Compile BasicBlockInterpreter
  0.000006 seconds (7 allocations: 960 bytes)
Compilation Overhead: Native / BasicBlockInterpreter: 96.9x
F Call Native 1st
  0.040432 seconds (46.67 k allocations: 2.885 MiB, 25.46% gc time, 99.96% compilation time)
F Call Native 2nd
  0.000006 seconds (4 allocations: 176 bytes)
F Call Interpreter 1st
  0.054123 seconds (121.41 k allocations: 7.389 MiB, 99.51% compilation time)
F Call Interpreter 2nd
  0.000062 seconds (37 allocations: 976 bytes)
F Call Julia 1st
  0.010146 seconds (16.20 k allocations: 966.099 KiB, 99.81% compilation time)
F Call Julia 2nd
  0.000005 seconds (4 allocations: 176 bytes)
F Call Overhead: BasicBlockInterpreter / julia = 28.9x
F Call Overhead: Native / julia = 0.6x
[ Info: Compilation + F Call
BasicBlockInterpreter: 171.7µs
Native: 714.2µs
Native / BasicBlockInterpreter = 4.2x
=#

function benchmark_julia(x, y)
    v = zeros(Float64, 3)
    v[1] = x
    v[2] = y

    result = (((-( v[1] - v[2])) + 1.0 ) * 2.0) / 1.0

    return v, result
end

bb = OIR.BasicBlock()
in1 = OIR.ImmutableVariable(:x)
in2 = OIR.ImmutableVariable(:y)
OIR.addinput!(bb, in1)
OIR.addinput!(bb, in2)

argvec = OIR.addinstruction!(bb, OIR.call(OP_ZEROS, OIR.constant(Float64), OIR.constant(3)))
var_vec = OIR.MutableVariable(:v)
OIR.assign!(bb, var_vec, argvec)

OIR.addinstruction!(bb, OIR.call(OP_SETINDEX, var_vec, in1, OIR.constant(1)))
OIR.addinstruction!(bb, OIR.call(OP_SETINDEX, var_vec, in2, OIR.constant(2)))

arg1 = OIR.addinstruction!(bb, OIR.call(OP_GETINDEX, var_vec, OIR.constant(1)))
arg2 = OIR.addinstruction!(bb, OIR.call(OP_GETINDEX, var_vec, OIR.constant(2)))

# (((-( x[1] - x[2])) + 1.0 ) * 2.0) / 1.0
arg3 = OIR.addinstruction!(bb, OIR.call(OP_SUB, arg1, arg2))
arg4 = OIR.addinstruction!(bb, OIR.call(OP_SUB, arg3))
arg5 = OIR.addinstruction!(bb, OIR.call(OP_SUM, arg4, OIR.constant(1.0)))
arg6 = OIR.addinstruction!(bb, OIR.call(OP_MUL, arg5, OIR.constant(2.0)))
arg7 = OIR.addinstruction!(bb, OIR.call(OP_DIV, arg6, OIR.constant(1.0)))

var_result = OIR.ImmutableVariable(:result)
OIR.assign!(bb, var_result, arg7)

OIR.addoutput!(bb, var_vec)
OIR.addoutput!(bb, var_result)

println()
@info("Benchmarks")

println("Compile Native:")
@time fnative = OIR.compile(OIR.Native, bb)

println("Compile BasicBlockInterpreter")
@time finterpreter = OIR.compile(OIR.BasicBlockInterpreter, bb)

println("Compile BasicBlockInterpreter{Float64}")
@time finterpreter_float64 = OIR.compile(OIR.BasicBlockInterpreter{Union{Nothing, Float64, Vector{Float64}}}, bb)

compilation_el_native = @elapsed OIR.compile(OIR.Native, bb)
compilation_el_interpreter = @elapsed OIR.compile(OIR.BasicBlockInterpreter, bb)
println("Compilation Overhead: Native / BasicBlockInterpreter: $(round(compilation_el_native / compilation_el_interpreter, digits=1))x")

input = (10.0, 20.0)

println("F Call Native 1st")
@time result_native = fnative(input...)

println("F Call Native 2nd")
@time result_native = fnative(input...)

println("F Call Interpreter 1st")
@time result_interpreter = finterpreter(input...)

println("F Call Interpreter 2nd")
@time result_interpreter = finterpreter(input...)

println("F Call Interpreter Float64 1st")
@time result_interpreter = finterpreter_float64(input...)

println("F Call Interpreter Float64 2nd")
@time result_interpreter = finterpreter_float64(input...)

println("F Call Julia 1st")
@time result_julia = benchmark_julia(input...)

println("F Call Julia 2nd")
@time result_julia = benchmark_julia(input...)

@test result_native == result_julia
@test result_interpreter == result_julia

execution_el_interpreter = @elapsed finterpreter(input...)
execution_el_native = @elapsed fnative(input...)
execution_el_julia = @elapsed benchmark_julia(input...)
println("F Call Overhead: BasicBlockInterpreter / julia = $(round(execution_el_interpreter / execution_el_julia, digits=1))x")
println("F Call Overhead: Native / julia = $(round(execution_el_native / execution_el_julia, digits=1))x")

sec2µs(secs) = "$(round(secs*1E6, digits=1))µs"

@info("Compilation + F Call")
total_el_interpreter = compilation_el_interpreter + execution_el_interpreter
total_native = compilation_el_native + execution_el_native
println("BasicBlockInterpreter: $(sec2µs(total_el_interpreter))")
println("Native: $(sec2µs(total_native))")
println("Native / BasicBlockInterpreter = $(round(total_native / total_el_interpreter, digits=1))x")

println()


text:
# load packages
using CSV
using Random

include(pwd()*"/src/SDEMEM OU neuron data/load_neuron_data.jl")

"""
    set_up(M::Int=312, seed::Int=100)

Returns data and prior distributions.
"""
function set_up(M::Int=312, seed::Int=100)

    # define priors
    μ_0_1 = log(1/10); M_0_1 = 1; α_1 = 2; β_1 = 1
    μ_0_2 = log(1.5); M_0_2 = 1; α_2 = 2; β_2 = 1
    μ_0_3 = log(0.5); M_0_3 = 1; α_3 = 2; β_3 = 1

    prior_parameters_η = [μ_0_1 M_0_1 α_1 β_1;μ_0_2 M_0_2 α_2 β_2;μ_0_3 M_0_3 α_3 β_3]
    prior_parameters_σ_ϵ = [-1; 1/1]

    # load data
    ISI_all, nbr_ISE = load_neuron_data()
    y = ISI_all[1:M]

    # return: data and parameteres for priors
    return y,prior_parameters_η,prior_parameters_σ_ϵ

end

max_stars_repo_path: src/SDEMEM OU neuron data/ou_sdemem.jl

max_issues_repo_path: src/SDEMEM OU neuron data/ou_sdemem.jl

max_forks_repo_path: src/SDEMEM OU neuron data/ou_sdemem.jl


text:
# Nguyen, Lam M., et al. "SARAH: A novel method for machine learning problems using stochastic recursive gradient."
# International Conference on Machine Learning. PMLR, 2017.

using LinearAlgebra
using ProximalOperators
using ProximalAlgorithms.IterationTools
using Printf
using Base.Iterators
using Random
using Flux

export solution

include("SARAH_basic.jl") #? what are Array and AbstractArray?
include("SARAH_prox.jl")
include("SARAH_prox_DNN.jl")

struct SARAH{R<:Real}
    γ::Maybe{R}
    maxit::Int
    verbose::Bool
    freq::Int
    m::Maybe{Int}
    ꞵ::Maybe{R}
    DNN::Bool
    function SARAH{R}(;
        γ::Maybe{R} = nothing,
        maxit::Int = 10000,
        verbose::Bool = false,
        freq::Int = 1000,
        m::Maybe{Int} = nothing, # number of inner loop updates
        ꞵ::Maybe{R} = nothing,
        DNN::Bool = false,
    ) where {R}
        @assert γ === nothing || γ > 0
        @assert maxit > 0
        @assert freq > 0
        new(γ, maxit, verbose, freq, m, ꞵ, DNN)
    end
end

function (solver::SARAH{R})(
    x0::AbstractArray{C};
    F = nothing,
    g = ProximalOperators.Zero(),
    L = nothing,
    μ = nothing,
    N = nothing,
) where {R,C<:RealOrComplex{R}}

    stop(state::SARAH_basic_state) = false
    disp(it, state) = @printf "%5d | %.3e  \n" it state.γ

    F === nothing && (F = fill(ProximalOperators.Zero(), (N,)))
    m = solver.m === nothing ? m = N : m = solver.m

    maxit = solver.maxit

    # dispatching the structure
    if g == ProximalOperators.Zero()
        iter = SARAH_basic_iterable(F, x0, N, L, μ, solver.γ, m)
    else
        iter = SARAH_prox_iterable(F, g, x0, N, L, μ, solver.γ, m, solver.ꞵ)
    end

    iter = take(halt(iter, stop), maxit)
    iter = enumerate(iter)
    num_iters, state_final = nothing, nothing
    for (it_, state_) in iter  # unrolling the iterator
        # see https://docs.julialang.org/en/v1/manual/interfaces/index.html
        if solver.verbose && mod(it_, solver.freq) == 0
            disp(it_, state_)
        end
        num_iters, state_final = it_, state_
    end
    if solver.verbose && mod(num_iters, solver.freq) !== 0
        disp(num_iters, state_final)
    end # for the final iteration
    return solution(state_final), num_iters
end

"""
    SARAH([γ, maxit, verbose, freq, m])

Instantiate the SARAH algorithm  for solving (strongly) convex and nonconvex optimization problems of the form

    minimize 1/N sum_{i=1}^N f_i(x)

If `solver = SARAH(args...)`, then the above problem is solved with

	solver(x0, [F, g, N, L, μ])

where F is an array containing f_i's, x0 is the initial point, and L, μ are arrays of
smoothness and strong convexity moduli of f_i's; they are optional when γ is provided.

Optional keyword arguments are:
* `γ`: stepsize
* `L`: an array of smoothness moduli of f_i's. Not considered if γ is provided.
* `μ`: an array of strong convexity moduli of f_i's. Not considered if γ is provided.
    - if provided, the problem is considered as strongly convex,
    - if it is zero, convex problem is considred,
    - if nothing is provided, a nonconvex problem is considred.
* `maxit::Integer` (default: `10000`), maximum number of iterations to perform.
* `verbose::Bool` (default: `true`), whether or not to print information during the iterations.
* `freq::Integer` (default: `100`), frequency of verbosity.

"""

SARAH(::Type{R}; kwargs...) where {R} = SARAH{R}(; kwargs...)
SARAH(; kwargs...) = SARAH(Float64; kwargs...)

"""
If `solver = SARAH(args...)`, then
    itr = iterator(solver, x0, [F, g, N, L, μ])
is an iterable object. Note that [maxit, verbose, freq] fields of the solver are ignored here.

The solution at any given state can be obtained using solution(state), e.g.,
for state in Iterators.take(itr, maxit)
    # do something using solution(state)
end

See https://docs.julialang.org/en/v1/manual/interfaces/index.html
and https://docs.julialang.org/en/v1/base/iterators/ for a list of iteration utilities
"""

function iterator(
    solver::SARAH{R},
    x0::Union{AbstractArray{C},Tp};
    F = nothing,
    g = ProximalOperators.Zero(),
    L = nothing,
    μ = nothing,
    N = nothing,
    data = nothing,
    DNN_config = Tdnn
) where {R,C<:RealOrComplex{R},Tdnn,Tp}
    F === nothing && (F = fill(ProximalOperators.Zero(), (N,)))
    m = solver.m === nothing ? m = N : m = solver.m

    # dispatching the structure
    if g == ProximalOperators.Zero()
        println("basic version")
        iter = SARAH_basic_iterable(F, x0, N, L, μ, solver.γ, m)
    else
        if solver.DNN
            println("DNN version")
            iter = SARAH_prox_DNN_iterable(F, g, x0, N, L, μ, solver.γ, m, solver.ꞵ, data, DNN_config)
        else
            println("prox version")
            iter = SARAH_prox_iterable(F, g, x0, N, L, μ, solver.γ, m, solver.ꞵ)
        end
    end

    return iter
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "x265"
version = v"3.0"

# Collection of sources required to build x265Builder
sources = [
    ArchiveSource("http://ftp.videolan.org/pub/videolan/x265/x265_3.0.tar.gz",
                  "c5b9fc260cabbc4a81561a448f4ce9cad7218272b4011feabc3a6b751b2f0662"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/x265*/

# We need `nasm` for x86_64
apk add nasm

mkdir bld && cd bld
cmake -DCMAKE_INSTALL_PREFIX="${prefix}" -DCMAKE_TOOLCHAIN_FILE="${CMAKE_TARGET_TOOLCHAIN}" -DENABLE_PIC=ON -DENABLE_SHARED=ON ../source
make -j${nproc}
make install
# Remove the large static archive
rm  ${prefix}/lib/libx265.a
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    ExecutableProduct("x265", :x265),
    LibraryProduct("libx265", :libx265)
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: waralex/Yggdrasil

max_issues_repo_name: waralex/Yggdrasil

max_forks_repo_name: waralex/Yggdrasil


text:
# @macroexpand @recipe(Coastlines) do scene
#     merge(
#         default_theme(scene, Lines),
#         Theme(crs = (source = LonLat(), dest = LonLat()),)
#     )
# end

"""
    coastlines(; kwargs...)

Plot coastlines along the world.
"""
coastlines(; kwargs...) = coastlines!(Scene(); kwargs...)

const Coastlines{ArgType} = AbstractPlotting.Combined{coastlines, ArgType}

(Base).show(var"#34#io"::IO, ::Type{<:Coastlines}) = begin
        #= /Users/anshul/.julia/dev/AbstractPlotting/src/recipes.jl:154 =#
        AbstractPlotting.print(var"#34#io", "Coastlines", "{...}")
end

function coastlines(var"#39#attributes"...)
    #= /Users/anshul/.julia/dev/AbstractPlotting/src/recipes.jl:13 =#
    var"#36#attr" = AbstractPlotting.Attributes(var"#39#attributes")
    #= /Users/anshul/.julia/dev/AbstractPlotting/src/recipes.jl:14 =#
    var"#37#kw" = AbstractPlotting.extract_scene_attributes!(var"#36#attr")
    #= /Users/anshul/.julia/dev/AbstractPlotting/src/recipes.jl:15 =#
    AbstractPlotting.plot!(AbstractPlotting.Scene(; var"#37#kw"...), Coastlines, var"#36#attr", 1)
end

#= /Users/anshul/.julia/dev/AbstractPlotting/src/recipes.jl:19 =#
function coastlines!(; var"#41#attributes"...)
    #= /Users/anshul/.julia/dev/AbstractPlotting/src/recipes.jl:20 =#
    AbstractPlotting.plot!(AbstractPlotting.current_scene(), Coastlines, AbstractPlotting.Attributes(var"#41#attributes"), 1)
end
#= /Users/anshul/.julia/dev/AbstractPlotting/src/recipes.jl:23 =#
function coastlines!(var"#42#scene"::AbstractPlotting.SceneLike; var"#44#attributes"...)
    #= /Users/anshul/.julia/dev/AbstractPlotting/src/recipes.jl:24 =#
    AbstractPlotting.plot!(var"#42#scene", Coastlines, AbstractPlotting.Attributes(var"#44#attributes"), 1)
end
#= /Users/anshul/.julia/dev/AbstractPlotting/src/recipes.jl:27 =#
function coastlines(var"#47#attributes"::AbstractPlotting.Attributes, var"#48#args"...; var"#49#kw_attributes"...)
    #= /Users/anshul/.julia/dev/AbstractPlotting/src/recipes.jl:28 =#
    var"#45#merged" = AbstractPlotting.merge!(AbstractPlotting.Attributes(var"#49#kw_attributes"), var"#47#attributes")
    #= /Users/anshul/.julia/dev/AbstractPlotting/src/recipes.jl:29 =#
    var"#46#kw" = AbstractPlotting.extract_scene_attributes!(var"#45#merged")
    #= /Users/anshul/.julia/dev/AbstractPlotting/src/recipes.jl:30 =#
    AbstractPlotting.plot!(AbstractPlotting.Scene(; var"#46#kw"...), Coastlines, var"#45#merged", 1)
end
#= /Users/anshul/.julia/dev/AbstractPlotting/src/recipes.jl:33 =#
function coastlines!(var"#50#attributes"::AbstractPlotting.Attributes, var"#51#args"...; var"#52#kw_attributes"...)
    #= /Users/anshul/.julia/dev/AbstractPlotting/src/recipes.jl:34 =#
    AbstractPlotting.plot!(AbstractPlotting.current_scene(), Coastlines, AbstractPlotting.merge!(AbstractPlotting.Attributes(var"#52#kw_attributes"), var"#50#attributes"), 1)
end
#= /Users/anshul/.julia/dev/AbstractPlotting/src/recipes.jl:37 =#
function coastlines!(var"#53999#scene"::AbstractPlotting.SceneLike, var"#54999#attributes"::AbstractPlotting.Attributes, var"#55999#args"...; var"#56999#kw_attributes"...)
    #= /Users/anshul/.julia/dev/AbstractPlotting/src/recipes.jl:38 =#
    AbstractPlotting.plot!(var"#53999#scene", Coastlines, AbstractPlotting.merge!(AbstractPlotting.Attributes(var"#56999#kw_attributes"), var"#54999#attributes"), 1)
end

AbstractPlotting.default_theme(var"#57999#scene", ::Type{<:Coastlines}) = begin
        (((scene,)->begin
                merge(default_theme(scene, Lines), Theme(crs = (source = LonLat(), dest = LonLat())))
            end))(var"#57999#scene")
    end
export Coastlines, coastlines, coastlines!

function AbstractPlotting.plot!(p::Coastlines)

    @extract p (crs,)

    new_linevec = Observable(COASTLINES_LINEVEC)

    lift(crs.source, crs.dest) do source, dest
        new_linevec[] = Proj4.transform.(source, dest, COASTLINES_LINEVEC)
    end

    lines!(p, new_linevec)
end

"""
    earth(; kwargs...)

Plots an image of the Earth in the given projection.
"""
earth(; kwargs...) = earth!(Scene(); kwargs...)

const Earth{ArgType} = AbstractPlotting.Combined{earth, ArgType}

(Base).show(var"#34#io"::IO, ::Type{<:Earth}) = begin
        #= /Users/anshul/.julia/dev/AbstractPlotting/src/recipes.jl:154 =#
        AbstractPlotting.print(var"#34#io", "Earth", "{...}")
end

function earth(var"#39#attributes"...)
    #= /Users/anshul/.julia/dev/AbstractPlotting/src/recipes.jl:13 =#
    var"#36#attr" = AbstractPlotting.Attributes(var"#39#attributes")
    #= /Users/anshul/.julia/dev/AbstractPlotting/src/recipes.jl:14 =#
    var"#37#kw" = AbstractPlotting.extract_scene_attributes!(var"#36#attr")
    #= /Users/anshul/.julia/dev/AbstractPlotting/src/recipes.jl:15 =#
    AbstractPlotting.plot!(AbstractPlotting.Scene(; var"#37#kw"...), Earth, var"#36#attr", 1)
end

#= /Users/anshul/.julia/dev/AbstractPlotting/src/recipes.jl:19 =#
function earth!(; var"#41#attributes"...)
    #= /Users/anshul/.julia/dev/AbstractPlotting/src/recipes.jl:20 =#
    AbstractPlotting.plot!(AbstractPlotting.current_scene(), Earth, AbstractPlotting.Attributes(var"#41#attributes"), 1)
end
#= /Users/anshul/.julia/dev/AbstractPlotting/src/recipes.jl:23 =#
function earth!(var"#42#scene"::AbstractPlotting.SceneLike; var"#44#attributes"...)
    #= /Users/anshul/.julia/dev/AbstractPlotting/src/recipes.jl:24 =#
    AbstractPlotting.plot!(var"#42#scene", Earth, AbstractPlotting.Attributes(var"#44#attributes"), 1)
end
#= /Users/anshul/.julia/dev/AbstractPlotting/src/recipes.jl:27 =#
function earth(var"#47#attributes"::AbstractPlotting.Attributes, var"#48#args"...; var"#49#kw_attributes"...)
    #= /Users/anshul/.julia/dev/AbstractPlotting/src/recipes.jl:28 =#
    var"#45#merged" = AbstractPlotting.merge!(AbstractPlotting.Attributes(var"#49#kw_attributes"), var"#47#attributes")
    #= /Users/anshul/.julia/dev/AbstractPlotting/src/recipes.jl:29 =#
    var"#46#kw" = AbstractPlotting.extract_scene_attributes!(var"#45#merged")
    #= /Users/anshul/.julia/dev/AbstractPlotting/src/recipes.jl:30 =#
    AbstractPlotting.plot!(AbstractPlotting.Scene(; var"#46#kw"...), Earth, var"#45#merged", 1)
end
#= /Users/anshul/.julia/dev/AbstractPlotting/src/recipes.jl:33 =#
function earth!(var"#50#attributes"::AbstractPlotting.Attributes, var"#51#args"...; var"#52#kw_attributes"...)
    #= /Users/anshul/.julia/dev/AbstractPlotting/src/recipes.jl:34 =#
    AbstractPlotting.plot!(AbstractPlotting.current_scene(), Earth, AbstractPlotting.merge!(AbstractPlotting.Attributes(var"#52#kw_attributes"), var"#50#attributes"), 1)
end
#= /Users/anshul/.julia/dev/AbstractPlotting/src/recipes.jl:37 =#
function earth!(var"#53999#scene"::AbstractPlotting.SceneLike, var"#54999#attributes"::AbstractPlotting.Attributes, var"#55999#args"...; var"#56999#kw_attributes"...)
    #= /Users/anshul/.julia/dev/AbstractPlotting/src/recipes.jl:38 =#
    AbstractPlotting.plot!(var"#53999#scene", Earth, AbstractPlotting.merge!(AbstractPlotting.Attributes(var"#56999#kw_attributes"), var"#54999#attributes"), 1)
end

AbstractPlotting.default_theme(var"#57999#scene", ::Type{<:Earth}) = begin
        (((scene,)->begin
        merge(
            default_theme(scene, Surface),
            Theme(
                crs = (
                    source = LonLat(),
                    dest = LonLat()
                ),
                shading = false
            )
        )
            end))(var"#57999#scene")
    end
export Earth, earth, earth!

function AbstractPlotting.plot!(p::Earth)

    @extract p (crs,)

    lons = LinRange(-180, 179, size(EARTH_IMG)[2])
    lats = LinRange(89.5, -89.5, size(EARTH_IMG)[1])

    xs = Observable([lon for lat in lats, lon in lons])
    ys = Observable([lat for lat in lats, lon in lons])

    lift(crs.source, crs.dest) do source, dest

        xs.val = [lon for lat in lats, lon in lons]
        ys.val = [lat for lat in lats, lon in lons]
        Proj4.transform!(source, dest, vec(xs.val), vec(ys.val))

        xs[] = xs[]
        ys[] = ys[]
    end

    surface!(p, xs, ys, zeros(size(xs[])); color = EARTH_IMG, shading = false, show_axis = false)
end


text:
@testset "Base Excel configuration tests" begin
    # File doesn't exist.
    mpSim, isOkay = initialiseFromExcel( joinpath( "base", "baseConfig0" ) )
    @test !isOkay

    # Excel file has no General tab.
    mpSim, isOkay = initialiseFromExcel( joinpath( "base", "baseConfig1" ) )
    @test !isOkay
end  # @testset "Base Excel configuration tests"

@testset "Database parameter test" begin
    mpSim = MPsim()
    filePath = joinpath( "base", "baseConfig" )
    fileName = string( filePath, ".xlsx" )

    XLSX.openxlsx( fileName ) do xf
        @test MPIO.readDBpars( mpSim, xf[ "General" ], filePath ) == "Excel"
    end  # XLSX.openXLSX( fileName ) do xf

    @test mpSim.dbName == joinpath( "base", "baseConfig", "simDB.sqlite" )
    @test mpSim.simName == "sim"
end  # @testset "Database parameter test"

@testset "General parameter test" begin
    mpSim = MPsim()
    filePath = joinpath( "base", "baseConfig" )
    fileName = string( filePath, ".xlsx" )

    XLSX.openxlsx( fileName ) do xf
        MPIO.readGeneralPars( mpSim, xf[ "General" ], filePath )
    end  # XLSX.openXLSX( fileName ) do xf

    @test mpSim.catFileName == joinpath( "base", "catalogue.xlsx" )
    @test mpSim.personnelTarget == 0
    @test mpSim.simLength == 300.0
end  # @testset "General parameter test"


text:
mutable struct Document
	"Document mutable struct"

	"terms:   A vector{Int} containing keys for the Corpus vocab Dict."
	"counts:  A Vector{Int} denoting the counts of each term in the Document."
	"readers: A Vector{Int} denoting the keys for the Corpus users Dict."
	"ratings: A Vector{Int} denoting the ratings for each reader in the Document."
	"title:   The title of the document (String)."

	terms::Vector{Int}
	counts::Vector{Int}
	readers::Vector{Int}
	ratings::Vector{Int}
	title::String

	function Document(;terms=[], counts=ones(length(terms)), readers=[], ratings=ones(length(readers)), title="")
		doc = new(terms, counts, readers, ratings, title)
		check_doc(doc)
		return doc
	end
end

### Document outer constructors.
Document(terms) = Document(terms=terms)

struct DocumentError <: Exception
    msg::String
end

Base.showerror(io::IO, e::DocumentError) = print(io, "DocumentError: ", e.msg)

function check_doc(doc::Document)
	"Check Document parameters."

	all(doc.terms .> 0)									|| throw(DocumentError("All terms must be positive integers."))
	all(doc.counts .> 0)								|| throw(DocumentError("All counts must be positive integers."))
	isequal(length(doc.terms), length(doc.counts))		|| throw(DocumentError("The terms and counts vectors must have the same length."))
	all(doc.readers .> 0)								|| throw(DocumentError("All readers must be positive integers."))
	all(doc.ratings .> 0)								|| throw(DocumentError("All ratings must be positive integers."))
	isequal(length(doc.readers), length(doc.ratings))	|| throw(DocumentError("The readers and ratings vectors must have the same length."))
 	nothing
 end

mutable struct Corpus
	"Corpus mutable struct."

	"docs:  A Vector{Document} containing the documents which belong to the Corpus."
	"vocab: A Dict{Int, String} containing a mapping term Int (key) => term String (value)."
	"users: A Dict{Int, String} containing a mapping user Int (key) => user String (value)."

	docs::Vector{Document}
	vocab::Dict{Int, String}
	users::Dict{Int, String}

	function Corpus(;docs=Document[], vocab=[], users=[])
		isa(vocab, Vector) && (vocab = Dict(vkey => term for (vkey, term) in enumerate(string.(vocab))))
		isa(users, Vector) && (users = Dict(ukey => user for (ukey, user) in enumerate(string.(users))))

		corp = new(docs, vocab, users)
		check_docs(corp)

		all(collect(keys(corp.vocab)) .> 0) || throw(CorpusError("All vocab keys must be positive integers."))
		all(collect(keys(corp.users)) .> 0) || throw(CorpusError("All user keys must be positive integers."))
		return corp
	end
end

### Corpus outer constructors.
Corpus(docs::Vector{Document}) = Corpus(docs=docs)
Corpus(docs::Vector{Document}; vocab=[], users=[]) = Corpus(docs=docs, vocab=vocab, users=users)
Corpus(doc::Document) = Corpus(docs=[doc])
Corpus(doc::Document; vocab=[], users=[]) = Corpus(docs=[doc], vocab=vocab, users=users)

struct CorpusError <: Exception
    msg::String
end

Base.showerror(io::IO, e::CorpusError) = print(io, "CorpusError: ", e.msg)

 function check_docs(corp::Corpus)
 	"Check Corpus documents."

	for (d, doc) in enumerate(corp)
		try
			check_doc(doc)
		catch
			throw(CorpusError("Document $d failed check."))
		end
	end
end

function check_corp(corp::Corpus)
	"Check Corpus parameters."

	check_docs(corp)

	all(collect(keys(corp.vocab)) .> 0) || throw(CorpusError("All vocab keys must be positive integers."))
	all(collect(keys(corp.users)) .> 0) || throw(CorpusError("All user keys must be positive integers."))

	issubset(vcat([doc.terms for doc in corp]...), keys(corp.vocab)) 		|| throw(CorpusError("Documents contain term keys not found in corpus vocabulary (see fixcorp! function)."))
	issubset(vcat([doc.readers for doc in corp]...), keys(corp.users)) 		|| throw(CorpusError("Documents contain user keys not found in corpus users (see fixcorp! function)."))
	(length(corp.vocab) == maximum(push!(collect(keys(corp.vocab)), 0)))	|| throw(CorpusError("Corpus vocab keys must form unit range starting at 1 (see fixcorp! function)."))
	(length(corp.users) == maximum(push!(collect(keys(corp.users)), 0)))	|| throw(CorpusError("Corpus user keys must form unit range starting at 1 (see fixcorp! function)."))
	nothing
end

Base.show(io::IO, doc::Document) = print(io, "Document with:\n * $(length(doc.terms)) terms\n * $(length(doc.readers)) readers")
Base.length(doc::Document) = length(doc.terms)
Base.size(doc::Document) = sum(doc.counts)
Base.in(doc::Document, corp::Corpus) = in(doc, corp.docs)
Base.isempty(doc::Document) = length(doc) == 0

Base.show(io::IO, corp::Corpus) = print(io, "Corpus with:\n * $(length(corp)) docs\n * $(length(corp.vocab)) vocab\n * $(length(corp.users)) users")
Base.iterate(corp::Corpus, d=1) = Base.iterate(corp.docs, d)
Base.push!(corp::Corpus, doc::Document) = push!(corp.docs, doc)
Base.pop!(corp::Corpus) = pop!(corp.docs)
Base.pushfirst!(corp::Corpus, doc::Document) = pushfirst!(corp.docs, doc)
Base.pushfirst!(corp::Corpus, docs::Vector{Document}) = pushfirst!(corp.docs, docs)
Base.popfirst!(corp::Corpus) = popfirst!(corp.docs)
Base.insert!(corp::Corpus, d::Int, doc::Document) = insert!(corp.docs, d, doc)
Base.deleteat!(corp::Corpus, d::Int) = deleteat!(corp.docs, d)
Base.deleteat!(corp::Corpus, doc_indices::Vector{Int}) = deleteat!(corp.docs, doc_indices)
Base.deleteat!(corp::Corpus, doc_indices::UnitRange{Int}) = deleteat!(corp.docs, doc_indices)
Base.getindex(corp::Corpus, d::Int) = getindex(corp.docs, d)
Base.getindex(corp::Corpus, doc_indices::Vector{Int}) = getindex(corp.docs, doc_indices)
Base.getindex(corp::Corpus, doc_indices::UnitRange{Int}) = getindex(corp.docs, doc_indices)
Base.getindex(corp::Corpus, doc_bool::Vector{Bool}) = getindex(corp.docs, doc_bool)
Base.setindex!(corp::Corpus, doc::Document, d::Int) = setindex!(corp.docs, doc, d)
Base.setindex!(corp::Corpus, docs::Vector{Document}, doc_indices::Vector{Int}) = setindex!(corp.docs, docs, doc_indices)
Base.setindex!(corp::Corpus, docs::Vector{Document}, doc_indices::UnitRange{Int}) = setindex!(corp.docs, docs, doc_indices)
Base.findfirst(corp::Corpus, doc::Document) = findfirst(corp.docs, doc)
Base.findall(corp::Corpus, docs::Vector{Document}) = findall((in)corp.docs, docs)
Base.findall(corp::Corpus, doc::Document) = findall(corp, [doc])
Base.length(corp::Corpus) = length(corp.docs)
Base.size(corp::Corpus) = (length(corp), length(corp.vocab), length(corp.users))
Base.copy(corp::Corpus) = Corpus(docs=copy(corp.docs), vocab=copy(corp.vocab), users=copy(corp.users))
Base.lastindex(corp::Corpus) = length(corp)
Base.enumerate(corp::Corpus) = enumerate(corp.docs)
Base.unique(corp::Corpus) = Corpus(docs=unique(corp.docs), vocab=corp.vocab, users=corp.users)

function showdocs(corp::Corpus, docs::Vector{Document})
	"Display document(s) in readable format."

	issubset(vcat([doc.terms for doc in docs]...), keys(corp.vocab)) || throw(DocumentError("Some documents contain term keys not found in Corpus vocabulary."))

	for (n, doc) in enumerate(docs)
		if !isempty(doc.title)
			@juliadots "$(doc.title)\n"

		else
			@juliadots "Document\n"
		end

		if !isempty(doc)
			println(Crayon(bold=false), join([corp.vocab[vkey] for vkey in doc.terms], " "))

		else
			println()
		end

		if n < length(docs)
			println()
		end
	end
end

function showdocs(corp::Corpus, doc_indices::Vector{<:Integer})
	"Display document(s) in readable format."

	issubset(doc_indices, 1:length(corp))											|| throw(CorpusError("Some document indices outside corpus range."))
	issubset(vcat([doc.terms for doc in corp[doc_indices]]...), keys(corp.vocab))	|| throw(DocumentError("Some documents contain term keys not found in Corpus vocabulary."))

	for (n, d) in enumerate(doc_indices)
		@juliadots "Document $d\n"

		if !isempty(corp[d].title)
			@juliadots "$(corp[d].title)\n"
		end

		if !isempty(corp[d])
			println(Crayon(bold=false), join([corp.vocab[vkey] for vkey in corp[d].terms], " "))

		else
			println()
		end

		if n < length(doc_indices)
			println()
		end
	end
end

showdocs(corp::Corpus, doc::Document) = showdocs(corp, [doc])
showdocs(corp::Corpus, doc_range::UnitRange{<:Integer}) = showdocs(corp, collect(doc_range))
showdocs(corp::Corpus, d::Integer) = showdocs(corp, [d])
showdocs(corp::Corpus) = showdocs(corp, 1:length(corp))

function showtitles(corp::Corpus, docs::Vector{Document})
	"Display document title(s) in readable format."

	issubset(vcat([doc.terms for doc in docs]...), keys(corp.vocab)) || throw(DocumentError("Some documents contain term keys not found in Corpus vocabulary."))

	for doc in docs
		print(Crayon(foreground=:yellow, bold=true), " • ")

		if !isempty(doc.title)
			println(Crayon(foreground=:white, bold=false), "$(doc.title)")

		else
			println(Crayon(foreground=:white, bold=true), "Document")
		end
	end
end

function showtitles(corp::Corpus, doc_indices::Vector{<:Integer})
	"Display document title(s) in readable format."

	issubset(doc_indices, 1:length(corp))											|| throw(CorpusError("Some document indices outside corpus range."))
	issubset(vcat([doc.terms for doc in corp[doc_indices]]...), keys(corp.vocab))	|| throw(DocumentError("Some documents contain term keys not found in Corpus vocabulary."))

	for d in doc_indices
		print(Crayon(foreground=:yellow, bold=true), " • ")

		if !isempty(corp[d].title)
			print(Crayon(foreground=:white, bold=true), "Document $d ")
			println(Crayon(foreground=:white, bold=false), "$(corp[d].title)")

		else
			println(Crayon(foreground=:white, bold=true), "Document $d")
		end
	end
end

showtitles(corp::Corpus, doc::Document) = showtitles(corp, [doc])
showtitles(corp::Corpus, doc_range::UnitRange{<:Integer}) = showtitles(corp, collect(doc_range))
showtitles(corp::Corpus, d::Integer) = showtitles(corp, [d])
showtitles(corp::Corpus) = showtitles(corp, 1:length(corp))

getvocab(corp::Corpus) = sort(collect(values(corp.vocab)))
getusers(corp::Corpus) = sort(collect(values(corp.users)))

function readcorp(;docfile::String="", vocabfile::String="", userfile::String="", titlefile::String="", delim::Char=',', counts::Bool=false, readers::Bool=false, ratings::Bool=false)
	"Load a Corpus object from text file(s)."

	(ratings <= readers) || (ratings = false; warn("Ratings require readers, ratings switch set to false."))
	(!isempty(docfile) | isempty(titlefile)) || warn("No docfile, titles will not be assigned.")

	corp = Corpus()
	if !isempty(docfile)
		docs = open(docfile)
		doc_kwargs = [:terms, :counts, :readers, :ratings]

		for (d, doc_block) in enumerate(Iterators.partition(readlines(docs), counts + readers + ratings + 1))
			try
				doc_lines = Vector{Int}[[parse(Int, p) for p in split(line, delim)] for line in doc_block]
				doc_input = zip(doc_kwargs[[true, counts, readers, ratings]], doc_lines)
				push!(corp, Document(;doc_input...))

			catch
				throw(CorpusError("Document $d beginning on line $((d - 1) * (counts + readers + ratings) + d) failed to load."))
			end
		end

	else
		warn("No docfile, topic models cannot be trained without documents.")
	end

	if !isempty(vocabfile)
		vocab = readdlm(vocabfile, '\t', comments=false)
		vkeys = vocab[:,1]
		terms = [string(j) for j in vocab[:,2]]
		corp.vocab = Dict{Int, String}(zip(vkeys, terms))
		all(collect(keys(corp.vocab)) .> 0) || throw(CorpusError("All vocab keys must be positive integers."))
	end

	if !isempty(userfile)
		users = readdlm(userfile, '\t', comments=false)
		ukeys = users[:,1]
		users = [string(u) for u in users[:,2]]
		corp.users = Dict{Int, String}(zip(ukeys, users))
		all(collect(keys(corp.users)) .> 0) || throw(CorpusError("All user keys must be positive integers."))
	end

	if !isempty(titlefile)
		titles = readdlm(titlefile, '\n', String)
		for (d, doc) in enumerate(corp)
			doc.title = titles[d]
		end
	end

	return corp
end

function readcorp(corp_symbol::Symbol)
	"Shortcuts for prepackaged corpora."

	datasets_path = joinpath((@__DIR__)[1:end-3], "datasets")

	if corp_symbol == :nsf
		docfile = joinpath(datasets_path, "nsf/nsfdocs.txt")
		vocabfile = joinpath(datasets_path, "nsf/nsfvocab.txt")
		titlefile = joinpath(datasets_path, "nsf/nsftitles.txt")
		corp = readcorp(docfile=docfile, vocabfile=vocabfile, titlefile=titlefile, counts=true)

	elseif corp_symbol == :citeu
		docfile = joinpath(datasets_path, "citeu/citeudocs.txt")
		vocabfile = joinpath(datasets_path, "citeu/citeuvocab.txt")
		userfile = joinpath(datasets_path, "citeu/citeuusers.txt")
		titlefile = joinpath(datasets_path, "citeu/citeutitles.txt")
		corp = readcorp(docfile=docfile, vocabfile=vocabfile, userfile=userfile, titlefile=titlefile, counts=true, readers=true)

	else
		println("Included corpora:\n:nsf\n:citeu")
		corp = nothing
	end

	return corp
end

function writecorp(corp::Corpus; docfile::String="", vocabfile::String="", userfile::String="", titlefile::String="", delim::Char=',', counts::Bool=false, readers::Bool=false, ratings::Bool=false)
	"Write a corpus."

	(ratings <= readers) || (ratings = false; warn("Ratings require readers, ratings switch set to false."))

	if !isempty(docfile)
		doc_kwargs = [:counts, :readers, :ratings]
		doc_kwargs = doc_kwargs[[counts, readers, ratings]]
		docfile = open(docfile, "w")
		for doc in corp
			write(docfile, join(doc.terms, delim), '\n')
			for arg in doc_kwargs
				write(docfile, join(doc.(arg), delim), '\n')
			end
		end
		close(docfile)
	end

	if !isempty(vocabfile)
		vkeys = sort(collect(keys(corp.vocab)))
		vocab = zip(vkeys, [corp.vocab[vkey] for vkey in vkeys])
		writedlm(vocabfile, vocab)
	end

	if !isempty(userfile)
		ukeys = sort(collect(keys(corp.users)))
		users = zip(ukeys, [corp.users[ukey] for ukey in ukeys])
		writedlm(userfile, users)
	end

	if !isempty(titlefile)
		writedlm(titlefile, [doc.title for doc in corp])
	end
	nothing
end

### The _corp and _docs functions are designed to provide safe methods for modifying corpora.
### The _corp functions only meaningfully modify the Corpus object.
### In so far as the _corp functions modify Document objects, it only amounts to a possible relabeling of the keys in the documents.
### The _doc functions only modify the Document objects attached the corpus, they do not modify the Corpus object.
### The exception to the above rule is the stop_corp! function, which removes stop words from both the Corpus vocab dictionary and associated keys in the documents.

function abridge_corp!(corp::Corpus, n::Integer=0)
	"All terms which appear less than n times in the corpus are removed from all documents."

	doc_vkeys = Set(vcat([doc.terms for doc in unique(corp)]...))
	vocab_count = Dict(Int(j) => 0 for j in doc_vkeys)

	for doc in unique(corp), (j, c) in zip(doc.terms, doc.counts)
		vocab_count[j] += c
	end

	for doc in unique(corp)
		keep = Bool[vocab_count[j] >= n for j in doc.terms]
		doc.terms = doc.terms[keep]
		doc.counts = doc.counts[keep]
	end
	nothing
end

function alphabetize_corp!(corp::Corpus; vocab::Bool=true, users::Bool=true)
	"Alphabetize vocab and/or user dictionaries."

	if vocab
		vkeys = sort(collect(keys(corp.vocab)))
		terms = sort(collect(values(corp.vocab)))

		vkey_map = Dict(vkey_old => vkey_new for (vkey_old, vkey_new) in zip(vkeys, vkeys[sortperm(sortperm([corp.vocab[vkey] for vkey in vkeys]))]))
		corp.vocab = Dict(vkey => term for (vkey, term) in zip(vkeys, terms))

		for doc in unique(corp)
			doc.terms = [vkey_map[j] for j in doc.terms]
		end
	end

	if users
		ukeys = sort(collect(keys(corp.users)))
		users = sort(collect(values(corp.users)))

		ukey_map = Dict(ukey_old => ukey_new for (ukey_old, ukey_new) in zip(ukeys, ukeys[sortperm(sortperm([corp.users[ukey] for ukey in ukeys]))]))
		corp.users = Dict(ukey => user for (ukey, user) in zip(ukeys, users))

		for doc in unique(corp)
			doc.readers = [ukey_map[r] for r in doc.readers]
		end
	end
	nothing
end

function clip_corp!(corp::Corpus; vocab::Bool=true, users::Bool=true)
	"Keep top n terms."
	nothing
end

function remove_terms!(corp::Corpus; terms::Vector{String}=[])
	"Vocab keys for specified terms are removed from all documents."

	remove_keys = filter(vkey -> lowercase(corp.vocab[vkey]) in terms, collect(keys(corp.vocab)))

	for doc in unique(corp)
		keep = Bool[!(j in remove_keys) for j in doc.terms]
		doc.terms = doc.terms[keep]
		doc.counts = doc.counts[keep]
	end
	nothing
end

remove_terms!(corp::Corpus, term::String) = remove_vocab!(corp, terms=[term])

function compact_corp!(corp::Corpus; vocab::Bool=true, users::Bool=true)
	"Relabel vocab and/or user keys so that they form a unit range starting at 1."

	if vocab
		vkeys = sort(collect(keys(corp.vocab)))
		vkey_map = Dict(vkey => v for (v, vkey) in enumerate(vkeys))
		corp.vocab = Dict(vkey_map[vkey] => corp.vocab[vkey] for vkey in keys(corp.vocab))

		for doc in unique(corp)
			doc.terms = [vkey_map[vkey] for vkey in doc.terms]
		end
	end

	if users
		ukeys = sort(collect(keys(corp.users)))
		ukey_map = Dict(ukey => u for (u, ukey) in enumerate(ukeys))
		corp.users = Dict(ukey_map[ukey] => corp.users[ukey] for ukey in keys(corp.users))

		for doc in unique(corp)
			doc.readers = [ukey_map[ukey] for ukey in doc.readers]
		end
	end
	nothing
end

function condense_corp!(corp::Corpus)
	"Ignore term order in documents."
	"Multiple seperate occurrences of terms are stacked and their associated counts increased."

	for doc in unique(corp)
		docdict = Dict(j => 0 for j in doc.terms)
		for (j, c) in zip(doc.terms, doc.counts)
			docdict[j] += c
		end

		doc.terms = collect(keys(docdict))
		doc.counts = collect(values(docdict))
	end
	nothing
end

function pad_corp!(corp::Corpus; vocab::Bool=true, users::Bool=true)
	"Enter generic values for vocab and/or user keys which appear in documents but not in the vocab/user dictionaries."

	if vocab
		doc_vkeys = Set(vcat([doc.terms for doc in corp]...))
		for vkey in setdiff(doc_vkeys, keys(corp.vocab))
			corp.vocab[vkey] = string(join(["#term", vkey]))
		end
	end

	if users
		doc_ukeys = Set(vcat([doc.readers for doc in corp]...))
		for ukey in setdiff(doc_ukeys, keys(corp.users))
			corp.users[ukey] = string(join(["#user", ukey]))
		end
	end
	nothing
end

function remove_empty_docs!(corp::Corpus)
	"Documents with no terms are removed from the corpus."

	keep = Bool[length(doc) > 0 for doc in corp]
	corp.docs = corp[keep]
	nothing
end

function remove_redundant!(corp::Corpus; vocab::Bool=true, users::Bool=true)
	"Remove vocab and/or user keys which map to redundant values."
	"Reassign Document term and/or reader keys."

	if vocab
		vkey_map = Dict()
		vocab_dict_inverse = Dict()
		vkeys = sort(collect(keys(corp.vocab)))
		for vkey in vkeys
			if corp.vocab[vkey] in keys(vocab_dict_inverse)
				vkey_map[vkey] = vocab_dict_inverse[corp.vocab[vkey]]
				delete!(corp.vocab, vkey)

			else
				vkey_map[vkey] = vkey
				vocab_dict_inverse[corp.vocab[vkey]] = vkey
			end
		end

		for doc in unique(corp)
			doc.terms = [vkey_map[vkey] for vkey in doc.terms]
		end
	end

	if users
		ukey_map = Dict()
		users_dict_inverse = Dict()
		ukeys = sort(collect(keys(corp.users)))
		for ukey in ukeys
			if corp.users[ukey] in keys(users_dict_inverse)
				ukey_map[ukey] = users_dict_inverse[corp.users[ukey]]
				delete!(corp.users, ukey)

			else
				ukey_map[ukey] = ukey
				users_dict_inverse[corp.users[ukey]] = ukey
			end
		end

		for doc in unique(corp)
			doc.readers = [ukey_map[ukey] for ukey in doc.readers]
		end
	end
	nothing
end

function stop_corp!(corp::Corpus)
	"Filter stop words in the associated corpus."

	datasets_path = joinpath((@__DIR__)[1:end-3], "datasets")

	stop_words = vec(readdlm(joinpath(datasets_path, "stopwords.txt"), String))
	stop_keys = filter(vkey -> lowercase(corp.vocab[vkey]) in stop_words, collect(keys(corp.vocab)))

	for doc in unique(corp)
		keep = Bool[!(j in stop_keys) for j in doc.terms]
		doc.terms = doc.terms[keep]
		doc.counts = doc.counts[keep]
	end
	nothing
end

function trim_corp!(corp::Corpus; vocab::Bool=true, users::Bool=true)
	"Those keys which appear in the corpus vocab and/or user dictionaries but not in any of the documents are removed from the corpus."

	if vocab
		doc_vkeys = Set(vcat([doc.terms for doc in corp]...))
		corp.vocab = Dict(vkey => corp.vocab[vkey] for vkey in intersect(keys(corp.vocab), doc_vkeys))
	end

	if users
		doc_ukeys = Set(vcat([doc.readers for doc in corp]...))
		corp.users = Dict(ukey => corp.users[ukey] for ukey in intersect(keys(corp.users), doc_ukeys))
	end
	nothing
end

function trim_docs!(corp::Corpus; terms::Bool=true, readers::Bool=true)
	"Those vocab and/or user keys which appear in documents but not in the corpus dictionaries are removed from the documents."

	if terms
		doc_vkeys = Set(vcat([doc.terms for doc in corp]...))
		bogus_vkeys = setdiff(doc_vkeys, keys(corp.vocab))
		for doc in unique(corp)
			keep = Bool[!(j in bogus_vkeys) for j in doc.terms]
			doc.terms = doc.terms[keep]
			doc.counts = doc.counts[keep]
		end
	end

	if readers
		doc_ukeys = Set(vcat([doc.readers for doc in corp]...))
		bogus_ukeys = setdiff(doc_ukeys, keys(corp.users))
		for doc in unique(corp)
			keep = Bool[!(u in bogus_ukeys) for u in doc.readers]
			doc.readers = doc.readers[keep]
			doc.ratings = doc.ratings[keep]
		end
	end
	nothing
end

function fixcorp!(corp::Corpus; vocab::Bool=true, users::Bool=true, abridge::Integer=0, alphabetize::Bool=false, condense::Bool=false, pad::Bool=false, remove_empty_docs::Bool=false, remove_redundant::Bool=false, remove_terms::Vector{String}=String[], stop::Bool=false, trim::Bool=false)
	"Generic function to ensure that a Corpus object can be loaded into a TopicModel object."
	"Either pad_corp! or trim_docs!."
	"compact_corp!."
	"Contains other optional keyword arguments."

	check_docs(corp)

	all(collect(keys(corp.vocab)) .> 0) || throw(CorpusError("All vocab keys must be positive integers."))
	all(collect(keys(corp.users)) .> 0) || throw(CorpusError("All user keys must be positive integers."))

	pad ? pad_corp!(corp) : trim_docs!(corp)

	remove_redundant			&& remove_redundant!(corp)
	condense 					&& condense_corp!(corp)
	abridge > 0 				&& abridge_corp!(corp, abridge)
	pad 						&& pad_corp!(corp, vocab=vocab, users=users)
	length(remove_terms) > 0	&& remove_terms!(corp, terms=remove_terms)
	stop 						&& stop_corp!(corp)
	trim 						&& trim_corp!(corp, vocab=vocab, users=users)
	alphabetize 				&& alphabetize_corp!(corp, vocab=vocab, users=users)
	remove_empty_docs 			&& remove_empty_docs!(corp)

	compact_corp!(corp)
	nothing
end


text:
using MPI
using Lazy
using Logging

export QXMPIContext
export get_rank_size, get_rank_start, get_rank_range

# Implementation of QXMPIContext which provides a Context that can be used with MPI
#
# The QXMPIContext struct contains a reference to another context struct which is used
# to perform contractions for individual sets of slice values
#
# Utility functions get_rank_size, get_rank_start and get_rank_range are provided as
# generic utility functions for working with multiple processing ranks

"""MPI context which can be used to perform distribted computations"""
struct QXMPIContext <: AbstractContext
    serial_ctx::AbstractContext
    comm::MPI.Comm # communicator containing all ranks
    sub_comm::MPI.Comm # communicator containing ranks in subgroup, used for partitions
    root_comm::MPI.Comm # communicator containing ranks of matching rank from other sub_comms
end

"""
    QXMPIContext(ctx::QXContext, comm::MPI.Comm, sub_comm_size::Int=1)

Constructor for QXMPIContext that initialises new sub-communicators for managing groups
of nodes.
"""
function QXMPIContext(ctx::QXContext,
                      comm::Union{MPI.Comm, Nothing}=nothing;
                      sub_comm_size::Int=1)
    if comm === nothing
        if !MPI.Initialized()
            MPI.Init()
            @info "MPI Initialised"
        end
        comm = MPI.COMM_WORLD
    end
    @info "Number processes $(MPI.Comm_size(MPI.COMM_WORLD))"
    @assert MPI.Comm_size(comm) % sub_comm_size == 0 "sub_comm_size must divide comm size evenly"
    sub_comm = MPI.Comm_split(comm, MPI.Comm_rank(comm) ÷ sub_comm_size, MPI.Comm_rank(comm) % sub_comm_size)
    root_comm = MPI.Comm_split(comm, MPI.Comm_rank(comm) % sub_comm_size, MPI.Comm_rank(comm) ÷ sub_comm_size)
    QXMPIContext(ctx, comm, sub_comm, root_comm)
end

######################################################################################
#  Forward each of the methods to work with serial_ctx from the QXMPIContext
######################################################################################
@forward QXMPIContext.serial_ctx gettensor
@forward QXMPIContext.serial_ctx settensor!
@forward QXMPIContext.serial_ctx Base.getindex
@forward QXMPIContext.serial_ctx Base.setindex!
@forward QXMPIContext.serial_ctx Base.haskey
@forward QXMPIContext.serial_ctx Base.zeros
@forward QXMPIContext.serial_ctx Base.zero
@forward QXMPIContext.serial_ctx Base.eltype
@forward QXMPIContext.serial_ctx set_open_bonds!
@forward QXMPIContext.serial_ctx set_slice_vals!

"""Make struct callable"""
(ctx::QXMPIContext)(args...; kwargs...) = ctx.serial_ctx(args...; kwargs...)

"""
    compute_amplitude!(ctx, bitstring::String; num_slices=nothing)

Calculate a single amplitude with the given context and bitstring. Involves a sum over
contributions from each slice. Can optionally set the number of bonds. By default all slices
are used.
"""
function compute_amplitude!(ctx::QXMPIContext, bitstring::String; max_slices=nothing)
    set_open_bonds!(ctx, bitstring)
    amplitude = nothing
    si = SliceIterator(ctx.serial_ctx.cg, max_slices=max_slices)
    r = get_comm_range(ctx.sub_comm, length(si))
    for p in SliceIterator(si, r.start, r.stop)
        set_slice_vals!(ctx, p)
        if amplitude === nothing
            amplitude = ctx()
        else
            amplitude += ctx()
        end
    end
    # reduce across sub_communicator
    # TODO: replace reduce with batched reduce or use one-sided communication to accumulate
    # to root of sub_comm
    MPI.Reduce!(amplitude, +, 0, ctx.sub_comm)
    if MPI.Comm_rank(ctx.sub_comm) != 0
        return nothing
    end
    if ndims(amplitude) == 0 amplitude = amplitude[] end
    amplitude
end

"""
    get_rank_size(n::Integer, size::Integer, rank::Integer)

Partition n items among processes of communicator of size size and return the size of the
given rank. Algorithm used is to:
1. Divide items equally among processes
2. Spread remainder over ranks in ascending order of rank number
"""
function get_rank_size(n::Integer, size::Integer, rank::Integer)
    (n ÷ size) + ((n % size) >= (rank + 1))
end

"""
    get_rank_start(n::Integer, size::Integer, rank::Integer)

Partition n items among processes of communicator of size size and return the starting of the
given rank. Algorithm used is to:
1. Divide items equally among processes
2. Spread remainder over ranks in ascending order of rank number
"""
function get_rank_start(n::Integer, size::Integer, rank::Integer)
    start = rank * (n ÷ size) + 1
    start + min(rank, n % size)
end

"""
    get_rank_range(n::Integer, size::Integer, rank::Integer)

Partition n items among processes of communicator of size size and return a range over these
indices.
"""
function get_rank_range(n::Integer, size::Integer, rank::Integer)
    start = get_rank_start(n, size, rank)
    UnitRange(start, start + get_rank_size(n, size, rank) - 1)
end

"""
    get_comm_range(comm::MPI.Comm, n::Integer)

Given a communicator and number of items, get the range for local rank
"""
function get_comm_range(comm::MPI.Comm, n::Integer)
    get_rank_range(n, MPI.Comm_size(comm), MPI.Comm_rank(comm))
end

"""
    ctxmap(f, ctx::QXMPIContext, items)

For each of the items in the local range apply the f function and return
the result.
"""
function ctxmap(f, ctx::QXMPIContext, items)
    map(f, items[get_comm_range(ctx.root_comm, length(items))])
end

"""
    ctxgather(ctx::QXMPIContext, items)

Gather local items to root rank
"""
function ctxgather(ctx::QXMPIContext, items)
    if MPI.Comm_rank(ctx.sub_comm) == 0
        return MPI.Gather(items, 0, ctx.root_comm)
    end
end

"""
    ctxreduce(f, ctx::QXMPIContext, items)

Reduce across items with funciton f
"""
function ctxreduce(f, ctx::QXMPIContext, items)
    if MPI.Comm_rank(ctx.sub_comm) == 0
        return MPI.Reduce(items, f, 0, MPI.root_comm)
    end
end