text:
using LinearAlgebra
using SparseArrays
using Plots
using ForwardDiff

uexact(x) = log(2+sin(pi*x))
dudx_exact(x) = ForwardDiff.derivative(uexact,x)
f(x) = -ForwardDiff.derivative(dudx_exact,x) # -d2u/d2x = f
α = uexact(-1)
β = uexact(1) - dudx_exact(1)

function solve(m,f,α,β)
    x = LinRange(-1,1,m+2) # add 2 endpoints
    xint = x[2:end]
    h = x[2]-x[1]
    A = (1/h^2) * spdiagm(0=>2*ones(m+1),-1=>-ones(m),1=>-ones(m))
    A[end,end] = (1-h)/h^2
    b = f.(xint)
    b[1] += α/h^2
    b[end] = .5*f(x[end]) - β/h
    u = A\b
    return u,x,xint,h
end

mvec = 2 .^ (2:9)
hvec = zeros(length(mvec))
err_inf = zeros(length(mvec))
for (i,m) in enumerate(mvec)
    u,x,xint,h = solve(m,f,α,β)
    hvec[i] = h
    err_inf[i] = maximum(@. abs(uexact(xint) - u))
end

plot(hvec,err_inf,marker=:dot,label="Max error")
plot!(hvec,hvec.^2,linestyle=:dash,label="O(h^2)")
plot!(xaxis=:log,yaxis=:log,legend=:bottomright)


text:
function f(x::Vector)
  (309.0 - 5.0 * x[1])^2 + (17.0 - x[2])^2
end
function g!(x::Vector, storage::Vector)
  storage[1] = -10.0 * (309.0 - 5.0 * x[1])
  storage[2] = -2.0 * (17.0 - x[2])
end

d = DifferentiableFunction(f, g!)

initial_x = [10.0, 10.0]
m = 10
store_trace, show_trace = false, false

results = Optim.l_bfgs(d, initial_x)
@assert length(results.trace.states) == 0
@assert results.converged
@assert norm(results.minimum - [309.0 / 5.0, 17.0]) < 0.01

ext: jl


text:
# GenSPGL
# spgOptions

export spgOptions

"""
GenSPGL input options composite type.

"""
mutable struct spgOptions
    fid::Integer
    verbosity::Integer
    iterations::Integer
    nPrevVals::Integer
    bpTol::Number
    lsTol::Number
    optTol::Number
    decTol::Number
    stepMin::Number
    stepMax::Number
    rootMethod::Integer
    activeSetIt::Number
    subspaceMin::Bool
    iscomplex::Nullable{Bool}
    maxMatvec::Number
    weights::AbstractArray
    quitPareto::Bool
    minPareto::Integer
    lineSrchIt::Integer
    feasSrchIt::Integer
    ignorePErr::Bool
    project::Function
    primal_norm::Function
    dual_norm::Function
    funPenalty::Function
    funCompositeR::Function
    proxy::Bool # Appropriate type?
    linear::Bool # ^
    restore::Bool # ^^
end

"""
A placeholder in the default spgl settings until all
functionality is ported over.

"""
function dummyfun()
    println("This doesn't do anything")
end

"""
#########################################################################################\n
spgOptions Outer Constructor\n
Arguments not specified will be set to their default values shown below.\n
#########################################################################################\n

    spgOptions(;fid = 1,
                    verbosity = 1,
                    iterations = 100000,
                    nPrevVals = 10,
                    bpTol = 1e-6,
                    lsTol = 1e-6,
                    optTol = 1e-4,
                    decTol = 1e-4,
                    stepMin = 1e-16,
                    stepMax = 1e+5,
                    rootMethod = 2,
                    activeSetIt = Inf,
                    subspaceMin = false,
                    iscomplex = Null,
                    maxMatvec = Inf,
                    weights = [1],
                    quitPareto = false,
                    minPareto = 3,
                    lineSrchIt = 1,
                    feasSrchIt = 10000,
                    ignorePErr = false,
                    project = dummyfun,
                    primal_norm = dummyfun,
                    dual_norm = NormL1_dual,
                    funPenalty = funLS,
                    funCompositeR = funCompR1,
                    proxy = false,
                    linear = false,
                    restore = false)

# Examples
```julia
# Default Options
julia> opts = spgOptions()

# Make optional changes
julia> opts = spgOptions(verbosity = 0, iterations = 20000)
```

Arguments:\n
    |___Field____||__Type__||_____________________Description____________________________|\n
    fid          |Integer|   ..... File ID for output\n
    verbosity    |Integer|   ..... Verbosity level, 0 = Silent, 1 = Loud, 2 = Debug\n
    iterations   |Integer|   ..... Max number of iterations\n
    nPrevVals    |Integer|   ..... Number previous func values for linesearch\n
    bpTol        |Number|    ..... Tolerance for basis pursuit solution \n
    lsTol        |Number|    ..... Least-squares optimality tolerance\n
    optTol       |Number|    ..... Optimality tolerance\n
    decTol       |Number|    ..... Reqd rel. change in primal obj. for Newton\n
    stepMin      |Number|    ..... Minimum spectral step\n
    stepMax      |Number|    ..... Maximum spectral step\n
    rootMethod   |Integer|   ..... Root finding method: 2=quad,1=linear (not used).\n
    activeSetIt, |Number|    ..... Exit with EXIT_ACTIVE_SET if nnz same for # its.\n
    subspaceMin, |Bool|      ..... Use subspace minimization\n
    iscomplex    |Nullable{Bool}|      ..... Flag set to indicate complex problem\n
    maxMatvec    |Number|    ..... Maximum matrix-vector multiplies allowed\n
    weights      |AbsrtactArray|    ..... Weights W in ||Wx||_1\n
    quitPareto   |Bool|      ..... Exits when pareto curve is reached\n
    minPareto    |Integer|   ..... If quitPareto is on, the minimum number of iterations
                                    before checking for quitPareto conditions\n
    lineSrchIt   |Integer|   ..... Maximum number of line search iterations for
                                    spgLineCurvy, originally 10\n
    feasSrchIt   |Integer|   ..... Maximum number of feasible direction line search
                                    iteraitons, originally 10\n
    ignorePErr   |Bool|      ..... Ignores projections error by issuing a warning instead
                                    of an error\n
    project      |Function|  ..... Projection function handle\n
    primal_norm  |Function|  ..... Primal Norm function handle\n
    dual_norm    |Function|  ..... Dual Norm function handle\n
    funPenalty   |Function|  ..... Penalty function handle\n
    proxy        |Bool|      ..... Advanced option that computes pareto curve in a user
                                    specified way. \n
    linear       |Bool|      ..... Advanced option that allows you to declare input
                                    functions to be linear\n
    restore      |Bool|      ..... Whether to restore best previous answer. for large
                                    problems, dont want to do this.\n
"""
function spgOptions(;fid::Integer = 1,
                        verbosity::Integer = 1,
                        iterations::Integer = 100000,
                        nPrevVals::Integer = 3,
                        bpTol::Number = 1e-6,
                        lsTol::Number = 1e-6,
                        optTol::Number = 1e-4,
                        decTol::Number = 1e-4,
                        stepMin::Number = 1e-16,
                        stepMax::Number = 1e+5,
                        rootMethod::Integer = 2,
                        activeSetIt::Number = Inf,
                        subspaceMin::Bool = false,
                        iscomplex::Nullable{Bool} = Nullable{Bool}(),
                        maxMatvec::Number = Inf,
                        weights::AbstractArray = [1],
                        quitPareto = false,
                        minPareto::Integer = 3,
                        lineSrchIt::Integer = 1,
                        feasSrchIt::Integer = 10000,
                        ignorePErr = false,
                        project::Function = GenSPGL.NormL1_project,
                        primal_norm::Function = GenSPGL.NormL1_primal,
                        dual_norm::Function = GenSPGL.NormL1_dual,
                        funPenalty::Function = GenSPGL.funLS,
                        funCompositeR::Function = GenSPGL.funCompR1,
                        proxy::Bool = false,
                        linear::Bool = false,
                        restore::Bool = false)

    opts = spgOptions(fid,
                        verbosity,
                        iterations,
                        nPrevVals,
                        bpTol,
                        lsTol,
                        optTol,
                        decTol,
                        stepMin,
                        stepMax,
                        rootMethod,
                        activeSetIt,
                        subspaceMin,
                        iscomplex,
                        maxMatvec,
                        weights,
                        quitPareto,
                        minPareto,
                        lineSrchIt,
                        feasSrchIt,
                        ignorePErr,
                        project,
                        primal_norm,
                        dual_norm,
                        funPenalty,
                        funCompositeR,
                        proxy,
                        linear,
                        restore)

    return opts
end

#DEVNOTE# Think of adding a method to the constructor so that non-nullable bools
        # could be passed as an optional input


text:
# Chemfiles.jl, a modern library for chemistry file reading and writing
# Copyright (C) Guillaume Fraux and contributors -- BSD license

export ChemfilesError

type ChemfilesError <: Exception
    message::AbstractString
end

Base.show(io::IO, e::ChemfilesError) = show(io, "Chemfiles error: $(e.message)")

"""
If there was an error try to get it and throw a ChemfilesError with it.
Show \"Unknown error\" in case it can't get the last one.
"""
function _check(result::Integer, message = "Unknown error")
    if result != 0
        str = message
        try
            str = last_error()
        end
        throw(ChemfilesError(str))
    end
    return nothing
end

"""
If there was an error try to get it and throw a ChemfilesError with it.
"""
function _check(result::Ptr)
    if Int(result) == 0
        throw(ChemfilesError(last_error()))
    end
    return nothing
end

"""
Get the last error message from the chemfiles runtime.
"""
function last_error()
    unsafe_string(lib.chfl_last_error())
end

"""
Clear any error message stored by the chemfiles runtime.
"""
function clear_errors()
    _check(lib.chfl_clear_errors())
end

"""
Define the warning callback.
"""
function _warning_callback(message::String)
    warn("[chemfiles] ", message)
end

"""
Set the global warning ``callback`` to be used for each warning event.

The ``callback`` function must take a ``String`` and return nothing.
"""
function set_warning_callback(callback::Function)
    function _cb_adaptor(message)
        try
            callback(unsafe_string(message))
        catch error
            warn("caught $error in warning callback")
        end
    end
    cb = cfunction(_cb_adaptor, Void, (Ptr{UInt8},))
    _check(lib.chfl_set_warning_callback(cb))
end

set_warning_callback(_warning_callback)

"""
Read configuration data from the file at `path`.

By default, chemfiles reads configuration from any file name `.chemfilesrc`
in the current directory or any parent directory. This function can be used
to add data from another configuration file.

This function will fail if there is no file at `path`, or if the file is
incorectly formatted. Data from the new configuration file will overwrite
any existing data.
"""
function add_configuration(path)
    lib.chfl_add_configuration(pointer(path))
end

"""
"""
function _strip_null(string)
    for i in 1:length(string)
        if string[i] == '\0'
            return string[1:i-1]
        end
    end
    throw(ChemfilesError("A C string is not NULL terminated"))
end

"""
"""
function _call_with_growing_buffer(callback::Function, initial_size=64)
    function _buffer_was_big_enough(buffer)
        if length(buffer) < 2
            return false
        else
            return buffer[end-2] == '\0'
        end
    end

    size = initial_size
    buffer = repeat("\0", size)

    callback(pointer(buffer), UInt64(size))
    while !_buffer_was_big_enough(buffer)
        # Grow the buffer and retry
        size *= 2
        buffer = repeat("\0", size)
        callback(pointer(buffer), UInt64(size))
    end

    return _strip_null(buffer)
end


text:
abstract type AbstractParallelogram <:AbstractLattice end

export Parallelogram, bravais_size, vec_a, vec_b, rectlattice

# square lattice
function Base.getindex(lt::AbstractParallelogram, i::Int)
    lt[CartesianIndices(bravais_size(lt))[i].I...]
