text:
module Rewrite

include("mkbuild.jl")

function rename_package!(filepath, newname)
    lines = readlines(filepath)
    for (ind,line) in emumerate(lines)
        if line[1:7] == "package"
            lines[ind] = string("package"," ", newname)
        end
    end

    open(filepath,"w") do file
        write(file, lines)
    end

    Union{}
end

function remove_comment!(filepath)
    lines = readlines(filepath)
    st = ed = 0
    state = false
    for (ind, line) in emumerate(lines)
        if line[1:2] == "/*"
            st = ind
            state = true
        elseif st != 0 && line[1:2] == "* "
            @assert state == true
        elseif line[1:2] == "*/"
            @assert state == true
            ed = ind
        end
    end
    lines = vcat(lines[1:st-1],lines[ed+1:end])
    write(filepath, lines)
end

end


text:
module RenderJay

const SMALL_NUM = 0.000001
const GREAT_NUM = 100000.0
const p = 1.0/(2.0 * pi)
const INV_PI = 1.0/pi

using CSV, LightXML, DataFrames, DelimitedFiles, LinearAlgebra, ProgressMeter, SharedArrays, ImageMagick, Images, Distributions

include("structs.jl")
export normalize
export Point
export Line
export Ray
export Coord
export Asset
export Hook
export Transform
export Bvh
export Camera
export Geometry
export Disk
export Ball
export Cylinder
export Cone
export Mesh
export Walls

include("info.jl")
export get_extremes

include("transform.jl")
export weibull
export carthesian2spherical
export spherical2carthesian
export relazi
export translate
export transform
export transform_forward
export transform_inverse
export rotate_forward
export rotate_inverse

include("create.jl")
export create_ray
export create_rays
export create_coords

include("bvh.jl")
export make_boxes
export produce_bvh_data
export median_cut
export build_containing_box
export produce_bounding_boxes
export make_bvh

include("light.jl")
export make_sky
export background_radiance
export sample_skymap

include("trace.jl")
export trace_back
export trace_forward
export record

include("render.jl")
export render_pixel
export render_image

include("brdf.jl")
export cross_scalars
export normalize_scalars
export project2normal
export sample_f
export propagate
export compute_propagation_of_reflectance
export compute_propagation_of_transmittance
export distr
export geom
export fresnel
export cooktorrance
export compute_half_vector
export ct_sampler
export pure_reflection
export rpv

include("intersect.jl")
export hit_box
export intersect_bvh
export intersect
export intersect_scene

include("shader.jl")
export lambertian_shader
export bilambertian_shader
export pure_reflection_shader
export rpv_shader
export lambertian_full_fate_shader
export bilambertian_full_fate_shader
export rpv_full_fate_shader
export lightsource_shader

include("io.jl")
export read_shaders
export read_items
export read_items_no_load_bvh
export read_items_load_bvh
export read_camera
export read_skymap
export read_scene
export read_cylinder_data
export read_cylinder
export read_cone_data
export read_cone
export read_mesh_data
export read_mesh
export write_bvh
export read_bvh

end # module

ext: jl

max_stars_repo_name: martinvanleeuwen/RenderJay

max_issues_repo_name: martinvanleeuwen/RenderJay

max_forks_repo_name: martinvanleeuwen/RenderJay


text:
abstract type LinearSolver end

abstract type AbstractSCSType end

Base.cconvert(::Type{Ptr{Cvoid}}, x::AbstractSCSType) = x

function Base.unsafe_convert(::Type{Ptr{Cvoid}}, x::AbstractSCSType)
    return pointer_from_objref(x)
end

mutable struct ScsMatrix{T} <: AbstractSCSType
    values::Ptr{Cdouble}
    rowval::Ptr{T}
    colptr::Ptr{T}
    m::T
    n::T
end

mutable struct ScsSettings{T} <: AbstractSCSType
    normalize::T # boolean, heuristic data rescaling
    scale::Cdouble # if normalized, rescales by this factor
    rho_x::Cdouble # x equality constraint scaling
    max_iters::T # maximum iterations to take
    eps::Cdouble # convergence tolerance
    alpha::Cdouble # relaxation parameter
    cg_rate::Cdouble # for indirect, tolerance goes down like (1/iter)^cg_rate
    verbose::T # boolean, write out progress
    warm_start::T # boolean, warm start (put initial guess in Sol struct)
    acceleration_lookback::T # acceleration memory parameter
    write_data_filename::Cstring
    ScsSettings{T}() where {T} = new{T}()
end

mutable struct ScsData{T} <: AbstractSCSType
    m::T
    n::T
    A::Ptr{Cvoid}
    b::Ptr{Cdouble}
    c::Ptr{Cdouble}
    stgs::Ptr{Cvoid}
end

mutable struct ScsSolution <: AbstractSCSType
    x::Ptr{Cdouble}
    y::Ptr{Cdouble}
    s::Ptr{Cdouble}
end

mutable struct ScsInfo{T} <: AbstractSCSType
    iter::T
    status::NTuple{32,Cchar}
    statusVal::T
    pobj::Cdouble
    dobj::Cdouble
    resPri::Cdouble
    resDual::Cdouble
    resInfeas::Cdouble
    resUnbdd::Cdouble
    relGap::Cdouble
    setupTime::Cdouble
    solveTime::Cdouble
    ScsInfo{T}() where {T} = new{T}()
end

mutable struct ScsCone{T} <: AbstractSCSType
    f::T
    l::T
    q::Ptr{T}
    qsize::T
    s::Ptr{T}
    ssize::T
    ep::T
    ed::T
    p::Ptr{Cdouble}
    psize::T
end

mutable struct Solution{T}
    x::Array{Float64,1}
    y::Array{Float64,1}
    s::Array{Float64,1}
    info::ScsInfo{T}
    ret_val::T
end

"""
    _ScsDataWrapper

A type for wrapping all data inputs to SCS to prevent them from being freed by
the garbage collector during a solve.

You should not construct this manually. Call `scs_solve` instead.
"""
struct _ScsDataWrapper{S,T}
    linear_solver::S
    m::T
    n::T
    A::ScsMatrix{T}
    values::Vector{Cdouble}
    rowval::Vector{T}
    colptr::Vector{T}
    b::Vector{Cdouble}
    c::Vector{Cdouble}
    f::T
    l::T
    q::Vector{T}
    s::Vector{T}
    ep::T
    ed::T
    p::Vector{Cdouble}
    primal::Vector{Cdouble}
    dual::Vector{Cdouble}
    slack::Vector{Cdouble}
    settings::ScsSettings
    options::Any
end

function _sanitize_options(options)
    option_dict = Dict{Symbol,Any}()
    fields = fieldnames(ScsSettings)
    for (key, value) in options
        if key == :linear_solver
            continue
        elseif !(key in fields)
            throw(ArgumentError("Unrecognized option $(key)"))
        end
        option_dict[key] = value
    end
    return option_dict
end