end
Base.lastindex(lt::AbstractParallelogram, i::Int) = size(lt, i)
bravais_size(sq::AbstractParallelogram) = size(sq)
bravais_size(sq::AbstractParallelogram, i::Int) = size(sq, i)
function Base.size(sq::AbstractParallelogram, i::Int)
    size(sq)[i]
end
unit(lt::AbstractParallelogram) = 1/(maximum(size(lt, 1) .* vec_a(lt) .+ size(lt, 2) .* vec_b(lt)))
vertices(sq::AbstractParallelogram) = 1:size(sq, 1)*size(sq, 2)

function Base.getindex(lt::AbstractParallelogram, i::Real, j::Real)
    a, b = vec_a(lt), vec_b(lt)
    step = unit(lt)
    res = (i-0.5)*step .* a .+ (j-0.5)*step .* b
    (res[1], res[2])
end

struct Parallelogram <: AbstractParallelogram
    Nx::Int
    Ny::Int
    a::Vector{Float64}
    b::Vector{Float64}
end

rectlattice(Nx::Int, Ny::Int, a, b) = Parallelogram(Nx, Ny, Float64[a, 0], Float64[0, b])

vec_a(p::Parallelogram) = p.a
vec_b(p::Parallelogram) = p.b
Base.size(sq::Parallelogram) = (sq.Nx, sq.Ny)

function isconnected(sq::AbstractParallelogram, i::Int, j::Int)
    CI = CartesianIndices(sq)
    ci = CI[i]
    cj = CI[j]
    (abs(ci.I[1] - cj[1]) == 1 && ci.I[2] == cj.I[2]) ||
    (abs(ci.I[2] - cj[2]) == 1 && ci.I[1] == cj.I[1])
end


text:
const Location2D = SVector{2, Float64}
const LatLonCoords = NamedTuple{(:lat, :lon), Tuple{Float64,Float64}}
LatLonCoords() = (lat = 0.0, lon = 0.0)

"""
Return the vector form of a latitude-longitude point [lat, lon]
"""
function convert_to_vector(c::LatLonCoords)
    return Location2D(c.lat, c.lon)
end

vector_type(::Type{LatLonCoords}) = Location2D

@with_kw struct RouteWaypoint
    stop_id::Int64
    arrival_time::Float64
end

"""
    OffTransitGraph{LOC}

Represents the set of locations not on the transit graph.
"""
@with_kw struct OffTransitGraph{LOC}
    depots::Vector{LOC}
    sites::Vector{LOC}
end

"""
    TransitGraph{LOC}

Stores information about the transit network. The integer keys are consistent across routes and metadata.
"""
@with_kw struct TransitGraph{LOC}
    stop_to_location::Dict{Int64,LOC}               # Maps Stop ID to location
    transit_trips::Vector{Vector{RouteWaypoint}}    # Always in chronological order
    transit_capacity::Vector{Int64}                 # Capacity for each trip - matches transit_trips
end

"""
    CityParams

Contains the bounding box (latitude-longitude) parameters for a city.
"""
@with_kw struct CityParams
    lat_start::Float64
    lat_end::Float64
    lon_start::Float64
    lon_end::Float64
end

"""
    parse_city_params(param_file::String)

Return the CityParams object from parsing the corresponding TOML file.
"""
function parse_city_params(param_file::String)

    params_dict = TOML.parsefile(param_file)

    return CityParams(lat_start = params_dict["LATSTART"],
                      lat_end = params_dict["LATEND"],
                      lon_start = params_dict["LONSTART"],
                      lon_end = params_dict["LONEND"])
end

"""
    DroneParams

Contains the operation parameters of a drone (system is homogeneous)
"""
@with_kw struct DroneParams
    max_distance::Float64
    max_air_time::Float64
    avg_speed::Float64
    stop_buffer_time::Float64
end

"""
    parse_drone_params(param_file::String)

Return the DroneParams object from parsing the corresponding TOML file.
"""
function parse_drone_params(param_file::String)

    params_dict = TOML.parsefile(param_file)

    return DroneParams(max_distance = params_dict["MAX_DISTANCE"],
                       max_air_time = params_dict["MAX_AIR_TIME"],
                       avg_speed = params_dict["AVG_SPEED"],
                       stop_buffer_time = params_dict["STOP_BUFFER_TIME"])
end

"""
    MAPFTransitState{LOC}

Denotes the time-stamped location of a transit stop.
"""
@with_kw struct MAPFTransitState{LOC}
    time::Float64       # Planned time to be at the vertex. For routewaypoint, should be before ETA
    location::LOC       # The actual vertex string ID (TG or OTG); used with get_location_or_routept
end

Base.isequal(s1::MAPFTransitState, s2::MAPFTransitState) = (s1.time, s1.location) == (s2.time, s2.location)

@enum ActionType Stay=1 Fly=2 Board=3
struct MAPFTransitAction <: MAPFAction
    action::ActionType
end

@enum ConflictType Transfer=1 Capacity=2

"""
Board Conflict - Two drones go to the same route waypoint at same time
Capacity Conflict - A transit route has more drones than it can accommodate
"""
@with_kw struct MAPFTransitConflict <: MAPFConflict
    type::ConflictType
    overlap_vertices::Set{Int64}          # Which vertex IDs are in the conflicting set
    agent_to_state_idx::Dict{Int64,Int64} # Maps overlapping agent ID to index along state vect
    cap::Int64                      = 1
end

@enum ConstraintSubPath towards=1 from=2

"""
Breaks up the transit vertex constraint into vertices on the path FROM depot -> package
and the vertices TO destination, i.e., from : depot->package, to : package->delivery
"""
struct MAPFTransitConstraints <: MAPFConstraints
    avoid_vertex_map::Dict{ConstraintSubPath,Set{Int64}}
end
MAPFTransitConstraints() = MAPFTransitConstraints(Dict{ConstraintSubPath,Set{Int64}}(towards => Set{Int64}(),
                                                                                from => Set{Int64}()))

Base.isempty(mtc::MAPFTransitConstraints) = isempty(mtc.avoid_vertex_map[towards]) && isempty(mtc.avoid_vertex_map[from])

"""
Stores the unique integer ID and string representation of an operation graph vertex.
Strings are "s-i" for package site i, "d-i" for depot i, "r-i-j" for waypoint j of transit route i.
"""
@with_kw mutable struct MAPFTransitVertexState{MTS <: MAPFTransitState} <: MAPFState
    idx::Int64
    state::MTS
    vertex_str::String
end

function reset_time(vs::MAPFTransitVertexState, time::Float64=0.0)
    vs.state = MAPFTransitState(time, vs.state.location)
end

function Graphs.vertex_index(g::SimpleVListGraph{MAPFTransitVertexState}, v::MAPFTransitVertexState)
    return v.idx
end

const AgentTask = NamedTuple{(:origin, :site, :dest)}

"""
Tracks the current state of the agent; current task; whether it has delivered its package, the distance flown already.
"""
@with_kw mutable struct AgentState
    task::AgentTask
    site_crossed::Bool              = false
    next_finish_time::Float64       = 0.0   # What is the time-stamp of the next event point
    dist_flown::Float64             = 0.0
end

"""
    MAPFTransitEnv{OTG <: OffTransitGraph, TG <: TransitGraph, MVTS <: MAPFTransitVertexState} <: MAPFEnvironment

The overarching structure to track the MAPF-TN contextual information.

N.B. - There are some hacks here to track drone state with time for the replanning. The MAPF library
is typically not well-suited to receding-horizon execution. Email me if you need to figure that out.
"""
@with_kw mutable struct MAPFTransitEnv{OTG <: OffTransitGraph,
                                       TG <: TransitGraph, MVTS <: MAPFTransitVertexState} <: MAPFEnvironment
    off_transit_graph::OTG
    transit_graph::TG
    state_graph::SimpleVListGraph{MVTS}                 # Vertex IDs are d-1 etc, s-1 etc, r-1-1 etc.
    agent_states::Vector{AgentState}
    depot_sites_to_vtx::Dict{String,Int64}              # Maps depot and site IDs to their vertex ID in graph - needed for start-goal IDXs
    trip_to_vtx_range::Vector{Tuple{Int64,Int64}}       # Maps the ID from NNTree to stop id
    stop_idx_to_trips::Dict{Int64,Set{Int64}}           # Maps the stop ID to trips passing through them
    aug_trips_fws_dists::Matrix{Float64}                # The distance matrix from the augmented trip metagraph
    drone_params::DroneParams
    dist_fn::Function                                   # Implements the distance metric between coordinates (lat-long/grid point)
    curr_site_points::Vector{Int64}                     # For each agent, which point on its solution vector is the intermediate goal, i.e. p
    curr_goal_idx::Int64                   = 0          # Tracks the current goal index for each low-level single-agent search
    plan_ref_times::Vector{Float64}        = Float64[]  # Tracks the system time from which each agent is being planned for
    # Solution Diagnostics
    valid_transit_options::Vector{Int64}   = Int64[]    # The number of transit options used by each agent
    any_invalid_path::Bool                 = false      # Is either subpath of an agent an invalid one, i.e. no weight-constrained path found
    valid_path_dists::Vector{Float64}      = Float64[]  # For the valid agent subpaths, what is the effective distance covered
    num_global_conflicts::Int64            = 0          # The number of high-level ECBS conflicts
    threshold_global_conflicts::Int64      = 10         # The threshold for high-level conflicts after which an error is thrown. TODO: Need to diagnose better
end

function get_vertex_location(env::MAPFTransitEnv, vtx_str::String)

    splitstr = split(vtx_str, "-")

    if splitstr[1] == "d"
        return env.off_transit_graph.depots[parse(Int64, splitstr[2])]
    elseif splitstr[1] == "s"
        return env.off_transit_graph.sites[parse(Int64, splitstr[2])]
    else
        rwp = env.transit_graph.transit_trips[parse(Int64, splitstr[2])][parse(Int64, splitstr[3])]
        return env.stop_to_location[rwp.stop_id]
    end

end


text:
module SorbetHelperReturnBerry

using ..Ahorn, Maple

@mapdef Entity "SorbetHelper/ReturnBerry" ReturnBerry(
    x::Integer,
    y::Integer,
    winged::Bool=false,
    delay::Number=0.3,
    checkpointID::Integer=-1,
    order::Integer=-1,
)

const placements = Ahorn.PlacementDict(
    "Strawberry (With Node Based Return) (Sorbet Helper)" => Ahorn.EntityPlacement(
        ReturnBerry,
        "point",
        Dict{String, Any}(),
        function(entity::ReturnBerry)
            entity.data["nodes"] = [
                (Int(entity.data["x"]) + 32, Int(entity.data["y"])),
                (Int(entity.data["x"]) + 64, Int(entity.data["y"]))
            ]
        end
    ),
    "Strawberry (Winged, With Node Based Return) (Sorbet Helper)" => Ahorn.EntityPlacement(
        ReturnBerry,
        "point",
        Dict{String, Any}(
            "winged" => true
        ),
        function(entity::ReturnBerry)
             entity.data["nodes"] = [
                 (Int(entity.data["x"]) + 32, Int(entity.data["y"])),
                 (Int(entity.data["x"]) + 64, Int(entity.data["y"]))
             ]
        end
    )
)

const sprites = Dict{Tuple{Bool, Bool}, String}(
    (false, false) => "collectables/strawberry/normal00",
    (true, false) => "collectables/strawberry/wings01",
    (false, true) => "collectables/ghostberry/idle00",
    (true, true) => "collectables/ghostberry/wings01",
)

Ahorn.nodeLimits(entity::ReturnBerry) = 2, -1

function Ahorn.selection(entity::ReturnBerry)
    x, y = Ahorn.position(entity)

    nodes = get(entity.data, "nodes", Tuple{Int, Int}[])
    winged = get(entity.data, "winged", false)
    seeded = length(nodes) > 2

    sprite = sprites[(winged, seeded)]
    seedSprite = "collectables/strawberry/seed00"
    bubbleSprite = "characters/player/bubble"

    res = Ahorn.Rectangle[Ahorn.getSpriteRectangle(sprite, x, y)]

    for (i, node) in enumerate(nodes)
        nx, ny = node
        if i < 3
            push!(res, Ahorn.getSpriteRectangle(bubbleSprite, nx, ny))
        else
            push!(res, Ahorn.getSpriteRectangle(seedSprite, nx, ny))
        end
    end

    return res
end

function Ahorn.renderSelectedAbs(ctx::Ahorn.Cairo.CairoContext, entity::ReturnBerry)
    x, y = Ahorn.position(entity)
    px, py = x, y
    nodes = get(entity.data, "nodes", Tuple{Int, Int}[])
    winged = get(entity.data, "winged", false)
    seeded = length(nodes) > 2

    bubbleSprite = "characters/player/bubble"

    for (i, node) in enumerate(nodes)
        nx, ny = Int.(node)

        if i < 3
            Ahorn.drawArrow(ctx, px, py, nx, ny, (1.0, 1.0, 1.0, 0.5), headLength=6)
            Ahorn.drawSprite(ctx, bubbleSprite, nx, ny)
            px, py = nx, ny
        else
            Ahorn.drawLines(ctx, Tuple{Number, Number}[(x, y), (nx, ny)], Ahorn.colors.selection_selected_fc)
        end
    end
end

function Ahorn.renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::ReturnBerry, room::Maple.Room)
    x, y = Ahorn.position(entity)
    px, py = x, y
    nodes = get(entity.data, "nodes", Tuple{Int, Int}[])
    winged = get(entity.data, "winged", false)
    seeded = length(nodes) > 2

    sprite = sprites[(winged, seeded)]
    seedSprite = "collectables/strawberry/seed00"
    bubbleSprite = "characters/player/bubble"

    for (i, node) in enumerate(nodes)
        nx, ny = Int.(node)
        if i < 3
            Ahorn.drawArrow(ctx, px, py, nx, ny, (1.0, 1.0, 1.0, 0.3), headLength=6)
            Ahorn.drawSprite(ctx, bubbleSprite, nx, ny, alpha=0.35)
            px, py = nx, ny
        else
            Ahorn.drawSprite(ctx, seedSprite, nx, ny)
        end
    end

    Ahorn.drawSprite(ctx, sprite, x, y)

end

end

max_stars_repo_name: Kitty-Cats/SorbetHelper

max_issues_repo_name: Kitty-Cats/SorbetHelper

max_forks_repo_name: Kitty-Cats/SorbetHelper


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, Pkg

name = "LASzip"
version = v"3.4.3"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/LASzip/LASzip.git", "1ab671e42ff1f086e29d5b7e300a5026e7b8d69b")
]

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

# Patch to find unordered_map from tr1, logic there doesn't work for MINGW32.
if [[ "${target}" == *-mingw* ]]; then
    sed -i '/add_definitions(-DUNORDERED)/d' LASzip/src/CMakeLists.txt;
fi

mkdir LASzip/build
cd LASzip/build/
cmake -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_CXX_FLAGS="-std=c++11" -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} -DCMAKE_BUILD_TYPE=RelWithDebInfo -DCMAKE_VERBOSE_MAKEFILE=OFF ..
cmake --build . --target install --config Release
"""

# 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 = [
    LibraryProduct(["liblaszip", "liblaszip3"], :liblaszip)
]

# 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: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
"""
    run_yaml(filename::String)
    run_yaml(filenames::Array{String,1})
This file runs (a) yaml file(s) if it includes the line `Editable: true`. If given a list of
input files, the function will print a list of yaml files that were not marked as editable.
A file might not be "editable" if SLiDE functionality cannot make all of the specified edits.

# Arguments
- `filename::String` or `filenames::Array{String,1}`: yaml file name (or list of propertynames) to run

# Return
- `filename::String` or `filenames::Array{String,1}`: yaml file name (or list of propertynames) that
    was/were not ran by the function because they were annotated with `Editable: false`
"""
function run_yaml(filename::String; save::Bool=true)
    println("\n$filename")
    println("  Reading yaml file...")

    y = read_file(filename)

    if haskey(y, "Editable") && y["Editable"]
        # Iteratively make all edits and update output data file in accordance with
        # operations defined in the input YAML file.
        println("  Parsing and standardizing...")
        df = unique(edit_with(y; print_status=true))
        _save_datastream(df, y; save=save)
        return df
    else
        println("  Skipping... (not editable)")
        return filename
    end
end

function run_yaml(lst::AbstractArray; save::Bool=true)
    return Dict(_inp_key(x) => run_yaml(x; save=save) for x in lst)
end

function run_yaml(files::Array{String,1}; save::Bool=true)
    @info("Standardize data files.")
    files = Dict(_inp_key(f) => run_yaml(f) for f in files)

    # Print message indicating which, if any, files could not be executed.
    unedited = [f for f in values(files) if typeof(f) == String]
    if length(unedited) > 0
        @warn(string("run_yaml() generated no output for:", string.("\n  ", unedited)...,
            "\nAdd \"Editable: true\" to yaml file to run automatically."))
    end

    # !!!! Would be nice to track output directories created (output from _save_datastream)
    # and delete after all of the edits in a directory have been made.

    return files
end

"""
    write_yaml(path, file::XLSXInput)
This function reads an XLSX file and writes a new yaml file containing the information in
each spreadsheet column. Sheet propertynames in the XLSX file correspond to the directory where new
files will be printed (`path/file.sheet/`). Yaml files will be named after the text in the
column's first row.

# Arguments
- `path::String` or `path::Array{String,1}`: XLSX file location. New yaml files will be
    printed here, as well.
- `file::XLSXInput` or `files::Array{XLSXInput,1}`: XLSX file information (file name, sheet name, sheet range) or list of multiple sheets.