function raw_status(info::ScsInfo)
    data = UInt8[info.status[i] for i in 1:findfirst(iszero, info.status)-1]
    return String(data)
end

function _to_sparse(::Type{T}, A::AbstractMatrix) where {T}
    return _to_sparse(T, SparseArrays.sparse(A))
end

function _to_sparse(::Type{T}, A::SparseArrays.SparseMatrixCSC) where {T}
    rowval = convert(Vector{T}, A.rowval .- 1)
    colptr = convert(Vector{T}, A.colptr .- 1)
    return A.nzval, rowval, colptr
end

"""
    scs_solve(args...)

SCS solves a problem of the form
```
minimize        c' * x
subject to      A * x + s = b
                s in K
```
where K is a product cone of
 * zero cone { x | x = 0 },
 * positive orthant { x | x >= 0 },
 * second-order cones (SOC) { (t,x) | ||x||_2 <= t },
 * semi-definite cones (SDC) { X | X psd }, and
 * exponential cones {(x,y,z) | y e^(x/y) <= z, y>0 },
 * power cone { (x,y,z) | x^a*y^(1-a) >= |z|, x>=0, y>=0 },
 * dual power cone { (u,v,w) | (u/a)^a * (v/(1-a))^(1-a) >= |w|, u>=0, v>=0}

Description of input argments:
 * A is the matrix with m rows and n cols
 * b is a vector of length m
 * c is a vector of length n

The rows of A correspond to cones in K and need to be specified in the order
above.

 * f (num primal zero / dual free cones, i.e. primal equality constraints)
 * l (num linear cones)
 * q (array of SOCs sizes)
 * s (array of SDCs sizes)
 * ep (num primal exponential cones)
 * ed (num dual exponential cones)
 * p (array of power cone params, must be in [-1, 1], negative values specify
   the dual cone)

Returns a Solution object.

!!! warning
    SCS expects the semi-definite cones to be scaled by a factor of √2. That is,
    the off-diagonal elements in the A matrix and b vector should be multiplied
    by √2, and then the corresponding rows in the dual and slack solution
    vectors should be multiplied by 1/√2.
"""
function scs_solve(
    linear_solver::Type{<:LinearSolver},
    m::Integer,
    n::Integer,
    A,
    b::Vector{Float64},
    c::Vector{Float64},
    f::Integer,
    l::Integer,
    q::Vector{<:Integer},
    s::Vector{<:Integer},
    ep::Integer,
    ed::Integer,
    p::Vector{Float64},
    primal_sol::Vector{Float64} = zeros(n),
    dual_sol::Vector{Float64} = zeros(m),
    slack::Vector{Float64} = zeros(m);
    warm_start::Bool = false,
    options...,
)
    if n <= 0
        throw(ArgumentError("The number of variables must be greater than 0"))
    elseif m <= 0
        throw(ArgumentError("The number of constraints must be greater than 0"))
    end
    if length(primal_sol) == n && length(dual_sol) == length(slack) == m
        if !warm_start
            fill!(primal_sol, 0.0)
            fill!(dual_sol, 0.0)
            fill!(slack, 0.0)
        end
    else
        if warm_start
            throw(ArgumentError("Warmstart doesn't match the problem size"))
        end
        primal_sol, dual_sol, slack = zeros(n), zeros(m), zeros(m)
    end
    T = scsint_t(linear_solver)
    m, n, f, l, ep, ed = T(m), T(n), T(f), T(l), T(ep), T(ed)
    values, rowval, colptr = _to_sparse(T, A)
    option_dict = _sanitize_options(options)
    if warm_start
        option_dict[:warm_start] = 1
    end
    model = _ScsDataWrapper(
        linear_solver,
        m,
        n,
        ScsMatrix(pointer(values), pointer(rowval), pointer(colptr), m, n),
        values,
        rowval,
        colptr,
        b,
        c,
        f,
        l,
        convert(Vector{T}, q),
        convert(Vector{T}, s),
        ep,
        ed,
        p,
        primal_sol,
        dual_sol,
        slack,
        ScsSettings{T}(),
        option_dict,
    )
    Base.GC.@preserve model begin
        return _unsafe_scs_solve(model)
    end
end

function _unsafe_scs_solve(model::_ScsDataWrapper{S,T}) where {S,T}
    scs_solution = ScsSolution(
        pointer(model.primal),
        pointer(model.dual),
        pointer(model.slack),
    )
    scs_cone = ScsCone{T}(
        model.f,
        model.l,
        pointer(model.q),
        length(model.q),
        pointer(model.s),
        length(model.s),
        model.ep,
        model.ed,
        pointer(model.p),
        length(model.p),
    )
    scs_info = ScsInfo{T}()
    scs_data = ScsData{T}(
        model.m,
        model.n,
        pointer_from_objref(model.A),
        pointer(model.b),
        pointer(model.c),
        pointer_from_objref(model.settings),
    )
    scs_set_default_settings(model.linear_solver, scs_data)
    for (key, value) in model.options
        setfield!(
            model.settings,
            key,
            convert(fieldtype(ScsSettings{T}, key), value),
        )
    end
    p = scs_init(model.linear_solver, scs_data, scs_cone, scs_info)
    status = scs_solve(
        model.linear_solver,
        p,
        scs_data,
        scs_cone,
        scs_solution,
        scs_info,
    )
    scs_finish(model.linear_solver, p)
    return Solution(model.primal, model.dual, model.slack, scs_info, status)
end


text:
@testset "1480.running-sum-of-1d-array.jl" begin
    @test running_sum([1,2,3,4]) == [1,3,6,10]
    @test running_sum([1,1,1,1,1]) == [1,2,3,4,5]
    @test running_sum([3,1,2,10,1]) == [3,4,6,16,17]
end


text:
# Utility #
# ========================== #
function eval_objective(objective::JuMP.GenericQuadExpr, x::AbstractVector)
    aff = objective.aff
    val = aff.constant
    for (i,var) in enumerate(aff.vars)
        val += aff.coeffs[i]*x[var.col]
    end
    return val
end

function fill_solution!(stochasticprogram::StochasticProgram)
    dep = DEP(stochasticprogram)
    # First stage
    first_stage = get_stage_one(stochasticprogram)
    nrows, ncols = length(first_stage.linconstr), first_stage.numCols
    first_stage.objVal = dep.objVal
    first_stage.colVal = dep.colVal[1:ncols]
    first_stage.redCosts = dep.redCosts[1:ncols]
    first_stage.linconstrDuals = dep.linconstrDuals[1:nrows]
    # Second stages
    fill_solution!(scenarioproblems(stochasticprogram), dep.colVal[ncols+1:end], dep.redCosts[ncols+1:end], dep.linconstrDuals[nrows+1:end])
    return nothing
end
function fill_solution!(scenarioproblems::ScenarioProblems, x::AbstractVector, μ::AbstractVector, λ::AbstractVector)
    cbegin = 0
    rbegin = 0
    for (i,subproblem) in enumerate(subproblems(scenarioproblems))
        snrows, sncols = length(subproblem.linconstr), subproblem.numCols
        subproblem.colVal = x[cbegin+1:cbegin+sncols]
        subproblem.redCosts = μ[cbegin+1:cbegin+sncols]
        subproblem.linconstrDuals = λ[rbegin+1:rbegin+snrows]
        subproblem.objVal = eval_objective(subproblem.obj, subproblem.colVal)
        cbegin += sncols
        rbegin += snrows
    end
    return nothing
end
function fill_solution!(scenarioproblems::DScenarioProblems, x::AbstractVector, μ::AbstractVector, λ::AbstractVector)
    cbegin = 0
    rbegin = 0
    @sync begin
        for w in workers()
            wncols = remotecall_fetch((sp)->sum([s.numCols::Int for s in fetch(sp).problems]), w, scenarioproblems[w-1])
            wnrows = remotecall_fetch((sp)->sum([length(s.linconstr)::Int for s in fetch(sp).problems]), w, scenarioproblems[w-1])
            crange = cbegin+1:cbegin+wncols
            rrange = rbegin+1:rbegin+wnrows
            @async remotecall_fetch((sp,x,μ,λ)->fill_solution!(fetch(sp),x,μ,λ),
                                    w,
                                    scenarioproblems[w-1],
                                    x[crange],
                                    μ[crange],
                                    λ[rrange])
            cbegin += wncols
            rbegin += wnrows
        end
    end
    return nothing
end

function calculate_objective_value!(stochasticprogram::StochasticProgram)
    first_stage = get_stage_one(stochasticprogram)
    objective_value = eval_objective(first_stage.obj, first_stage.colVal)
    objective_value += calculate_subobjectives(scenarioproblems(stochasticprogram))
    first_stage.objVal = objective_value
    return nothing
end
function calculate_subobjectives(scenarioproblems::ScenarioProblems)
    return sum([(probability(scenario)*eval_objective(subprob.obj,subprob.colVal))::Float64 for (scenario,subprob) in zip(scenarios(scenarioproblems),subproblems(scenarioproblems))])
end
function calculate_subobjectives(scenarioproblems::DScenarioProblems)
    partial_subobjectives = Vector{Float64}(undef, nworkers())
    @sync begin
        for (i,w) in enumerate(workers())
            @async partial_subobjectives[i] = remotecall_fetch((sp) -> calculate_subobjectives(fetch(sp)),
                                                               w,
                                                               scenarioproblems[w-1])
        end
    end
    return sum(partial_subobjectives)
end

function invalidate_cache!(stochasticprogram::StochasticProgram)
    cache = problemcache(stochasticprogram)
    delete!(cache, :evp)
    delete!(cache, :dep)
    return nothing
end

function remove_stage!(stochasticprogram::StochasticProgram{N}, s::Integer) where N
    1 <= s <= N || error("Stage $s not in range 1 to $N.")
    if s == 1
        haskey(stochasticprogram.problemcache, :stage_1) || return nothing
        delete!(stochasticprogram.problemcache, :stage_1)
    else
        haskey(stochasticprogram.problemcache, :stage_1) || return nothing
        remove_subproblems!(stochasticprogram, s)
    end
    s < N && clear_parent!(scenarioproblems(stochasticprogram, s+1))
    return nothing
end

function remove_stages!(stochasticprogram::StochasticProgram{N}, s::Integer) where N
    1 <= s <= N || error("Stage $s not in range 1 to $N.")
    for i = s:N
        remove_stage!(stochasticprogram, i)
    end
end

function _clear_model!(model::JuMP.Model)
    empty!(model.colLower)
    empty!(model.colUpper)
    empty!(model.colVal)
    empty!(model.colCat)
    empty!(model.colNames)
    empty!(model.colNamesIJulia)
    empty!(model.objDict)
    model.numCols = 0
    return nothing
end
function clear_parent!(scenarioproblems::ScenarioProblems)
    _clear_model!(scenarioproblems.parent)
    return nothing
end
function clear_parent!(scenarioproblems::DScenarioProblems)
    _clear_model!(scenarioproblems.parent)
    for w in workers()
        remotecall_fetch((sp) -> clear_parent!(fetch(sp)), w, scenarioproblems[w-1])
    end
    return nothing
end

function remove_scenarios!(stochasticprogram::StochasticProgram)
    remove_scenarios!(stochasticprogram.scenarioproblems)
    return nothing
end

function remove_subproblems!(stochasticprogram::StochasticProgram, s::Integer)
    remove_subproblems!(scenarioproblems(stochasticprogram, s))
    return nothing
end

function remove_subproblems!(stochasticprogram::StochasticProgram{2}, s::Integer = 2)
    remove_subproblems!(stochasticprogram.scenarioproblems)
    return nothing
end

function transfer_model!(dest::StochasticProgram, src::StochasticProgram)
    empty!(dest.generator)
    merge!(dest.generator, src.generator)
    return dest
end

function masterterms(scenarioproblems::ScenarioProblems, i::Integer)
    model = scenarioproblems.problems[i]
    parent = parentmodel(scenarioproblems)
    masterterms = Vector{Tuple{Int,Int,Float64}}()
    for (i,constr) in enumerate(model.linconstr)
        for (j,var) in enumerate(constr.terms.vars)
            if var.m == parent
                push!(masterterms,(i,var.col,-constr.terms.coeffs[j]))
            end
        end
    end
    return masterterms
end

function masterterms(scenarioproblems::DScenarioProblems, i::Integer)
    j = 0
    for w in workers()
        n = scenarioproblems.scenario_distribution[w-1]
        if i <= n+j
            return remotecall_fetch((sp,idx) -> masterterms(fetch(sp),idx), w, scenarioproblems[w-1], i-j)
        end
        j += n
    end
    throw(BoundsError(scenarioproblems, i))
end

function Base.copy(src::TwoStageStochasticProgram; procs = workers())
    dest = StochasticProgram(stage_parameters(src, 1), stage_parameters(src, 2), scenariotype(src); procs = procs)
    dest.spsolver.solver = src.spsolver.solver
    merge!(dest.generator, src.generator)
    return dest
end

function supports_zero(types::Vector, provided_def::Bool)
    for vartype in types
        if !hasmethod(zero, (Type{vartype}, ))
            !provided_def && @warn "Zero not defined for $vartype. Cannot generate zero function."
            return false
        end
    end
    return true
end

function supports_expected(types::Vector, provided_def::Bool)
    for vartype in types
        if !hasmethod(+, (vartype, vartype))
            !provided_def && @warn "Addition not defined for $vartype. Cannot generate expectation function."
            return false
        end
        if !hasmethod(*, (Float64, vartype)) || Base.code_typed(*, (Float64, vartype))[1].second != vartype
            !provided_def && @warn "Scalar multiplication with Float64 not defined for $vartype. Cannot generate expectation function."
            return false
        end
    end
    return true
end