# Returns
- `filenames::Array{String,1}`: List of yaml files
"""
function write_yaml(path, file::XLSXInput; overwrite::Bool=true)
    @info("Generating yaml files from " * file.name * "\n\tSheet: " * file.sheet)

    # Aesthetics: Define a message to print at the top of the YAML file and
    # list all key words in the yaml file and use to add (purely aesthetic) spacing.
    HEADER = "# Autogenerated from " * joinpath(path, file.name) * ", sheet: " * file.sheet
    KEYS = string.([IU.subtypes.(IU.subtypes(DataStream))...; "PathIn"], ":")

    # Read the XLSX file, identify relevant (not "missing") columns,
    # and generate list of resultant yaml file propertynames.
    println("  Reading ", file.name)
    df = read_file(path, file)
    df = df[:, .!occursin.(:missing, propertynames(df))]

    # Make sure the path exists and save the names of yaml files to generaate.
    # path = joinpath(SLIDE_DIR, path, file.sheet)
    path = joinpath(path, "readfiles", file.sheet)
    !isdir(path) && mkpath(path)

    files = [joinpath(path, "$k.yml") for k in propertynames(df)]

    # Iterate through columns. For each column, create a new yaml file and fill it with
    # the column text. Mark the yaml file as "Autogenerated" and ensure one space only
    # between each input file element that corresponds with a SLiDE Datastream subtype.
    for (col, yamlfile) in zip(propertynames(df), files)

        if overwrite == false && isfile(yamlfile)
            println("  Skipping overwrite. Found: $yamlfile")
        else
            println("  Generating $yamlfile")

            lines = dropmissing(df, col)[:,col]
            lines = lines[match.(r"\S.*", lines) .!== nothing]

            open(yamlfile, "w") do f
                println(f, HEADER)
                for line in lines
                    any(occursin.(KEYS, line)) && (println(f, ""))
                    println(f, line)
                end
            end
        end
    end
    return files
end

function write_yaml(path::Array{String,1}, file::XLSXInput; overwrite::Bool=true)
    write_yaml(joinpath(path...), file; overwrite=overwrite)
end

function write_yaml(path, files::Array{XLSXInput,1}; overwrite::Bool=true)
    if length(unique(get_descriptor.(files))) < length(files)
        @error("Input XLSX files must have unique descriptors.")
    end
    # files = Dict(_inp_key(file) => file for file in files)
    # files = Dict(k => write_yaml(path, file; overwrite=overwrite) for (k,file) in files)
    return [write_yaml(path, file; overwrite=overwrite) for file in files]
end

"""
"""
function _save_datastream(df::DataFrame, y::Dict; save=true)
    save_path = save == true ? "" : save
    haskey(y, "Temporary")   && (y["FileOut"] = _temporary(y["FileOut"]))
    haskey(y, "Development") && (y["FileOut"] = _development(y["FileOut"]))

    path = joinpath(save_path, ensurearray(y["PathOut"])...)
    file = joinpath(path, y["FileOut"])
    !isdir(path) && mkpath(path)

    println("  Saving $file...")
    CSV.write(file, df)
    return path
end

_temporary(x::String) = "__" * x
_development(x::String) = "_" * x

_isdevelopment(x::String) = match(r"^_[^_]", x) !== nothing
_istemporary(x::String) = match(r"^_{2}[^_]", x) !== nothing

"""
Delete temporary files. These are pre-pended "_".
"""
function _delete_temporary(path::String)
    curr_dir = pwd()
    files = readdir(path)
    files = files[SLiDE._istemporary.(files)]
    cd(path)
    [rm(file) for file in files]
    cd(curr_dir)
    return nothing
end

max_stars_repo_name: NREL/SLiDE

max_issues_repo_name: NREL/SLiDE

max_forks_repo_name: NREL/SLiDE


text:
@cache struct MidpointSplittingCache{uType,rateType,WType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  uprev2::uType
  tmp::uType
  fsalfirst::rateType
  W::WType
  k::rateType
end

function alg_cache(alg::MidpointSplitting,u,rate_prototype,uEltypeNoUnits,
                   tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Val{true})
  W = false .* vec(rate_prototype) .* vec(rate_prototype)' # uEltype?
  k = zero(rate_prototype); fsalfirst = zero(rate_prototype)
  MidpointSplittingCache(u,uprev,uprev2,similar(u),fsalfirst,W,k)
end

struct MidpointSplittingConstantCache <: OrdinaryDiffEqConstantCache
end

function alg_cache(alg::MidpointSplitting,u,rate_prototype,uEltypeNoUnits,
                   tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Val{false})
  MidpointSplittingConstantCache()
end

struct LinearExponentialConstantCache <: OrdinaryDiffEqConstantCache end

alg_cache(alg::LinearExponential,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,
  tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Val{false}) = LinearExponentialConstantCache()

@cache struct LinearExponentialCache{uType,rateType,KsType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  tmp::uType
  rtmp::rateType
  KsCache::KsType # different depending on alg.krylov
end

function alg_cache(alg::LinearExponential,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,
  tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Val{true})
  tmp = similar(u)
  rtmp = zero(rate_prototype)
  n = length(u)
  T = eltype(u)
  m = min(alg.m, n)

  if alg.krylov == :off
    KsCache = nothing
  elseif alg.krylov == :simple
    Ks = KrylovSubspace{T}(n, m)
    expv_cache = ExpvCache{T}(m)
    KsCache = (Ks, expv_cache)
  elseif alg.krylov == :adaptive
    KsCache = _phiv_timestep_caches(u, m, 0)
  else
    throw(ArgumentError("Unknown krylov setting $(alg.krylov). Can be :off, :simple or :adaptive."))
  end
  LinearExponentialCache(u,uprev,tmp,rtmp,KsCache)
end


text:
#=

We start with the cable equation, spatially discretize it, then temporal discretize it. We have two methods of temporally discretizing the cable equation, backward Euler and Crank-Nicolson. Luckily, these are easily related, so we only need to derive the backward Euler cable equation.

The main matrix equation to go from one time step to another for an array of voltage values for all nodes, V_n

 A \ (V_n - sum(i)) = delta_V

Where A is basically the matrix of resistances between nodes, and i is the current from each channel

Only the diagonals of A change from iteration to iteration, by a factor proportional to the conductances of the active components (di/dv).

=#

function fillA!(neuron::Neuron)

    for i=1:length(neuron.nodes)

        if neuron.nodes[i].parent != 0

            r_p=neuron.nodes[neuron.nodes[i].parent].ri.r+neuron.nodes[i].ri.l
            area_p=neuron.nodes[neuron.nodes[i].parent].area.t
            area_n=neuron.nodes[i].area.t

            neuron.b[i]=100/(r_p*area_n)
            neuron.a[i]=100/(r_p*area_p)

            neuron.diag_old[i]=100/(r_p*area_n)+.001*neuron.Cm/neuron.dt

            for j in neuron.nodes[i].children
                r_c=neuron.nodes[i].ri.r+neuron.nodes[j].ri.l
                neuron.diag_old[i] += 100/(r_c*area_n)
            end
        else

            neuron.diag_old[i]=.001*neuron.Cm/neuron.dt
            area_n=neuron.nodes[i].area.t

            for j in neuron.nodes[i].children
                r_c=neuron.nodes[i].ri.r+neuron.nodes[j].ri.l
                neuron.diag_old[i] += 100/(r_c*area_n)
            end
        end

        neuron.par[i]=neuron.nodes[i].parent
    end
    nothing
end

function initialcon!(neuron::Neuron, vi=-65.0,dt=.025)

    neuron.dt=dt

    neuron.v[:]=vi

    for i=1:4
        init!(getfield(neuron,i),neuron.v,neuron.internal_nodes[i])
    end
end

function main(neuron::Neuron)

    #t=tentry+dt for euler, t=tentry+dt/2 for CN
    #ext=false

    #reset membrane current
    @inbounds @simd for i=1:length(neuron.v)
        neuron.i_vm[i] = 0.0
    end

    for ind=1:4

        cur!(getfield(neuron,ind),neuron.v,neuron.i_vm,neuron.internal_nodes[ind])
        cur!(getfield(neuron,ind),neuron.v1,neuron.i2,neuron.internal_nodes[ind])

    end

    @inbounds @simd for i=1:length(neuron.v)
        neuron.divm[i] = (neuron.i2[i]-neuron.i_vm[i])/.001
    end

    rhs_diag!(neuron)

    #solve A \ rhs to get delta_v
    hines_solve!(neuron)

    #=
    #update voltages v_new = delta_v + v_old for euler, v_new  = 2*delta_v + v_old for CN
    if ext==true
        neuron.vext[:] += neuron.delta_vext
        neuron.v[:] += neuron.delta_v + neuron.delta_vext
    else
        add_delta!(neuron)
    end
    =#
    add_delta!(neuron)

    #find non voltage states (like gate variables for conductances)
    con!(neuron.soma,neuron.v,neuron.internal_nodes[1])
    con!(neuron.axon,neuron.v,neuron.internal_nodes[2])
    con!(neuron.dendrite,neuron.v,neuron.internal_nodes[3])
    con!(neuron.apical,neuron.v,neuron.internal_nodes[4])

    nothing
end

function add_delta!(neuron::Neuron)
    @fastmath @inbounds @simd for i=1:length(neuron.v)
        neuron.v[i] += neuron.rhs[i]
        neuron.rhs[i] = 0.0
        neuron.v1[i]=neuron.v[i]+.001
        neuron.i2[i]=0.0
    end
    nothing
end

function rhs_diag!(neuron::Neuron)

    dv=0.0

    @fastmath @inbounds for i=1:(length(neuron.v)-2)
        dv=neuron.v[neuron.par[i]]-neuron.v[i]
        neuron.rhs[i] += neuron.b[i]*dv
        neuron.rhs[neuron.par[i]] -= neuron.a[i]*dv
        neuron.d[i] = neuron.diag_old[i] + neuron.divm[i]
        neuron.rhs[i] -= neuron.i_vm[i]
    end

    i=length(neuron.v)-1

    neuron.d[i] = neuron.diag_old[i] + neuron.divm[i]
    neuron.rhs[i] -= neuron.i_vm[i]

    i=length(neuron.v)

    dv=neuron.v[neuron.par[i]]-neuron.v[i]
    neuron.rhs[i] += neuron.b[i]*dv
    neuron.rhs[neuron.par[i]] -= neuron.a[i]*dv
    neuron.d[i] = neuron.diag_old[i] + neuron.divm[i]
    neuron.rhs[i] -= neuron.i_vm[i]

    nothing
end

function hines_solve!(neuron::Neuron)

    p=0.0

    @fastmath @inbounds for i=(length(neuron.d)-3):-1:1

	p = -neuron.a[i] / neuron.d[i]
	neuron.d[neuron.par[i]] -= p * -neuron.b[i]
	neuron.rhs[neuron.par[i]] -= p * neuron.rhs[i]

    end

    i=length(neuron.d)-2

    p = -neuron.a[i] / neuron.d[i]
    neuron.d[neuron.par[i]] -= p * -neuron.b[i]
    neuron.rhs[neuron.par[i]] -= p * neuron.rhs[i]

    i=length(neuron.d)-2

    p = -neuron.a[i] / neuron.d[i]
    neuron.d[neuron.par[i]] -= p * -neuron.b[i]
    neuron.rhs[neuron.par[i]] -= p * neuron.rhs[i]

    neuron.rhs[length(neuron.d)-1] /= neuron.d[length(neuron.d)-1]

    i=length(neuron.d)

    neuron.rhs[i] -= -neuron.b[i] * neuron.rhs[neuron.par[i]]
    neuron.rhs[i] /= neuron.d[i]

    @fastmath @inbounds for i=1:(length(neuron.d)-2)
	neuron.rhs[i] -= -neuron.b[i] * neuron.rhs[neuron.par[i]]
	neuron.rhs[i] /= neuron.d[i]
    end
    nothing
end

max_stars_repo_name: paulmthompson/JNeuron

max_issues_repo_name: paulmthompson/JNeuron

max_forks_repo_name: paulmthompson/JNeuron


text:
using NCDatasets
using Statistics
using Interpolations
using LinearAlgebra
using Glob
using JLD
using JSON
using Random
# EKP modules
using EnsembleKalmanProcesses.ParameterDistributionStorage

"""
    run_SCAMPy(
        u::Array{FT, 1},
        u_names::Array{String, 1},
        y_names::Union{Array{String, 1}, Array{Array{String,1},1}},
        scampy_dir::String,
        scm_data_root::String,
        scm_names::Array{String, 1},
        ti::Union{Array{FT,1}, Array{Array{FT,1},1}},
        tf::Union{Array{FT,1}, Array{Array{FT,1},1}, Nothing} = nothing;
        norm_var_list = nothing,
        P_pca_list = nothing,
    ) where {FT<:AbstractFloat}

Run call_SCAMPy.sh using a set of parameters u and return
the value of outputs defined in y_names, possibly after
normalization and projection onto lower dimensional space
using PCA.

Inputs:
 - u                :: Values of parameters to be used in simulations.
 - u_names          :: SCAMPy names for parameters `u`.
 - y_names          :: Name of outputs requested for each flow configuration.
 - scampy_dir       :: Path to SCAMPy directory
 - scm_data_root    :: Path to input data for the SCM model.
 - scm_names        :: Names of SCAMPy cases
 - ti               :: Vector of starting times for observation intervals.
                        If `tf=nothing`, snapshots at `ti` are returned.
 - tf               :: Vector of ending times for observation intervals.
 - norm_var_list    :: Pooled variance vectors. If given, use to normalize output.
 - P_pca_list       :: Vector of projection matrices `P_pca` for each flow configuration.
Outputs:
 - sim_dirs         :: Vector of simulation output directories
 - g_scm            :: Vector of model evaluations concatenated for all flow configurations.
 - g_scm_pca        :: Projection of `g_scm` onto principal subspace spanned by eigenvectors.
"""
function run_SCAMPy(
        u::Array{FT, 1},
        u_names::Array{String, 1},
        y_names::Union{Array{String, 1}, Array{Array{String,1},1}},
        scampy_dir::String,
        scm_data_root::String,
        scm_names::Array{String, 1},
        ti::Union{Array{FT,1}, Array{Array{FT,1},1}},
        tf::Union{Array{FT,1}, Array{Array{FT,1},1}, Nothing} = nothing;
        norm_var_list = nothing,
        P_pca_list = nothing,
    ) where {FT<:AbstractFloat}

    # Check parameter dimensionality
    @assert length(u_names) == length(u)

    # run SCAMPy and get simulation dirs
    sim_dirs = run_SCAMPy_handler(u, u_names, scampy_dir, scm_names, scm_data_root)

    # Check consistent time interval dims
    @assert length(ti) == length(sim_dirs)

    g_scm = zeros(0)
    g_scm_pca = zeros(0)

    if typeof(ti) == Array{FT,1} # 1 interval per simulation
        for (i, sim_dir) in enumerate(sim_dirs)
            ti_ = ti[i]
            tf_ = !isnothing(tf) ? tf[i] : nothing
            y_names_ = typeof(y_names)==Array{Array{String,1},1} ? y_names[i] : y_names

            g_scm_flow = get_profile(sim_dir, y_names_, ti = ti_, tf = tf_)
            if !isnothing(norm_var_list)
                g_scm_flow = normalize_profile(g_scm_flow, y_names_, norm_var_list[i])
            end
            append!(g_scm, g_scm_flow)
            if !isnothing(P_pca_list)
                append!(g_scm_pca, P_pca_list[i]' * g_scm_flow)
            end
        end
    elseif typeof(ti) == Array{Array{FT,1},1} # multiple intervals per simulation
        config_num = 1
        for (i, sim_dir) in enumerate(sim_dirs)
            y_names_ = typeof(y_names)==Array{Array{String,1},1} ? y_names[i] : y_names
            for (j, ti_j) in enumerate(ti[i]) # Loop on time intervals per sim
                tf_j = !isnothing(tf) ? tf[i][j] : nothing
                g_scm_flow = get_profile(sim_dir, y_names_, ti = ti_j, tf = tf_j)
                if !isnothing(norm_var_list)
                    g_scm_flow = normalize_profile(g_scm_flow, y_names_, norm_var_list[config_num])
                end
                append!(g_scm, g_scm_flow)
                if !isnothing(P_pca_list)
                    append!(g_scm_pca, P_pca_list[config_num]' * g_scm_flow)
                end
                config_num += 1
            end
        end
    end
    # penalize nan-values in output
    for i in eachindex(g_scm)
        g_scm[i] = isnan(g_scm[i]) ? 1.0e5 : g_scm[i]
    end
    if !isnothing(P_pca_list)
        for i in eachindex(g_scm_pca)
            g_scm_pca[i] = isnan(g_scm_pca[i]) ? 1.0e5 : g_scm_pca[i]
        end
        println("LENGTH OF G_SCM_ARR : ", length(g_scm))
        println("LENGTH OF G_SCM_ARR_PCA : ", length(g_scm_pca))
        return sim_dirs, g_scm, g_scm_pca
    else
        return sim_dirs, g_scm
    end
end

"""
    function run_SCAMPy_handler(
        u::Array{FT, 1},
        u_names::Array{String, 1},
        scampy_dir::String,
        scm_names::String,
        scm_data_root::String,
    ) where {FT<:AbstractFloat}

Run a list of cases using a set of parameters `u_names` with values `u`,
and return a list of directories pointing to where data is stored for
each simulation run.

Inputs:
 - u :: Values of parameters to be used in simulations.
 - u_names :: SCAMPy names for parameters `u`.
 - scampy_dir :: Path to SCAMPy directory
 - scm_names :: Names of SCAMPy cases to run
 - scm_data_root :: Path to SCAMPy case data (<scm_data_root>/Output.<scm_name>.00000)
Outputs:
 - output_dirs :: list of directories containing output data from the SCAMPy runs.
"""
function run_SCAMPy_handler(
        u::Array{FT, 1},
        u_names::Array{String, 1},
        scampy_dir::String,
        scm_names::Array{String, 1},
        scm_data_root::String,
    ) where {FT<:AbstractFloat}
    # create temporary directory to store SCAMPy data in
    tmpdir = mktempdir(pwd())

    # output directories
    output_dirs = String[]

    for simname in scm_names
        # For each scm case, fetch namelist and paramlist
        inputdir = joinpath(scm_data_root, "Output.$simname.00000")
        namelist = JSON.parsefile(joinpath(inputdir, "$simname.in"))
        paramlist = JSON.parsefile(joinpath(inputdir, "paramlist_$simname.in"))

        # update parameter values
        for (pName, pVal) in zip(u_names, u)
            paramlist["turbulence"]["EDMF_PrognosticTKE"][pName] = pVal
        end
        # write updated paramlist to `tmpdir`
        paramlist_path = joinpath(tmpdir, "paramlist_$simname.in")
        open(paramlist_path, "w") do io
            JSON.print(io, paramlist, 4)
        end

        # generate random uuid
        uuid_end = randstring(5)
        uuid_start = namelist["meta"]["uuid"][1:end-5]
        namelist["meta"]["uuid"] = "$uuid_start$uuid_end"
        # set output dir to `tmpdir`
        namelist["output"]["output_root"] = string(tmpdir, "/")
        # write updated namelist to `tmpdir`
        namelist_path = joinpath(tmpdir, "$simname.in")
        open(namelist_path, "w") do io
            JSON.print(io, namelist, 4)
        end

        # run SCAMPy with modified parameters
        main_path = joinpath(scampy_dir, "main.py")
        run(`python $main_path $namelist_path $paramlist_path`)

        push!(output_dirs, joinpath(tmpdir, "Output.$simname.$uuid_end"))
    end  # end `simnames` loop
    return output_dirs
end

"""
    obs_LES(y_names, sim_dir, ti, tf;
            z_scm = nothing, normalize = false)

Get LES output for observed variables y_names, interpolated to
z_scm (if given), and possibly normalized with respect to the pooled variance.

Inputs:
 - y_names :: Name of outputs requested for each flow configuration.
 - ti :: Vector of starting times for observation intervals. If tf=nothing,
         snapshots at ti are returned.
 - tf :: Vector of ending times for observation intervals.
 - z_scm :: If given, interpolate LES observations to given levels.
 - normalize :: If true, normalize observations and cov matrix by pooled variances.
Outputs:
 - y_ :: Mean of observations, possibly interpolated to z_scm levels.
 - y_tvar :: Observational covariance matrix, possibly pool-normalized.
"""
function obs_LES(y_names::Array{String, 1},
                    sim_dir::String,
                    ti::FT,
                    tf::FT;
                    z_scm::Union{Array{FT, 1}, Nothing} = nothing,
                    normalize = true,
                    ) where {FT<:AbstractFloat}

    y_names_les = get_les_names(y_names, sim_dir)
    y_tvar, poolvar_vec = get_time_covariance(sim_dir, y_names_les,
        ti = ti, tf = tf, z_scm=z_scm, normalize=normalize)
    y_highres = get_profile(sim_dir, y_names_les, ti = ti, tf = tf)
    if normalize
        y_highres = normalize_profile(y_highres, y_names, poolvar_vec)
    end
    if !isnothing(z_scm)
        y_ = zeros(0)
        z_les = get_profile(sim_dir, ["z_half"])
        num_outputs = Integer(length(y_highres)/length(z_les))
        for i in 1:num_outputs
            y_itp = interpolate( (z_les,),
                y_highres[1 + length(z_les)*(i-1) : i*length(z_les)],
                Gridded(Linear()) )
            append!(y_, y_itp(z_scm))
        end
    else
        y_ = y_highres
    end
    return y_, y_tvar, poolvar_vec
end

"""
    obs_PCA(y_mean, y_var, allowed_var_loss = 1.0e-1)

Perform dimensionality reduction using principal component analysis on
the variance y_var. Only eigenvectors with eigenvalues that contribute
to the leading 1-allowed_var_loss variance are retained.
Inputs:
 - y_mean :: Mean of the observations.
 - y_var :: Variance of the observations.
 - allowed_var_loss :: Maximum variance loss allowed.
Outputs:
 - y_pca :: Projection of y_mean onto principal subspace spanned by eigenvectors.
 - y_var_pca :: Projection of y_var on principal subspace.
 - P_pca :: Projection matrix onto principal subspace, with leading eigenvectors as columns.
"""
function obs_PCA(y_mean, y_var, allowed_var_loss = 1.0e-1)
    eig = eigen(y_var)
    eigvals, eigvecs = eig; # eigvecs is matrix with eigvecs as cols
    # Get index of leading eigenvalues, eigvals are ordered from low to high in julia
    # This expression recovers 1 extra eigenvalue compared to threshold
    leading_eigs = findall(<(1.0-allowed_var_loss), -cumsum(eigvals)/sum(eigvals).+1)
    P_pca = eigvecs[:, leading_eigs]
    λ_pca = eigvals[leading_eigs]
    # Check correct PCA projection
    @assert Diagonal(λ_pca) ≈ P_pca' * y_var * P_pca
    # Project mean
    y_pca = P_pca' * y_mean
    y_var_pca = Diagonal(λ_pca)
    return y_pca, y_var_pca, P_pca
end

function interp_padeops(padeops_data,
                    padeops_z,
                    padeops_t,
                    z_scm,
                    t_scm
                    )
    # Weak verification of limits for independent vars
    @assert abs(padeops_z[end] - z_scm[end])/padeops_z[end] <= 0.1
    @assert abs(padeops_z[end] - z_scm[end])/z_scm[end] <= 0.1

    # Create interpolating function
    padeops_itp = interpolate( (padeops_t, padeops_z), padeops_data,
                ( Gridded(Linear()), Gridded(Linear()) ) )
    return padeops_itp(t_scm, z_scm)
end

function padeops_m_σ2(padeops_data,
                    padeops_z,
                    padeops_t,
                    z_scm,
                    t_scm,
                    dims_ = 1)
    padeops_snapshot = interp_padeops(padeops_data,padeops_z,padeops_t,z_scm, t_scm)
    # Compute variance along axis dims_
    padeops_var = cov(padeops_data, dims=dims_)
    return padeops_snapshot, padeops_var
end

function get_profile(sim_dir::String,
                     var_name::Array{String,1};
                     ti::Float64=0.0,
                     tf=nothing,
                     getFullHeights=false)

    if length(var_name) == 1 && occursin("z_half", var_name[1])
        prof_vec = nc_fetch(sim_dir, "profiles", var_name[1])
    else
        t = nc_fetch(sim_dir, "timeseries", "t")
        dt = length(t) > 1 ? abs(t[2]-t[1]) : 0.0
        # Check that times are contained in simulation output
        ti_diff, ti_index = findmin( broadcast(abs, t.-ti) )
        if !isnothing(tf)
            tf_diff, tf_index = findmin( broadcast(abs, t.-tf) )
        end
        prof_vec = zeros(0)
        # If simulation does not contain values for ti or tf, return high value
        if ti_diff > dt
            println("ti_diff > dt ", "ti_diff = ", ti_diff, "dt = ", dt, "ti = ", ti,
                 "t[1] = ", t[1], "t[end] = ", t[end])
            for i in 1:length(var_name)
                var_ = nc_fetch(sim_dir, "profiles", "z_half")
                append!(prof_vec, 1.0e5*ones(length(var_[:])))
            end
        else
            for i in 1:length(var_name)
                if occursin("horizontal_vel", var_name[i])
                    u_ = nc_fetch(sim_dir, "profiles", "u_mean")
                    v_ = nc_fetch(sim_dir, "profiles", "v_mean")
                    var_ = sqrt.(u_.^2 + v_.^2)
                else
                    var_ = nc_fetch(sim_dir, "profiles", var_name[i])
                    # LES vertical fluxes are per volume, not mass
                    if occursin("resolved_z_flux", var_name[i])
                        rho_half=nc_fetch(sim_dir, "reference", "rho0_half")
                        var_ = var_.*rho_half
                    end
                end
                if !isnothing(tf)
                    append!(prof_vec, mean(var_[:, ti_index:tf_index], dims=2))
                else
                    append!(prof_vec, var_[:, ti_index])
                end
            end
        end
    end
    return prof_vec
end

"""
    normalize_profile(profile_vec, var_name, var_vec)

Perform normalization of profiles contained in profile_vec
using the standard deviation associated with each variable in
var_name. Variances for each variable are contained
in var_vec.
"""
function normalize_profile(profile_vec, var_name, var_vec)
    prof_vec = deepcopy(profile_vec)
    dim_variable = Integer(length(profile_vec)/length(var_name))
    for i in 1:length(var_name)
        prof_vec[dim_variable*(i-1)+1:dim_variable*i] =
            prof_vec[dim_variable*(i-1)+1:dim_variable*i] ./ sqrt(var_vec[i])
    end
    return prof_vec
end

"""
    get_time_covariance(sim_dir::String,
                     var_name::Array{String,1};
                     ti::Float64=0.0,
                     tf=0.0,
                     getFullHeights=false,
                     z_scm::Union{Array{Float64, 1}, Nothing} = nothing,
                     normalize=false)

Obtain the covariance matrix of a group of profiles, where the covariance
is obtained in time.
Inputs:
 - sim_dir :: Name of simulation directory.
 - var_name :: List of variable names to be included.
 - ti, tf :: Initial and final times defining averaging interval.
 - z_scm :: If given, interpolates covariance matrix to this locations.
 - normalize :: Boolean specifying variable normalization.
"""
function get_time_covariance(sim_dir::String,
                     var_name::Array{String,1};
                     ti::Float64=0.0,
                     tf=0.0,
                     getFullHeights=false,
                     z_scm::Union{Array{Float64, 1}, Nothing} = nothing,
                     normalize=false)

    t = nc_fetch(sim_dir, "timeseries", "t")
    # Find closest interval in data
    ti_index = argmin( broadcast(abs, t.-ti) )
    tf_index = argmin( broadcast(abs, t.-tf) )
    ts_vec = zeros(0, length(ti_index:tf_index))
    num_outputs = length(var_name)
    poolvar_vec = zeros(num_outputs)

    for i in 1:num_outputs
        var_ = nc_fetch(sim_dir, "profiles", var_name[i])
        # LES vertical fluxes are per volume, not mass
        if occursin("resolved_z_flux", var_name[i])
            rho_half=nc_fetch(sim_dir, "reference", "rho0_half")
            var_ = var_.*rho_half
        end
        # Store pooled variance
        poolvar_vec[i] = mean(var(var_[:, ti_index:tf_index], dims=2))  # vertically averaged time-variance of variable
        ts_var_i = normalize ? var_[:, ti_index:tf_index]./ sqrt(poolvar_vec[i]) : var_[:, ti_index:tf_index]
        # Interpolate in space
        if !isnothing(z_scm)
            z_les = getFullHeights ? get_profile(sim_dir, ["z"]) : get_profile(sim_dir, ["z_half"])
            # Create interpolant
            ts_var_i_itp = interpolate(
                (z_les, 1:tf_index-ti_index+1),
                ts_var_i,
                ( Gridded(Linear()), NoInterp() )
            )
            # Interpolate
            ts_var_i = ts_var_i_itp(z_scm, 1:tf_index-ti_index+1)
        end
        ts_vec = cat(ts_vec, ts_var_i, dims=1)  # dims: (Nz*num_outputs, Nt)
    end
    cov_mat = cov(ts_vec, dims=2)  # covariance, w/ samples across time dimension (t_inds).
    return cov_mat, poolvar_vec
end

function get_les_names(scm_y_names::Array{String,1}, sim_dir::String)
    y_names = deepcopy(scm_y_names)
    if "thetal_mean" in y_names
        if occursin("GABLS",sim_dir) || occursin("Soares",sim_dir)
            y_names[findall(x->x=="thetal_mean", y_names)] .= "theta_mean"
        else
            y_names[findall(x->x=="thetal_mean", y_names)] .= "thetali_mean"
        end
    end
    if "total_flux_qt" in y_names
        y_names[findall(x->x=="total_flux_qt", y_names)] .= "resolved_z_flux_qt"
    end
    if "total_flux_h" in y_names && (occursin("GABLS",sim_dir) || occursin("Soares",sim_dir))
        y_names[findall(x->x=="total_flux_h", y_names)] .= "resolved_z_flux_theta"
    elseif "total_flux_h" in y_names
        y_names[findall(x->x=="total_flux_h", y_names)] .= "resolved_z_flux_thetali"
    end
    if "u_mean" in y_names
        y_names[findall(x->x=="u_mean", y_names)] .= "u_translational_mean"
    end
    if "v_mean" in y_names
        y_names[findall(x->x=="v_mean", y_names)] .= "v_translational_mean"
    end
    if "tke_mean" in y_names
        y_names[findall(x->x=="tke_mean", y_names)] .= "tke_nd_mean"
    end
    return y_names
end

function nc_fetch(dir, nc_group, var_name)
    find_prev_to_name(x) = occursin("Output", x)
    split_dir = split(dir, ".")
    sim_name = split_dir[findall(find_prev_to_name, split_dir)[1]+1]
    ds = NCDataset(string(dir, "/stats/Stats.", sim_name, ".nc"))
    ds_group = ds.group[nc_group]
    ds_var = deepcopy( Array(ds_group[var_name]) )
    close(ds)
    return Array(ds_var)
end

"""
agg_clima_ekp(n_params::Integer, output_name::String="ekp_clima")