problemcache(stochasticprogram::StochasticProgram) = stochasticprogram.problemcache
function get_problem(stochasticprogram::StochasticProgram, key::Symbol)
    haskey(stochasticprogram.problemcache, key)|| error("No $key in problem cache")
    return stochasticprogram.problemcache[key]
end
function get_stage_one(stochasticprogram::StochasticProgram)
    haskey(stochasticprogram.problemcache, :stage_1) || error("First-stage problem not generated.")
    return stochasticprogram.problemcache[:stage_1]
end
function get_stage(stochasticprogram::StochasticProgram, stage::Integer)
    stage_key = Symbol(:stage_, stage)
    haskey(stochasticprogram.problemcache, stage_key) || error("Stage problem $stage not generated.")
    return stochasticprogram.problemcache[stage_key]
end
function pick_solver(stochasticprogram::StochasticProgram, supplied_solver::SPSolverType)
    if supplied_solver isa JuMP.UnsetSolver
        return stochasticprogram.spsolver.solver
    end
    return supplied_solver
end

internal_solver(solver::MathProgBase.AbstractMathProgSolver) = solver

solverstr(solver::MathProgBase.AbstractMathProgSolver) = split(split(string(solver), "Solver")[1], ".")[2]

typename(dtype::UnionAll) = dtype.body.name.name
typename(dtype::DataType) = dtype.name.name

function set_decision!(stochasticprogram::StochasticProgram, x::AbstractVector)
    first_stage = get_stage_one(stochasticprogram)
    length(x) == first_stage.numCols || error("Incorrect length of given decision vector, has ", length(x), " should be ", first_stage.numCols)
    first_stage.colVal = copy(x)
    return nothing
end
function set_first_stage_redcosts!(stochasticprogram::StochasticProgram, μ::AbstractVector)
    first_stage = get_stage_one(stochasticprogram)
    length(μ) == first_stage.numCols || error("Incorrect length of given reduced costs, has ", length(μ), " should be ", first_stage.numCols)
    first_stage.redCosts = copy(μ)
    return nothing
end
function set_first_stage_duals!(stochasticprogram::StochasticProgram, λ::AbstractVector)
    first_stage = get_stage_one(stochasticprogram)
    length(λ) == length(first_stage.linconstr) || error("Incorrect length of given constraint duals, has ", length(μ), " should be ", first_stage.numCols)
    first_stage.linconstrDuals = copy(λ)
    return nothing
end

function add_subscript(src::AbstractString, subscript::Integer)
    return @sprintf("%s%s", src, unicode_subscript(subscript))
end
add_subscript(src::Symbol, subscript::Integer) = add_subscript(String(src), subscript)

function unicode_subscript(subscript::Integer)
    if subscript < 0
        error("$subscript is negative")
    end
    return join('₀'+d for d in reverse(digits(subscript)))
end
# ========================== #


text:
project"""
[deps]
Configurations="0.15"
"""


text:
"""
 WeightSpace represents a space that weights another space.
 Overload weight(S,x).
"""
abstract WeightSpace{S,T,DD,d} <: Space{T,DD,d}

domain(S::WeightSpace)=domain(S.space)

points(sp::WeightSpace,n)=points(sp.space,n)

immutable WeightSpacePlan{S,P,T,V}
    space::S
    plan::P
    points::Vector{T}
    weights::Vector{V}
end

for TYP in (:plan_transform,:plan_itransform)
    @eval function $TYP(S::WeightSpace,vals::Vector)
        pts=points(S,length(vals))
        WeightSpacePlan(S,$TYP(S.space,vals),pts,weight(S,pts))
    end
end

transform(sp::WeightSpace,vals::Vector,plan::WeightSpacePlan) =
    transform(sp.space,vals./(sp==plan.space?plan.weights:weight(sp,plan.points)),plan.plan)
itransform(sp::WeightSpace,cfs::Vector,plan::WeightSpacePlan) =
    itransform(sp.space,cfs,plan.plan).*(sp==plan.space?plan.weights:weight(sp,plan.points))

function evaluate(f::AbstractVector,S::WeightSpace,x...)
    tol=1.0E-14
    fv=Fun(f,S.space)(x...)
    if isa(fv,Number)&&abs(fv)<tol
        #TODO: Why this special case??
        zero(eltype(fv))
    else
        weight(S,x...).*fv
    end
end

include("divide_singularity.jl")
include("JacobiWeight.jl")
include("JacobiWeightOperators.jl")
include("JacobiWeightChebyshev.jl")
include("HeavisideSpace.jl")
include("DiracSpace.jl")
include("LogWeight.jl")
include("ExpWeight.jl")


text:
mutable struct TempField

    _NEWX_    :: AbstractArray{Float64, 2}
    _INTMX_    :: AbstractArray{Float64, 2}
    _ΔX_       :: AbstractArray{Float64, 2}
    _WKRSTΔX_   :: AbstractArray{Float64, 2}

    _CHKX_   :: AbstractArray{Float64, 2}
    _CHKX    :: AbstractArray{Float64, 1}

    _TMP_CHKX_   :: AbstractArray{Float64, 2} # Used to store the ∫ X dz before steping
    _TMP_CHKX    :: AbstractArray{Float64, 1}

    _TMP_SUBSTEP_BUDGET_   :: AbstractArray{Float64, 2} # Used to store substep tracer budget

    op_vdiff :: Union{AbstractArray, Nothing}
    sv       :: Union{Dict, Nothing}

    datastream   :: Union{Dict, Nothing}

    check_usage :: Dict

    function TempField(
        ev :: Env,
    )

        Nz, Nx, Ny = ev.Nz, ev.Nx, ev.Ny

        sT_pts = Nx * Ny * 1
        T_pts = Nx * Ny * Nz
        U_pts = Nx * Ny * Nz
        V_pts = Nx * (Ny+1) * Nz
        W_pts = Nx * Ny * (Nz+1)

        _NEWX_  = zeros(Float64, T_pts, 2)
        _INTMX_ = zeros(Float64, T_pts, 2)
        _ΔX_    = zeros(Float64, T_pts, 2)
        _WKRSTΔX_    = zeros(Float64, T_pts, 2)

        _CHKX_ = zeros(Float64, sT_pts, 2)
        _CHKX  = view(_CHKX_, :)

        _TMP_CHKX_ = zeros(Float64, sT_pts, 2)
        _TMP_CHKX  = view(_TMP_CHKX_, :)

        _TMP_SUBSTEP_BUDGET_ = zeros(Float64, sT_pts, 2)

        check_usage = Dict(
            :Ks_H_U => zeros(Nz, Nx, Ny),
            :Ks_H_V => zeros(Nz, Nx, Ny+1),
        )

        tmpfi = new(

            _NEWX_,
            _INTMX_,
            _ΔX_,
            _WKRSTΔX_,

            _CHKX_,
            _CHKX,

            _TMP_CHKX_,
            _TMP_CHKX,

            _TMP_SUBSTEP_BUDGET_,

            nothing,
            nothing,
            nothing,

            check_usage,
        )

        tmpfi.sv = getSugarView(ev, tmpfi)

        return tmpfi

    end