Aggregate all iterations of the parameter ensembles and write to file.
"""
function agg_clima_ekp(n_params::Integer, output_name::String="ekp_clima")
    # Get versions
    version_files = glob("versions_*.txt")
    # Recover parameters of last iteration
    last_up_versions = readlines(version_files[end])

    ens_all = Array{Float64, 2}[]
    for (it_num, file) in enumerate(version_files)
        versions = readlines(file)
        u = zeros(length(versions), n_params)
        for (ens_index, version_) in enumerate(versions)
            if it_num == length(version_files)
                open("../../../ClimateMachine.jl/test/Atmos/EDMF/$(version_)", "r") do io
                    u[ens_index, :] = [parse(Float64, line) for (index, line) in enumerate(eachline(io)) if index%3 == 0]
                end
            else
                open("$(version_).output/$(version_)", "r") do io
                    u[ens_index, :] = [parse(Float64, line) for (index, line) in enumerate(eachline(io)) if index%3 == 0]
                end
            end
        end
        push!(ens_all, u)
    end
    save( string(output_name,".jld"), "ekp_u", ens_all)
    return
end

"""
    precondition_ensemble!(params::Array{FT, 2}, priors,
        param_names::Vector{String}, ::Union{Array{String, 1}, Array{Array{String,1},1}},
        ti::Union{FT, Array{FT,1}}, tf::Union{FT, Array{FT,1}};
        lim::FT=1.0e3,) where {IT<:Int, FT}

Substitute all unstable parameters by stable parameters drawn from
the same prior.
"""
function precondition_ensemble!(params::Array{FT, 2}, priors,
    param_names::Vector{String}, y_names::Union{Array{String, 1}, Array{Array{String,1},1}},
    scampy_dir::String, scm_data_root::String, scm_names::Array{String, 1},
    ti::Union{FT, Array{FT,1}};
    tf::Union{FT, Array{FT,1}, Nothing}=nothing, lim::FT=1.0e4,) where {IT<:Int, FT}

    # Check dimensionality
    @assert length(param_names) == size(params, 1)
    # Wrapper around SCAMPy in original output coordinates
    g_(x::Array{Float64,1}) = run_SCAMPy(
        x, param_names, y_names, scampy_dir, scm_data_root, scm_names, ti, tf,
    )

    params_cons_i = deepcopy(transform_unconstrained_to_constrained(priors, params))
    params_cons_i = [row[:] for row in eachrow(params_cons_i')]
    N_ens = size(params_cons_i, 1)
    g_ens_arr = pmap(g_, params_cons_i) # [N_ens N_output]
    @assert size(g_ens_arr, 1) == N_ens
    N_out = size(g_ens_arr, 2)
    # If more than 1/4 of outputs are over limit lim, deemed as unstable simulation
    uns_vals_frac = sum(count.(x->x>lim, g_ens_arr), dims=2)./N_out
    unstable_point_inds = findall(x->x>0.25, uns_vals_frac)
    println(string("Unstable parameter indices: ", unstable_point_inds))
    # Recursively eliminate all unstable parameters
    if !isempty(unstable_point_inds)
        println(length(unstable_point_inds), " unstable parameters found:" )
        for j in 1:length(unstable_point_inds)
            println(params[:, unstable_point_inds[j]])
        end
        println("Sampling new parameters from prior...")
        new_params = construct_initial_ensemble(priors, length(unstable_point_inds))
        precondition_ensemble!(new_params, priors, param_names,
            y_names, scampy_dir, scm_data_root, scm_names, ti, tf=tf, lim=lim)
        params[:, unstable_point_inds] = new_params
    end
    println("\nPreconditioning finished.")
    return
end

"""
    compute_errors(g_arr, y)

Computes the L2-norm error of each elmt of g_arr
wrt vector y.
"""
function compute_errors(g_arr, y)
    diffs = [g - y for g in g_arr]
    errors = map(x->dot(x,x), diffs)
    return errors
end

"""
    cov_from_cov_list(cov_list::Array{Array{FT,2},1}; indices=nothing)

Returns a block-diagonal covariance matrix constructed from covariances
within cov_list given by the indices. If isempty(indices), use all
covariances to construct block-diagonal matrix.
"""
function cov_from_cov_list(cov_list::Array{Array{FT,2},1};
         indices=[]) where {FT<:AbstractFloat}
    size_ = isempty(indices) ? sum([length(cov[1,:]) for cov in cov_list]) :
        sum([length(cov[1,:]) for (i, cov) in enumerate(cov_list) if i in indices])

    cov_ = zeros(size_, size_)
    vars_num = 1
    for (index, small_cov) in enumerate(cov_list)
        if index in indices
            vars = length(small_cov[1,:])
            cov_[vars_num:vars_num+vars-1, vars_num:vars_num+vars-1] = small_cov
            vars_num = vars_num+vars
        end
    end
    return cov_
end

"""
    vec_from_vec_list(vec_list::Array{Array{FT,1},1}; indices=[], return_mapping=false)

Returns a vector constructed from vectors within vec_list given by the
indices. If isempty(indices), use all vectors to construct returned vector.
If return_mapping, function returns the positions of all the elements used
to construct the returned vector.
"""
function vec_from_vec_list(vec_list::Array{Array{FT,1},1};
         indices=[], return_mapping=false) where {FT<:AbstractFloat}
    vector_ = zeros(0)
    elmt_num = []
    chosen_elmt_num = []
    for (index, small_vec) in enumerate(vec_list)
        index < 2 ? append!(elmt_num, 1:length(small_vec)) :
                    append!(elmt_num, elmt_num[end]+1:elmt_num[end]+length(small_vec))
        if index in indices
            append!(vector_, small_vec)
            append!(chosen_elmt_num, elmt_num[end] - length(small_vec) + 1 : elmt_num[end])
        end
    end
    if return_mapping
        return vector_, chosen_elmt_num
    else
        return vector_
    end
end

"""
    logmean_and_logstd(μ, σ)

Returns the lognormal parameters μ and σ from the mean μ and std σ of the
lognormal distribution.
"""
function logmean_and_logstd(μ, σ)
    σ_log = sqrt(log(1.0 + σ^2/μ^2))
    μ_log = log(μ / (sqrt(1.0 + σ^2/μ^2)))
    return μ_log, σ_log
end

"""
    logmean_and_logstd_from_mode_std(μ, σ)

Returns the lognormal parameters μ and σ from the mode and the std σ of the
lognormal distribution.
"""
function logmean_and_logstd_from_mode_std(mode, σ)
    σ_log = sqrt( log(mode)*(log(σ^2)-1.0)/(2.0-log(σ^2)*3.0/2.0) )
    μ_log = log(mode) * (1.0-log(σ^2)/2.0)/(2.0-log(σ^2)*3.0/2.0)
    return μ_log, σ_log
end

"""
    mean_and_std_from_ln(μ, σ)

Returns the mean and variance of the lognormal distribution
from the lognormal parameters μ and σ.
"""
function mean_and_std_from_ln(μ_log, σ_log)
    μ = exp(μ_log + σ_log^2/2)
    σ = sqrt( (exp(σ_log^2) - 1)* exp(2*μ_log + σ_log^2) )
    return μ, σ
end

log_transform(a::AbstractArray) = log.(a)
exp_transform(a::AbstractArray) = exp.(a)


text:
module Backtest

export xbacktest

include("xbacktest.jl")

end


text:
import GLMakie.GLAbstraction: Texture, render, std_renderobject, LazyShader
import GLMakie.GLAbstraction: @frag_str, @vert_str
import ModernGL: GL_TEXTURE_2D

function FrameBuffer(context::Context, t::Texture{T,2}) where {T}
    frame_buffer = ContextCreateFramebufferFromGLTexture2D(context, GL_TEXTURE_2D, 0, t.id)
    x = FrameBuffer(frame_buffer)
    finalizer(release, x)
    return x
end

"""
Creates an interactive context with the help of glvisualize to view the texture
"""
function interactive_context(glwindow)
    context = Context(RPR_CREATION_FLAGS_ENABLE_GPU0 | RPR_CREATION_FLAGS_ENABLE_GL_INTEROP)
    w, h = widths(glwindow)
    texture = Texture(RGBA{Float16}, (w, h))
    view_tex(texture, glwindow)
    g_frame_buffer = FrameBuffer(context, texture)
    set!(context, RPR_AOV_COLOR, g_frame_buffer)
    return context, g_frame_buffer
end

"""
blocking renderloop that uses tiling
"""
function tiledrenderloop(glwindow, context, framebuffer)
    ti = TileIterator(widths(glwindow), (256, 256))
    tile_state = iterate(ti)
    while isopen(glwindow)
        glBindTexture(GL_TEXTURE_2D, 0)

        if isnothing(tile_state)
            # reset iterator
            tile_state = iterate(ti)
        end

        render(context, tile_state[1])

        tile_state = iterate(ti, tile_state[2])
        GLWindow.render_frame(glwindow)
    end
end

const tex_frag = frag"""
{{GLSL_VERSION}}
in vec2 o_uv;
out vec4 fragment_color;
out uvec2 fragment_groupid;

uniform sampler2D image;

void main(){
    vec4 color = texture(image, o_uv);
    fragment_color = color/color.a;
    fragment_groupid = uvec2(0);
}
"""

const tex_vert = vert"""
{{GLSL_VERSION}}

in vec2 vertices;
in vec2 texturecoordinates;

out vec2       o_uv;

void main(){
    o_uv        = texturecoordinates;
    gl_Position = vec4(vertices, 0, 1);
}
"""

"""
A matrix of colors is interpreted as an image
"""
function view_tex(tex::Texture{T,2}, window) where {T}
    data = Dict{Symbol,Any}(:image => tex,
                            :primitive => GLUVMesh2D(SimpleRectangle{Float32}(-1.0f0, -1.0f0, 2.0f0, 2.0f0)))
    robj = std_renderobject(data, LazyShader(tex_vert, tex_frag))
    return view(robj, window; camera=:nothing)
end


text:
using BlackBoxOptim, Gadfly

# run Borg MOAE
schaffer1(x) = (sumabs2(x), sumabs2(x .- 2.0))
res = bboptimize(schaffer1; Method=:borg_moea,
                FitnessScheme=ParetoFitnessScheme{2}(is_minimizing=true),
                SearchRange=(-10.0, 10.0), NumDimensions=3, ϵ=0.1,
                MaxSteps=15000, TraceInterval=1.0, TraceMode=:verbose);

# the parameterized version of the exact Pareto frontier
# N is the number of the problem (not fitness) dimensions
pareto_curve_func(t, ::Type{Val{N}}) where {N} = (N*t[1]^2, N*(2-t[1])^2)
pareto_curve = BlackBoxOptim.Hypersurface(pareto_curve_func,
                                          symmetric_search_space(1, (0.0, 2.0)))

# generate the set of ϵ-indexed points on the exact Pareto frontier
pareto_pts = BlackBoxOptim.generate(pareto_curve,
                                   fitness_scheme(res), Val{numdims(res)});
# calculate the distance between the solution and the exact frontier
BlackBoxOptim.IGD(pareto_curve, pareto_frontier(res), fitness_scheme(res), Val{numdims(res)})

# draw the results
plot(layer(x = [x.orig[1] for x in values(pareto_pts)], # Draw the exact Pareto frontier
           y = [x.orig[2] for x in values(pareto_pts)],
           Geom.line),
    layer(x= [fitness(x).orig[1] for x in pareto_frontier(res)], # the result of the method
          y= [fitness(x).orig[2] for x in pareto_frontier(res)],
          color= [haskey(pareto_pts, fitness(x).index) for x in pareto_frontier(res)],
          Geom.point))


text:
export Wprint, W2Mstr, W2Tex
export HasGraphicsHead, HasRecursiveGraphicsHead

function Wprint(WExpr)
    weval(W"Print"(WExpr))
    return
end

###Adda function that prints the form back to Mathematica Input-language

W2Mstr(x::Number) = "$x"
W2Mstr(z::Complex) = W2Mstr(WComplex(z))
W2Mstr(x::Rational) = "($(x.num)/$(x.den))"
function W2Mstr(x::Array)
    Dim = length(size(x))
    Str="{"
    for j in 1:size(x)[1]
        if j>1
            Str*=","
        end
        if Dim == 1
            Str*=W2Mstr(x[j])
        elseif Dim == 2
            Str*=W2Mstr(x[j,:])
        else
            ###Arrays with larger dimension: $Dim != 1,2 not implemented yet
        end
    end
    Str*="}"
    return Str
end

W2Mstr(x::MathLink.WSymbol) = x.name
function W2Mstr(x::MathLink.WExpr)
    #println("W2Mstr::",x.head.name)
    if x.head.name == "Plus"
        Str = W2Mstr_PLUS(x.args)
    elseif x.head.name == "Times"
        Str = W2Mstr_TIMES(x.args)
    elseif x.head.name == "Power"
        Str = W2Mstr_POWER(x.args)
    elseif x.head.name == "Complex"
        Str = W2Mstr_COMPLEX(x.args)
    else
        Str=x.head.name*"["
        for j in 1:length(x.args)
            if j>1
                Str*=","
            end
            Str*=W2Mstr(x.args[j])
        end
        Str*="]"
    end
    return Str
end

function W2Mstr_PLUS(x::Union{Array,Tuple})
    #println("W2Mstr_PLUS:",x)
    Str="("
    for j in 1:length(x)
        if j>1
            Str*=" + "
        end
        Str*=W2Mstr(x[j])
    end
    Str*=")"
end

function W2Mstr_TIMES(x::Union{Array,Tuple})
    #println("W2Mstr_TIMES:",x)
    Str="("
    for j in 1:length(x)
        if j>1
            Str*="*"
        end
        Str*=W2Mstr(x[j])
    end
    Str*=")"
end

function W2Mstr_POWER(x::Union{Array,Tuple})
    #println("W2Mstr_POWER:",x)
    if length(x) != 2
        error("Power takes two arguments")
    end
    Str="("*W2Mstr(x[1])*"^"*W2Mstr(x[2])*")"
end

function W2Mstr_COMPLEX(x::Union{Tuple,Array})
    #println("W2Mstr_COMPLEX:",x)
    if length(x) != 2
        error("Complex takes two arguments")
    end
    if x[1] == 0
        ###Imaginary
        Str="("*W2Mstr(x[2])*"*I)"
    elseif x[2] == 0
        ### Real
        ###Complex
        Str=W2Mstr(x[1])
    else
        ###Complex
        Str="("*W2Mstr(x[1])*"+"*W2Mstr(x[2])*"*I)"
    end
end

#### Code to produce LaTex strings
W2Tex(x::MLTypes) = weval(W`ToString@TeXForm[#]&`(x))

#### Allow latex string to be shown when supported. Relevant for the jupyter notebook.

HasRecursiveGraphicsHead(w::MathLink.WSymbol) = false

function HasRecursiveGraphicsHead(w::MathLink.WExpr)
    if HasGraphicsHead(w)
        return true
    end
    for arg in w.args
        if typeof(arg) == MathLink.WExpr
            ##Only check for MathLink Expressions
            if HasRecursiveGraphicsHead(arg)
                return true
            end
        end
    end
    return false
end

function HeadsEndsWith(HeadString,Target)
    ###Check if name ends with $Target
    if length(HeadString) >= length(Target)
        return  HeadString[end-(length(Target)-1):end] == Target
    end
    return false
end

HasGraphicsHead(w::MathLink.WSymbol) = false

function HasGraphicsHead(w::MathLink.WExpr)
    HeadString = w.head.name
    GraphicsHeadsList = ["Graphics","GeoGraphics","Graphics3D","Legended"]
    ###Check for names not based on ending on Plot or Plot3D
    if HeadString in GraphicsHeadsList
        return  true
    end

    HeadsEndsWith(HeadString,"Plot") && return true
    HeadsEndsWith(HeadString,"Chart") && return true
    HeadsEndsWith(HeadString,"Plot3D") && return true
    HeadsEndsWith(HeadString,"Chart3D") && return true
    return false
end

import Base.show
Base.show(io,::MIME"text/latex",x::MathLink.WSymbol) = print(io,"\$"*W2Tex(x)*"\$")

import Base.Multimedia.showable
function Base.Multimedia.showable(::MIME"text/latex", w::MathLink.WExpr)
    if HasRecursiveGraphicsHead(w)
        return false
    else
        return true
    end
end
function Base.show(io,::MIME"text/latex",x::MathLink.WExpr)
    print(io,"\$"*W2Tex(x)*"\$")
end

#Base.show(io,::MIME"text/latex",x::MLTypes) = print(io,"\$"*W2Tex(x)*"\$")

max_stars_repo_name: fremling/MathLinkExtras

max_issues_repo_name: fremling/MathLinkExtras

max_forks_repo_name: fremling/MathLinkExtras


text:
import Base.BaseDocs: @kw_str

"""
**欢迎来到Julia $(string(VERSION)).** 完整的手册（英文）可以在这里找到

    https://docs.julialang.org/

同样下面的网址也列出了很多的入门教程和课程

    https://julialang.org/learning/

更多中文资料和教程，也请关注Julia中文社区

    https://juliacn.github.io (境内域名 juliacn.com 还在备案中)

输入 `?`， 然后输入你想要查看帮助文档的函数或者宏名称就可以查看它们的文档。例如`?cos`, 或者 `?@time`
然后按回车键即可。

退出REPL之后，重新进入将恢复英文模式。
"""
kw"help", kw"?", kw"julia", kw""

"""
    using

`using Foo` 将会加载一个名为 `Foo` 的模块（module）或者一个包，然后其[`export`](@ref)的名称将可以直接使用。
不论是否被`export`，名称都可以通过点来访问（例如，输入`Foo.foo`来访问到`foo`）。查看[手册中关于模块的部分](@ref modules)以获取更多细节。
"""
kw"using"

"""
    import

`import Foo` 将会加载一个名为 `Foo` 的模块（module）或者一个包。
`Foo`模块中的名称可以通过点来访问到（例如，输入`Foo.foo`可以获取到`foo`）。
查看[手册中关于模块的部分](@ref modules)以获取更多细节。
"""
kw"import"

"""
    export

`export`被用来在模块中告诉Julia哪些函数或者名字可以由用户使用。例如`export foo`将在[`using`](@ref)这个module的时候使得
`foo`可以直接被访问到。查看[手册中关于模块的部分](@ref modules)以获取更多细节。
"""
kw"export"

"""
    abstract type

`abstract type`声明来一个不能实例化的类型，它将仅仅作为类型图中的一个节点存在，从而能够描述一系列相互关联的具体类型（concrete type）：
这些具体类型都是抽象类型的子节点。抽象类型在概念上使得Julia的类型系统不仅仅是一系列对象的集合。例如：

```julia
abstract type Number end
abstract type Real <: Number end
```

[`Number`](@ref)没有父节点（父类型）, 而 [`Real`](@ref) 是 `Number` 的一个抽象子类型。
"""
kw"abstract type"

"""
    module

`module` 会声明一个 `Module` 类型的实例用于描述一个独立的变量名空间。在一个模块（module）里，你可以控制
来自于其它模块的名字是否可见（通过载入，import），你也可以决定你的名字有哪些是可以公开的（通过暴露，export）。
模块使得你在在创建上层定义时无需担心命名冲突。查看[手册中关于模块的部分](@ref modules)以获取更多细节。

# 例子
```julia
module Foo
import Base.show
export MyType, foo

struct MyType
    x
end

bar(x) = 2x
foo(a::MyType) = bar(a.x) + 1
show(io::IO, a::MyType) = print(io, "MyType \$(a.x)")
end
```
"""
kw"module"

"""
    baremodule

`baremodule` 将声明一个不包含`using Base`或者`eval`定义的模块。但是它将仍然载入`Core`模块。
"""
kw"baremodule"

"""
    primitive type

`primitive type`声明了一个其数据仅仅由一系列二进制数表示的具体类型。比较常见的例子是整数类型和浮点类型。下面是一些内置
的原始类型（primitive type）：

```julia
primitive type Char 32 end
primitive type Bool <: Integer 8 end
```

名称后面的数字表达了这个类型存储所需的比特数目。目前这个数字要求是8 bit的倍数。[`Bool`](@ref)类型的声明展示了一个原始类型如何
选择成为另一个类型的子类型。
"""
kw"primitive type"

"""
    macro

`macro`定义了一种会将生成的代码包含在最终程序体中的方法，这称之为宏。一个宏将一系列输入映射到一个表达式，然后所返回的表达式将会被
直接进行编译而不需要在运行时调用`eval`函数。宏的输入可以包括表达式，字面量，符号。例如：

# 例子

```jldoctest
julia> macro sayhello(name)
           return :( println("Hello, ", \$name, "!") )
       end
@sayhello (macro with 1 method)

julia> @sayhello "小明"
Hello, 小明!
```
"""
kw"macro"

"""
    local

`local`将会定义一个新的局部变量。

查看[手册：变量作用域](@ref scope-of-variables)以获取更详细的信息。

# 例子
```jldoctest
julia> function foo(n)
           x = 0
           for i = 1:n
               local x # introduce a loop-local x
               x = i
           end
           x
       end
foo (generic function with 1 method)

julia> foo(10)
0
```
"""
kw"local"

"""
    global

`global x`将会使得当前作用域和当前作用所包含的作用域里的`x`指向名为`x`的全局变量。
查看[手册：变量作用域](@ref scope-of-variables)以获取更多信息。

# 例子
```jldoctest
julia> z = 3
3

julia> function foo()
           global z = 6 # use the z variable defined outside foo
       end
foo (generic function with 1 method)

julia> foo()
6

julia> z
6
```
"""
kw"global"

"""
    let

`let`会在每次被运行时声明一个新的变量绑定。这个新的变量绑定将拥有一个新的地址。这里的不同只有当
变量通过闭包生存在它们的作用域外时才会显现。`let`语法接受逗号分割的一系列赋值语句和变量名：

```julia
let var1 = value1, var2, var3 = value3
    code
end
```

这些赋值语句是按照顺序求值的，等号右边的表达式将会首先求值，然后才绑定给左边的变量。因此这使得 `let x = x`
这样的表达式有意义，因为这两个`x`变量将具有不同的地址。
"""
kw"let"

"""
    quote