end

function getSugarView(
    ev :: Env,
    tmpfi :: TempField,
)

    Nx, Ny, Nz = ev.Nx, ev.Ny, ev.Nz

    return sv = Dict(
        :CHKTEMP  => reshape(view(tmpfi._CHKX_, :, 1), 1, Nx, Ny),
        :CHKSALT  => reshape(view(tmpfi._CHKX_, :, 2), 1, Nx, Ny),
        :NEWTEMP  => reshape(view(tmpfi._NEWX_, :, 1), Nz, Nx, Ny),
        :NEWSALT  => reshape(view(tmpfi._NEWX_, :, 2), Nz, Nx, Ny),
        :INTMTEMP => reshape(view(tmpfi._INTMX_, :, 1), Nz, Nx, Ny),
        :INTMSALT => reshape(view(tmpfi._INTMX_, :, 2), Nz, Nx, Ny),
    )

end

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
using CUDD
using CUDD.Compat.Test

@testset "addind vars" begin
    for i = 1:10
        manager = initialize_cudd()
        rnd = rand(1:100)
        for i = 1:rnd
            add_var(manager)
        end
        @test node_count(manager)-2 == rnd # the constant 0 and 1 nodes are created automatically
    end

    for i = 1:10
        manager = initialize_cudd()
        rnd = rand(i*10:i*10+9)
        indexed_var = add_ith_var(manager, rnd)
        @test read_index(indexed_var) == rnd
    end

    manager = initialize_cudd()
    for i = 1:10
        new_nonconst = add_var(manager)
        @test is_nonconst(new_nonconst) == 1
        new_nonconst = add_level_var(manager, i)
        @test is_nonconst(new_nonconst) == 1
        new_nonconst = add_ith_var(manager, i)
        @test is_nonconst(new_nonconst) == 1
        new_const = add_const(manager, i)
        @test is_const(new_const) == 1
    end
end

@testset "outputting files" begin
    manager = initialize_cudd()
    g = add_var(manager)
    @test output_dot(manager, g, "test.dot") == 1
    rm("test.dot")
    @test output_stats(manager, "stats") == 1
    rm("stats")
end

@testset "applying constant functions" begin

    manager = initialize_cudd()
    x1, x2 = rand(1:50), rand(51:100)
    f = add_const(manager, x1)
    ref(f)
    g = add_const(manager, x2)
    ref(g)
    res = add_apply(manager, add_plus_c, f, g)
    @test get_value(res) == x1+x2

    res = add_apply(manager, add_times_c, f, g)
    @test get_value(res) == x1*x2

    res = add_apply(manager, add_minus_c, f, g)
    @test get_value(res) == x1-x2

    res = add_apply(manager, add_divide_c, f, g)
    @test get_value(res) == x1/x2

    res = add_apply(manager, add_diff_c, f, g)
    @test get_value(res) == x1

    smaller = add_leq(manager, f, g)
    @test smaller == 1

    bigger = add_leq(manager, g, f)
    @test bigger == 0

    deref(f)
    deref(g)
    recursive_deref(manager, res)
end

@testset "addition testing" begin
    manager = initialize_cudd()
    f = add_ith_var(manager, 0)
    ref(f)
    for i = 1:100
        x = add_ith_var(manager, i)
        ref(x)
        tmp = add_apply(manager, add_plus_c, f, x)
        ref(tmp)
        recursive_deref(manager, x)
        recursive_deref(manager, f)
        f = tmp
    end
    assignment = Int32[]
    res = 0
    for i = 0:100
        generated = i%2
        push!(assignment, generated)
        res += generated
    end
    cudd_res = evaluate(manager, f, assignment)
    value = convert(Int, get_value(cudd_res))
    @test value == res

    for i = 1:10
        assignment = Int32[]
        res = 0
        for i = 1:101
            generated = rand(0:1)
            push!(assignment, generated)
            res += generated
        end
        cudd_res = evaluate(manager, f, assignment)
        value = convert(Int, get_value(cudd_res))
        @test value == res
    end
end

@testset "multiplication testing" begin
    manager = initialize_cudd()
    f = add_ith_var(manager, 1)
    ref(f)
    for i = 2:101
        x = add_ith_var(manager, i)
        ref(x)
        tmp = add_apply(manager, add_times_c, f, x)
        ref(tmp)
        recursive_deref(manager, x)
        recursive_deref(manager, f)
        f = tmp
    end
    assignment = Int32[]
    res = 0
    for i = 1:101
        generated = i%2
        push!(assignment, generated)
        res *= generated
    end
    cudd_res = evaluate(manager, f, assignment)
    value = convert(Int, get_value(cudd_res))
    @test value == res
end


text:
export
    Instance,
    VariableInstance,
    PlaceholderInstance,
    Env,
    Runtime,
    InstantRuntime,
    DynamicRuntime,
    collect_messages,
    current_instance,
    current_time,
    get_all_instances,
    get_belief,
    get_placeholder_beliefs,
    get_definition,
    get_env,
    get_evidence,
    get_instance,
    get_intervention,
    get_message,
    get_name,
    get_network,
    get_range,
    get_sfunc,
    get_state,
    get_time,
    get_value,
    get_node,
    has_belief,
    has_evidence,
    has_instance,
    has_intervention,
    has_previous_instance,
    has_range,
    has_state,
    has_value,
    previous_instance,
    rng,
    clear_state!,
    concat!,
    delete_evidence!,
    delete_state!,
    delete_value!,
    distribute_messages!,
    instantiate!,
    ensure_all!,
    ensure_instance!,
    post_belief!,
    post_evidence!,
    post_intervention!,
    latest_instance_before,
    remove_messages!,
    set_message!,
    set_range!,
    set_state!,
    set_time!,
    set_value!

using DataStructures
using Dates
using Random

using ..Models

import Base.isless

"""
    Instance{O}

An abstract base type for variable and placeholder instances.

This is instance can have values associated with it in the runtime.
O is the output type of the node.
"""
abstract type Instance{O} end

"""
    VariableInstance

An instance of variable node at `time`.   `sf` is an sfunc generated from
the variable's model.
"""
struct VariableInstance{O} <: Instance{O}
    node :: Variable{I,J,O} where {I,J}
    sf :: Union{SFunc{I,O}, SFunc{J,O}} where {I,J}
    time :: Any
end

"""
    PlaceholderInstance

An instance of placeholder `node` at `time`.  This instance has no sfunc, but can still take
values in the runtime.
"""
struct PlaceholderInstance{O} <: Instance{O}
    node :: Placeholder{O}
    time :: Any
end

"""
    get_sfunc(i::VariableInstance)::SFunc

Get the instance's sfunc.
"""
get_sfunc(i::VariableInstance)::SFunc = i.sf