`quote` 会将其包含的代码扩变成一个多重的表达式对象，而无需显示调用`Expr`的构造器。这称之为引用，比如说

```julia
ex = quote
    x = 1
    y = 2
    x + y
end
```
和其它引用方式不同的是，`:( ... )`形式的引用（被包含时）将会在表达式树里引入一个在操作表达式树时必须要考虑的`QuoteNode`元素。
而在其它场景下，`:( ... )`和 `quote .. end` 代码块是被同等对待的。
"""
kw"quote"

"""
    '

厄米算符（共轭转置），参见[`adjoint`](@ref)

# 例子
```jldoctest
julia> A = [1.0 -2.0im; 4.0im 2.0]
2×2 Array{Complex{Float64},2}:
 1.0+0.0im  -0.0-2.0im
 0.0+4.0im   2.0+0.0im

julia> A'
2×2 Array{Complex{Float64},2}:
  1.0-0.0im  0.0-4.0im
 -0.0+2.0im  2.0-0.0im
```
"""
kw"'"

"""
    .'

转置算符，参见[`transpose`](@ref)

# 例子
```jldoctest
julia> A = [1.0 -2.0im; 4.0im 2.0]
2×2 Array{Complex{Float64},2}:
 1.0+0.0im  -0.0-2.0im
 0.0+4.0im   2.0+0.0im

julia> A.'
2×2 Array{Complex{Float64},2}:
  1.0+0.0im  0.0+4.0im
 -0.0-2.0im  2.0+0.0im
```
"""
kw".'"

"""
    const

`const`被用来声明常数全局变量。在大部分（尤其是性能敏感的代码）全局变量应当被声明为常数。

```julia
const x = 5
```

可以使用单个`const`声明多个常数变量。
```julia
const y, z = 7, 11
```

注意`const`只会作用于一个`=`操作，因此 `const x = y = 1` 声明了 `x` 是常数，而 `y` 不是。在另一方面，
`const x = const y = 1`声明了`x`和`y`都是常数。

注意 “常数性质” 并不会强制容器内部变成常数，所以如果`x`是一个数组或者字典（举例来讲）你仍然可以给它们添加
或者删除元素。

严格来讲，你甚至可以重新定义 `const` （常数）变量，尽管这将会让编译器产生一个警告。唯一严格的要求是这个变量的
**类型**不能改变，这也是为什么常数变量会比一般的全局变量更快的原因。
"""
kw"const"

"""
    function

函数由`function`关键词定义：

```julia
function add(a, b)
    return a + b
end
```

或者是更短的形式：

```julia
add(a, b) = a + b
```

[`return`](@ref)关键词的使用方法和其它语言完全一样，但是常常是不使用的。一个没有显示声明`return`的函数将返回函数体最后一个表达式。
"""
kw"function"


text:
# module for defining world coordinate system

using Plots
pyplot()

mutable struct World
    x_min
    x_max
    y_min
    y_max

    # init
    function World(x_min::Float64, x_max::Float64,
                   y_min::Float64, y_max::Float64)
        self = new()
        self.x_min = x_min
        self.x_max = x_max
        self.y_min = y_min
        self.y_max = y_max
        return self
    end
end

function draw(self::World)
    plot([], [], aspect_ratio=true,
         xlabel="X", ylabel="Y",
         xlims=(self.x_min, self.x_max),
         ylims=(self.y_min, self.y_max),
         legend=false)
end

max_stars_repo_name: ShisatoYano/JuliaAutonomy

max_issues_repo_name: ShisatoYano/JuliaAutonomy

max_forks_repo_name: ShisatoYano/JuliaAutonomy


text: # TODO: Move SimpleLogger in here

max_stars_repo_name: YingboMa/julia

max_issues_repo_name: YingboMa/julia

max_forks_repo_name: YingboMa/julia


text:
function check_size(x, t)
    xs, ts = size(x, 2), length(t)
    if xs != ts
        throw(DimensionMismatch("first dimension of the input data is $xs, and teaching data is $ts. these dimensions must match."))
    end
end

expand(x::AbstractArray) = vcat(x, ones(1, size(x, 2)))


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Pre-defined formats for LaTeX tables.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

export tf_latex_default, tf_latex_simple

const tf_latex_default = LatexTableFormat()

const tf_latex_simple = LatexTableFormat(
    top_line       = "\\hline",
    header_line    = "\\hline",
    mid_line       = "\\hline",
    bottom_line    = "\\hline",
    left_vline     = "|",
    mid_vline      = "|",
    right_vline    = "|",
    header_envs    = ["textbf"],
    subheader_envs = ["texttt"]
   )


text:
using Gen
using Distributions
using LinearAlgebra
using OffsetArrays
include("structs.jl")
include("distributions.jl")

"""
`confusion_matrix` must be a 0:num_object_categories x 0:num_object_categories array
where V_{ij} = P(j observed | object i). The 0th object category
represents a hallucination. The column V_{:, 0} is 0
and the row V_{0, :} is a probability distribution with
support 1, ..., num_object categories.
"""
@gen function gen_scene_and_observations(
    max_objects::Int64,
    num_object_categories::Int64,
    num_observations::Int64,
    room::RoomParams,
    camera::CameraParams,
    detection_sd::Float64,
    confusion_matrix::OffsetArray{Float64, 2, Array{Float64, 2}},
    p_hallucination:: Float64
    )

    # Generates scene
    num_objects = @trace(uniform_discrete(1, max_objects), :num_objects)
    println("num_objects: ", num_objects)
    categories = [@trace(uniform_discrete(1, num_object_categories),
        :objects=>i=>:category) for i in 1:num_objects]
    locations = [@trace(mvuniform(room.bottom_left_lim, room.top_right_lim),
        :objects=>i=>:location) for i in 1:num_objects]
    # Generates observations
    observations = []
    object_dist = OffsetArray([k==0 ? p_hallucination : (1-p_hallucination)/num_objects
        for k in 0:num_objects], 0:num_objects)
    for i in 1:num_observations
        camera_position = @trace(mvuniform(room.bottom_left_lim, room.top_right_lim),
            :observations=>i=>:camera_location)

        object = @trace(zero_categorical(object_dist), :observations=>i=>:object)

        object_category = object == 0 ? 0 : categories[object]

        category_dist = confusion_matrix[object_category, :]

        category = @trace(zero_categorical(category_dist), :observations=>i=>:category)

        if object == 0
            # using this mvuniform as a placeholder until I switch to spherical coordinates
            # note this placeholder would break inference.
            direction = @trace(mvuniform([-1.0, -1.0, -1.0], [1.0, 1.0, 1.0]), :observations=>i=>:direction)
        else
            location = locations[object]
            direction = @trace(choose_direction(camera_position, location, detection_sd), :observations=>i=>:direction)
        end

        append!(observations, Dict(:camera_location => camera_position,
            :category => category, :direction => direction,
            :object => object))
    end
    return num_objects, categories, locations, observations
end

confusion_matrix = OffsetArray([0 0.5 0.5; 0 0.9 0.1;0 0.1 0.9], 0:2, 0:2)
args = (5, 2, 20, RoomParams(), CameraParams(), 0.5, confusion_matrix, 0.2)

(trace, _) = Gen.generate(gen_scene_and_observations, args)

Gen.get_choices(trace)

max_stars_repo_name: zhangir-azerbayev/MetaGen-New

max_issues_repo_name: zhangir-azerbayev/MetaGen-New

max_forks_repo_name: zhangir-azerbayev/MetaGen-New


text:
using GLMakie
using Makie.Colors

# Change these
Δθ = 1//16 # Must be rational. Multiplied by π later to avoid floating point errors
max_prime = 500 # List will choose all numbers less than or equal to this
framerate = 60
speed_multiplier = 50//1 # Must be Rational. Multiplier to change dot speed
total_time = 180 # seconds
set_theme!(theme_black())
resolution = (1920, 1080)
tail_length = 1000
show_dots = false

# Don't change these
primes = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97,
            101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199,
            211, 223, 227, 229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293,
            307, 311, 313, 317, 331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 389, 397,
            401, 409, 419, 421, 431, 433, 439, 443, 449, 457, 461, 463, 467, 479, 487, 491, 499,
            503, 509, 521, 523, 541, 547, 557, 563, 569, 571, 577, 587, 593, 599,
            601, 607, 613, 617, 619, 631, 641, 643, 647, 653, 659, 661, 673, 677, 683, 691,
            701, 709, 719, 727, 733, 739, 743, 751, 757, 761, 769, 773, 787, 797,
            809, 811, 821, 823, 827, 829, 839, 853, 857, 859, 863, 877, 881, 883, 887,
            907, 911, 919, 929, 937, 941, 947, 953, 967, 971, 977, 983, 991, 997]
primes_list = collect(Iterators.takewhile(<=(max_prime), primes))
y_height = ceil(1.05 * primes_list[end])
x_width = ceil(y_height * resolution[1] // resolution[2])
speed_multiplier *= (framerate//60)
Δθ /= denominator(speed_multiplier)

# Function to get next frame
function step_data!(points, current_step, cycle, top_half, Δθ, scale)
    current_step += 1
    progress = current_step * Δθ
    if progress >= 1
        push!(points, (2*scale*(cycle+1), 0))
        current_step -= 1//Δθ
        top_half = !top_half
        cycle += 1
    end

    if current_step != 0
        if top_half
            θ = π*(1 - progress)
        else
            θ = π*(1 + progress)
        end
        push!(points, (2*scale*cycle + scale*(cos(θ)+1), scale*sin(θ)))
    end

    return current_step, cycle, top_half
end

function filter_points(list, tail_length)

    num_to_remove = 0
    if length(list) > tail_length
        num_to_remove = length(list) - tail_length
    end

    counter = 1
    x_lim = x_center[] - (3*x_width/4)
    while list[counter][1] < x_lim
        num_to_remove += 1
        counter +=1
    end

    if num_to_remove > 0
        list = list[1+num_to_remove:end]
    end

    return num_to_remove
end

#######################################################################
#   ANIMATION
#######################################################################
begin
    x_center = Node(0//1)
    tick_step = ceil(x_width / 4.4)

    # Make scene
    fig = Figure(resolution = resolution)
    ax = fig[1,1] = Axis(fig,
                        xticks = MultiplesTicks(4, tick_step, string(" * ", tick_step)),
                        yticks = LinearTicks(4))

    ax.autolimitaspect = 1

    # Add scroll
    on(x_center) do x_center
        # Screen limits
        screen_begin = x_center[] - 2*x_width/3
        limits!(ax, screen_begin, screen_begin + x_width, -y_height/2, y_height/2)
    end
    notify(x_center) # Manually do first update

    # Add time
    timestamps = range(0, total_time, step=1/framerate)

    # Initialize lists
    points, current_step, cycle, top_half = Node[], [], [], []
    for p in primes_list
        push!(points, Point2[(0.0,0.0)])
        append!(current_step, [0//1])
        append!(cycle, [0])
        append!(top_half, [true])
    end

    # Add lines
    for i in 1:length(points)
        lines!(ax, points[i])
    end
    scat_points = Node([a[end] for a in [b[] for b in points]])
    if show_dots
        scatter!(scat_points)
    end

    # Record it
    record(fig, "append_animation.mp4", timestamps;
            framerate = framerate) do frame
        for i in 1:length(points)
            for j in 1:numerator(speed_multiplier)
                current_step[i], cycle[i], top_half[i] = step_data!(points[i][], current_step[i], cycle[i], top_half[i], Δθ//primes_list[i], primes_list[i]/2)
            end

            num_to_remove = filter_points(points[i][], tail_length)
            if num_to_remove > 0
                deleteat!(points[i][], 1:num_to_remove)
            end
            notify(points[i]) # Force notify update
        end

        # Draw dots
        if show_dots
            scat_points[] = [a[end] for a in [b[] for b in points]]
        end
        x_center[] += Δθ * numerator(speed_multiplier)
    end
end