"""
    get_node(i::Instance)::Node

Get the instance's node, whether it is a placeholder or a variable.
"""
get_node(i::Instance)::Node = i.node

"""
    get_time(i::Instance)

Get the instance's time.
"""
get_time(i::Instance) = i.time

"""
    get_name(i::Instance)::Symbol

Get the name of the instance's node.
"""
get_name(i::Instance):: Symbol = get_node(i).name

"""
    get_definition(i::VariableInstance)::D where {D<:ValueTyped}

Get the instance's variable's underlying model.
"""
get_definition(i::VariableInstance)::D where {D<:ValueTyped} = get_variable(i).model

"""
    get_model(i::VariableInstance)::D where {D<:ValueTyped}

Get the instance's variable's underlying model.
"""
get_model(i::VariableInstance)::D where {D<:ValueTyped} = get_definition(i)

output_type(i::Instance) = output_type(get_node(i))

"""
    struct Env

Holds all external state of a Runtime.  The `Env` supports the following methods:

```
get_state(env::Env) :: Dict{Symbol, Any}
has_state(env::Env, key::Symbol) :: Bool
get_state(env::Env, key::Symbol)
set_state!(env::Env, key::Symbol, value)
delete_state!(env::Env, key::Symbol)
clear_state!(env::Env)
clone(env::Env)
```
"""
struct Env
    state::Dict{Symbol, Any}
    Env() = new(Dict{Symbol, Any}())
end

get_state(env::Env) :: Dict{Symbol, Any} = env.state
has_state(env::Env, key::Symbol) :: Bool = haskey(get_state(env), key)
function get_state(env::Env, key::Symbol)
    has_state(env, key) || ArgumentError("symbol :$key does not exist in the environment")
    get_state(env)[key]
end

set_state!(env::Env, key::Symbol, value) = get_state(env)[key] = value
delete_state!(env::Env, key::Symbol) = pop!(get_state(env), key, nothing)
clear_state!(env::Env) = empty!(get_state(env))

clone(env::Env) :: Env = deepcopy(env)

# TODO what should happen if a global and local key conflicts?
function concat!(local_env::Env, global_env::Env) :: Env
    for (k,v) in get_state(clone(global_env))
        set_state!(local_env, k, v)
    end
    local_env
end

"""
    abstract type Runtime

A struct that contains the state of the compute graph. This code makes the assumption
that values are associated with instances but messages are passed between
variables and applied to the relevant instances later.  It has to be this way
because the receiving instance might not exist at the time the message is sent.
"""
abstract type Runtime end

"""
    struct InstantRuntime <: Runtime

A runtime that represents a network whose variables take on a single instance.  As
a convenience, the following methods create an `InstantRuntime`:

```
Runtime()
Runtime(net :: InstantNetwork)
Runtime(net :: DynamicNetwork)
```
"""
struct InstantRuntime <: Runtime
    env :: Env
    name :: Symbol
    network :: InstantNetwork
    instances :: Dict{Node, Instance}
    values :: Dict{Tuple{Instance, Symbol}, Any}
    messages :: Dict{Tuple{Node, Symbol}, Dict{Node, Any}}
end

"""
    struct DynamicRuntime <: Runtime

A runtime that represents a network whose variables take on many instances at different times `T`.
"""
struct DynamicRuntime{T} <: Runtime
    env :: Env
    name :: Symbol
    network :: DynamicNetwork
    instances :: Dict{Node,
        SortedDict{T, Instance, Base.Order.ReverseOrdering}}
    values :: Dict{Tuple{Instance, Symbol}, Any}
    messages :: Dict{Tuple{Node, Symbol}, Dict{Node, Any}}
end

Runtime() = InstantRuntime(InstantNetwork(Variable[], Placeholder[], Placeholder[],
    VariableGraph()))
Runtime(net :: InstantNetwork) = InstantRuntime(Env(), gensym(), net,
    Dict(), Dict(), Dict())
Runtime(net :: DynamicNetwork) = Runtime(net, 0)
function Runtime(net :: DynamicNetwork, time::T) where {T}
    rt = DynamicRuntime{T}(Env(), gensym(), net, Dict(), Dict(), Dict())
    set_time!(rt, time)
    return rt
end

function rng(r::Runtime)
    return Random.GLOBAL_RNG
end

get_env(runtime::Runtime) :: Env = runtime.env
get_name(runtime::Runtime) :: Symbol = runtime.name
get_network(runtime::Runtime) :: Network = runtime.network
get_nodes(runtime::Runtime) :: Set{Node} = get_nodes(get_network(runtime))
function get_node(runtime::Runtime, name::Symbol) :: Union{Placeholder, Variable, Nothing}
    get_node(get_network(runtime), name)
end

#=
    state functions
=#

get_state(runtime::Runtime) :: Dict{Symbol, Any} = get_state(get_env(runtime))
has_state(runtime::Runtime, key::Symbol) :: Bool = has_state(get_env(runtime), key)
get_state(runtime::Runtime, key::Symbol) = get_state(get_env(runtime), key)
set_state!(runtime::Runtime, key::Symbol, value) = set_state!(get_env(runtime), key, value)
delete_state!(runtime::Runtime, key::Symbol) = delete_state!(get_env(runtime), key)
clear_state!(runtime::Runtime) = empty!(get_state(runtime))

#=
    Managing Runtime time
=#
const TIME = :__simulated_time__

"""
    current_time(runtime::Runtime) -> T

Returns the currently set time of the given `Runtime`
"""
function current_time(runtime::DynamicRuntime{T}) :: T where {T}
    return get_state(runtime, TIME)
end

"""
    set_time(runtime::Runtime{T}, newtime::T) -> T

Sets the current time for the given `Runtime`
"""
set_time!(runtime::DynamicRuntime{T}, newtime::T) where {T} = set_state!(runtime, TIME, newtime)

#=
    Instantiating instances
=#

isless(::Nothing, ::Nothing) = false

"""
    instantiate!(runtime::InstantRuntime,variable::Variable,time = 0)
    instantiate!(runtime::InstantRuntime,placeholder::Placeholder,time = 0)
    instantiate!(runtime::DynamicRuntime{T}, node::Node,time::T = current_time(runtime))::Instance where {T}

Instantiate and return an instance for the given runtime at the given time; the default
time is the current time of the runtime in unix time (an `Int64`). For an `InstantRuntime`,
there is only a single instance for each variable.
"""
function instantiate!(runtime::InstantRuntime,variable::Variable,time = 0)
        # time argument is provided for uniformity but ignored
    return get!(runtime.instances, variable,
        VariableInstance(variable, make_initial(variable.model, time), time))
end

function instantiate!(runtime::InstantRuntime,placeholder::Placeholder,time = 0)
    return get!(runtime.instances, placeholder, PlaceholderInstance(placeholder, time))
end

function instantiate!(runtime::DynamicRuntime{T}, node::Node,time::T = current_time(runtime))::Instance where {T}
    if haskey(runtime.instances, node)
        curr = (runtime.instances[node])
        @assert !in(time, keys(curr)) "variable $(variable.name) at time $(time) is already instantiated"
        parents = get_transition_parents(runtime.network, node)
        parenttimes = Vector{T}()
        for p in parents
            parinst = latest_instance_before(runtime, p, time, p != node)
            if isnothing(parinst)
                error("In instantiate! for ", variable.name, ": parent ", p.name, " not instantiated at time ", time)
            end
            push!(parenttimes, get_time(parinst))
        end
        parenttimes = tuple(parenttimes...)
        if isa(node, Variable)
            sf = make_transition(node.model, parenttimes, time)
            inst = VariableInstance(node, sf, time)
        else
            inst = PlaceholderInstance(node, time)
        end
        curr[time] = inst
        return inst
    else
        curr =
            SortedDict{T, Instance, Base.Order.ReverseOrdering}
                (Base.Order.ForwardOrdering())
        local inst::Instance
        if isa(node, Variable)
            sf = make_initial(node.model, time)
            inst = VariableInstance(node, sf, time)
        else
            inst = PlaceholderInstance(node, time)
        end
        instvec = Pair{T, Instance}[time => inst]
        runtime.instances[node] =
            SortedDict{T, Instance, Base.Order.ReverseOrdering}(Base.Order.ReverseOrdering(), instvec)
        return inst
    end
end

"""
    ensure_all!(runtime::InstantRuntime, time=0) :: Dict{Symbol, Instance}
    ensure_all!(runtime::DynamicRuntime, time = current_time(runtime)) :: Dict{Symbol, Instance}

Instantiate all the variables in the network at the given time and returns them as a
dict from variable names to their corresponding instance; the default time is the current
time of the runtime in unix time.
"""
function ensure_all!(runtime::InstantRuntime, time=0) :: Dict{Symbol, Instance}
    insts = Dict{Symbol, Instance}()
    for node in get_nodes(runtime)
        insts[node.name] = instantiate!(runtime, node)
    end
    return insts
end

function ensure_all!(runtime::DynamicRuntime, time = current_time(runtime)) :: Dict{Symbol, Instance}
    insts = Dict{Symbol, Instance}()
    for node in get_nodes(runtime)
        insts[node.name] = ensure_instance!(runtime, node, time)
    end
    return insts
end

"""
    ensure_instance!(runtime::Runtime{T}, node::Node{O}, time::T = current_time(runtime))::Instance{O} where {O,T}

Returns an instance for the given variable at the given time, either by using an existing one or creating a new one.
"""
function ensure_instance!(runtime::Runtime, node::Node{O}, time = current_time(runtime))::Instance{O} where O
    if has_instance(runtime, node, time)
        return get_instance(runtime, node, time)
    else
        return instantiate!(runtime, node, time)
    end
end

"""
    current_instance(runtime::InstantRuntime, node::Node)
    current_instance(runtime::DynamicRuntime, node::Node)

Returns the current (last) instance for the given runtime and node; this method will
throw an exception if there is no current instance for the given node
"""
current_instance(runtime::InstantRuntime, node::Node) = runtime.instances[node]

function current_instance(runtime::DynamicRuntime, node::Node) :: Instance
    first(runtime.instances[node]).second
end

"""
    previous_instance(runtime::DynamicRuntime, node::Node)

Returns the previous instance for the given runtime and node.  This will throw
and exception if there is no previous instance.
"""
function previous_instance(runtime::DynamicRuntime, node::Node)
    insts = runtime.instances[node]
    key = keys[insts][2]
    return insts[key]
end

"""
    get_all_instances(runtime::DynamicRuntime, variable::Variable)
    get_all_instances(runtime::InstantRuntime, variable::Variable)

Returns all instances in the `runtime` for the given variable, in order, as a
`Vector{Instance}`.
"""
function get_all_instances(runtime::DynamicRuntime, variable::Variable)
    has_instance(runtime, variable) ? collect(values(runtime.instances[variable])) : Instance[]
end

function get_all_instances(runtime::InstantRuntime, variable::Variable)
    has_instance(runtime, variable) ? [runtime.instances[variable]] : Instance[]
end

"""
    function get_instance(runtime::DynamicRuntime{T}, node::Node, t::T)::Instance

Returns instance for the given variable at time `T`; throws an error if no
such instance exists.
"""
function get_instance(runtime::DynamicRuntime{T}, node::Node, t::T) :: Instance where {T}
    return runtime.instances[node][t]
end

"""
    latest_instance_before(runtime::DynamicRuntime{T}, node::Node, t::T, allow_equal::Bool) :: Union{Instance, Nothing} where T

Return the latest instance of the node before the given time.  The `allow_equal` flag indicates
whether instances at a time equal to `t`` are allowed.

If there is no such instance, returns `nothing`.

"""
function latest_instance_before(runtime::DynamicRuntime{T}, node::Node, t::T, allow_equal::Bool) :: Union{Instance, Nothing} where T
    if !(node in keys(runtime.instances))
        return nothing
    end
    smd = runtime.instances[node]
    for u in keys(smd)
        if allow_equal ? u <= t : u < t
            return smd[u]
        end
    end
    return nothing
end

"""
    has_instance(runtime::DynamicRuntime, node::Node, time = current_time(runtime))
    has_instance(runtime::InstantRuntime, node::Node)

Returns true if the given node has an instance in the given runtime at a time greater than or equal to the given time.
"""
function has_instance(runtime::DynamicRuntime, node::Node, time = current_time(runtime))
    haskey(runtime.instances, node) && time in keys(runtime.instances[node])
end

has_instance(runtime::InstantRuntime, node::Node) = haskey(runtime.instances, node)

"""
    has_previous_instance(runtime::DynamicRuntime, node::Node)

Checks if the specified node has an instance prior to the current one.
"""
function has_previous_instance(runtime::DynamicRuntime, node::Node)
    haskey(runtime.instances, node) && length(keys(runtime.instances, node)) > 1
end

#=
    Setting and getting values associated with instances
=#

"""
    set_value!(runtime::Runtime, instance::Instance, key::Symbol, value)

Set the value on an instance for the given key
"""
function set_value!(runtime::Runtime, instance::Instance, key::Symbol, value)
    runtime.values[(instance, key)] = value
end

"""
    get_value(runtime::Runtime, instance::Instance, key::Symbol)

Get the value on an instance for the given key; this will throw an exception
if the instance does not contain the given key
"""
function get_value(runtime::Runtime, instance::Instance, key::Symbol)
    result = runtime.values[(instance, key)]
    return result
end

"""
    delete_value!(runtime::Runtime{T}, instance::Instance, key::Symbol) where {T}

Deletes the mapping for the given instance and key in the runtime and returns it
"""
function delete_value!(runtime::Runtime, instance::Instance, key::Symbol)
    return pop!(runtime.values, (instance, key), nothing)
end

function has_value(runtime::Runtime, instance::Instance, key::Symbol)
    (instance, key) in keys(runtime.values)
end

#=
    Sending and receiving messages between variables
=#

function set_message!(runtime::Runtime, sender::Node, recipient::Node, key::Symbol, value)
    msgs = get(runtime.messages, (sender, key), Dict())
    msgs[recipient] = value
    runtime.messages[(sender, key)] = msgs
end

function get_message(runtime::Runtime, sender::Node, recipient::Node, key::Symbol)
    runtime.messages[(sender, key)][recipient]
end

function distribute_messages!(runtime::Runtime, sender::Node, recipients::Vector{Node},
     key::Symbol, values)
    msgs = get(runtime.messages, (sender, key), Dict())
    for (rcp, val) in zip(recipients, values)
        msgs[rcp] = val
    end
    runtime.messages[(sender, key)] = msgs
end

function collect_messages(runtime::Runtime, senders::Vector{Node}, recipient::Node,
        key::Symbol)
    [get_message(runtime, sdr, recipient, key) for sdr in senders]
end

function remove_messages!(runtime::Runtime, sender::Node, key::Symbol)
    delete!(runtime.messages, (sender, key))
end

#=
    Invoking an operation on an instance
=#

function _unwrap(p)
    if p isa Tuple && length(p) == 1
        p[1]
    else
        p
    end
end

"""
    RANGE

The constant key used to store the range of a specific variable instance
"""
const RANGE = :__range__

"""
    set_range!(runtime::Runtime, inst::Instance{O}, range::Vector{<:O}, depth::Int = 1) where O

Sets the range value for the given instance. Defaults to depth of 1.
"""
function set_range!(runtime::Runtime, inst::Instance{O}, range::Vector{<:O}, depth::Int = 1) where O
    if has_value(runtime, inst, RANGE)
        curr = get_value(runtime, inst, RANGE)
        s = Tuple{Int, Vector{O}}[]
        i = 1
        while i <= length(curr)
            pair = curr[i]
            d = pair[1]
            if d > depth
                push!(s, pair)
                i += 1
            end
            push!(s, (depth, range))
            i = d == depth ? i+1 : i
            for j = i:length(curr)
                push!(s, curr[j])
            end
            set_value!(runtime, inst, RANGE, s)
            return
        end
        push!(s, (depth, range))
        set_value!(runtime, inst, RANGE, s)
    else
        set_value!(runtime, inst, RANGE, [(depth, range)])
    end
end

"""
    get_range(runtime::Runtime, inst::Instance, depth = max_value(Int))

Returns the range value for the given instance; this will return
`nothing` if no range has been set.

The depth specifies the maximum depth of range desired.
"""
function get_range(runtime::Runtime, inst::Instance, depth = typemax(Int))
    has_range(runtime, inst, depth) || return nothing
    rng = get_value(runtime, inst, RANGE)
    for i in 1:length(rng)
        (d,r) = rng[i]
        if d <= depth
            return r
        end
    end
    return nothing
end

function has_range(runtime::Runtime, inst::Instance, depth::Int = typemax(Int))
    has_value(runtime, inst, RANGE) || return false
    r = get_value(runtime, inst, RANGE)
    (d,_) = r[length(r)]
    return d <= depth
end

"""
    BELIEF

The constant key used to store belief for a specific variable instance
"""
const BELIEF = :__belief__

"""
    post_belief!(runtime::Runtime, inst::Instance, belief)

Posts the given evidence for the given instance.
"""
function post_belief!(runtime::Runtime, inst::Instance, belief)
    set_value!(runtime, inst, BELIEF, belief)
end

"""
    get_belief(runtime::Runtime, inst::Instance)

Returns the last posted belief for the given instance; this will return
`nothing` if no belief has been posted
"""
function get_belief(runtime::Runtime, inst::Instance)
    has_value(runtime, inst, BELIEF) || return nothing
    return get_value(runtime, inst, BELIEF)
end

has_belief(runtime::Runtime, inst::Instance) = has_value(runtime, inst, BELIEF)

function get_placeholder_beliefs(runtime::Runtime)::Dict{Symbol,Dist}
    result = Dict{Symbol,Dist}()
    for ph in get_placeholders(get_network(runtime))
        i = current_instance(runtime, ph)
        if has_belief(runtime, i)
            result[ph.name] = get_belief(runtime, i)
        end
    end
    return result
end
"""
    EVIDENCE

The constant key used to store evidence for a specific variable instance
"""
const EVIDENCE = :__evidence__

"""
    post_evidence!(runtime::Runtime, inst::Instance, evidence)

Posts the given evidence for the given instance.
"""
function post_evidence!(runtime::Runtime, inst::Instance, evidence::Score)
    set_value!(runtime, inst, EVIDENCE, evidence)
end

"""
    get_evidence(runtime::Runtime, inst::Instance)

Returns the last posted evidence for the given instance; this will return
`nothing` if no evidence has been posted
"""
function get_evidence(runtime::Runtime, inst::Instance)
    has_value(runtime, inst, EVIDENCE) || return nothing
    return get_value(runtime, inst, EVIDENCE)
end

function delete_evidence!(runtime::Runtime, inst::Instance)
    if has_value(runtime, inst, EVIDENCE)
        return delete_value!(runtime, inst, EVIDENCE)
    end
end

has_evidence(runtime::Runtime, inst::Instance) = has_value(runtime, inst, EVIDENCE)

"""
    INTERVENTION

The constant key used to store interventions for a specific variable instance
"""
const INTERVENTION = :__intervention__

"""
    post_intervention!(runtime::Runtime, inst::Instance, intervention::Dist)

Posts the given intervention for the given instance.
"""
function post_intervention!(runtime::Runtime, inst::Instance, intervention::Dist)
    set_value!(runtime, inst, INTERVENTION, intervention)
end

"""
    get_intervention(runtime::Runtime, inst::Instance)

Returns the last posted intervention for the given instance; this will return
`nothing` if no intervention has been posted
"""
function get_intervention(runtime::Runtime, inst::Instance)
    has_value(runtime, inst, INTERVENTION) || return nothing
    return get_value(runtime, inst, INTERVENTION)
end

has_intervention(runtime::Runtime, inst::Instance) = has_value(runtime, inst, INTERVENTION)

#=
    Define the interface for Algorithms
=#

"""
    initialize_algorithm(alg_fun::Function, runtime::Runtime)

Initializes the given algorithm's [`Runtime`](@ref).  By default, this does nothing.
"""
initialize_algorithm(alg_fun::Function, runtime::Runtime) = nothing

"""
    start_algorithm(alg_fun::Function, runtime::Runtime, args...)

Starts the algorithm in the current thread
"""
start_algorithm(alg_fun::Function, runtime::Runtime, args...) = nothing

# TODO implement ; should take a list of instance/evidence pairs
# (or rather, variable-name/time/evidence tuples)
function post_evidence!(alg_fun::Function, runtime::Runtime, inst::Instance, evidence)
end