text:
IS.@scoped_enum(
    InputCategory,
    BRANCH = 1,
    BUS = 2,
    DC_BRANCH = 3,
    GENERATOR = 4,
    LOAD = 5,
    RESERVE = 6,
    SIMULATION_OBJECTS = 7,
    STORAGE = 8,
)

const ENUMS = (
    AngleUnits,
    BusTypes,
    GeneratorCostModels,
    InputCategory,
    LoadModels,
    PrimeMovers,
    StateTypes,
    ThermalFuels,
    UnitSystem,
)

const ENUM_MAPPINGS = Dict()

for enum in ENUMS
    ENUM_MAPPINGS[enum] = Dict()
    for value in instances(enum)
        ENUM_MAPPINGS[enum][lowercase(string(value))] = value
    end
end

"""Get the enum value for the string. Case insensitive."""
function get_enum_value(enum, value::String)
    if !haskey(ENUM_MAPPINGS, enum)
        throw(ArgumentError("enum=$enum is not valid"))
    end

    val = lowercase(value)
    if !haskey(ENUM_MAPPINGS[enum], val)
        throw(ArgumentError("enum=$enum does not have value=$val"))
    end

    return ENUM_MAPPINGS[enum][val]
end

Base.convert(::Type{AngleUnits}, val::String) = get_enum_value(AngleUnits, val)
Base.convert(::Type{BusTypes}, val::String) = get_enum_value(BusTypes, val)
Base.convert(::Type{GeneratorCostModels}, val::String) =
    get_enum_value(GeneratorCostModels, val)
Base.convert(::Type{LoadModels}, val::String) = get_enum_value(LoadModels, val)
Base.convert(::Type{PrimeMovers}, val::String) = get_enum_value(PrimeMovers, val)
Base.convert(::Type{StateTypes}, val::String) = get_enum_value(StateTypes, val)
Base.convert(::Type{ThermalFuels}, val::String) = get_enum_value(ThermalFuels, val)


text:
using PerceptualHashes
using Base.Test

if VERSION > v"0.4.5"
  using Combinatorics
end

if VERSION >= v"0.7"
  using FFTW
end

arc1n = perceptual_hash("images/architecture1.bmp")
arc1c = perceptual_hash("images/architecture1-compressed.jpg")
arc2n = perceptual_hash("images/architecture_2.bmp")
arc2c = perceptual_hash("images/architecture_2-compressed.jpg")

@test hdist(arc1n, arc1c) == 4
@test hdist(arc2n, arc2c) == 0
@test hdist(arc1n, arc2n) == 130

dlist = readdir("images")
@test length(dlist) == 16

compressed  = map(x -> "images/"*x, filter(x -> contains(x, "compressed.jpg"), dlist))
normal      = map(x -> "images/"*x, filter(x -> contains(x, ".bmp"), dlist))

bmv_compressed  = map(perceptual_hash, compressed)
bmv_normal      = map(perceptual_hash, normal)
bmv_cmbo        = zip(bmv_compressed, bmv_normal)
bmv_diff        = map(x -> hdist(x[1],x[2]), bmv_cmbo)
@test sum(bmv_diff) <= 8

dct_compressed  = map(x -> perceptual_hash(x, dct_hash), compressed)
dct_normal      = map(x -> perceptual_hash(x, dct_hash), normal)
dct_cmbo        = zip(dct_compressed, dct_normal)
dct_diff        = map(x -> hdist(x[1],x[2]), dct_cmbo)
@test sum(dct_diff) <= 2

resdict = resolve_images("images", 5)

#16 mappings
@test length(resdict) == 16

#Each unique
for (k,v) in resdict
  @test length(v) == 1
end

ext: jl


text:
const DEFAULT_TYPE_MAP = Dict{Symbol,Symbol}(
    :Float => :Float64, # graphql Float is Float64
    :Int => :Int, # graphql Int is Int
    :Boolean => :Bool,
    :ID => :String, # serialized as a String but not intended to be human-readable
)

const BUILTIN_GRAPHQL_TYPES = Set{Symbol}([:ID, :Int, :String, :Float, :Boolean])

function get_schema_types(doc::Document)
    schema_types = Dict{Symbol,Symbol}()
    for d in doc
        if d isa SchemaDefinition
            merge!(schema_types, get_schema_types(d))
        end
    end
    return schema_types
end

function get_schema_types(sd::SchemaDefinition)
    schema_types = Dict{Symbol,Symbol}()
    for o in sd.operation_type_definitions
        if o.operation_type === "query" ||
           o.operation_type === "mutation" ||
           o.operation_type === "subscription"
            schema_types[jltype(o.named_type.type)] = Symbol(o.operation_type)
        end
    end
    return schema_types
end

function tojl(doc::Document, scalar_type_map::Dict)
    schema_types = get_schema_types(doc)
    doc, revisited_graph = dagify(doc)
    types = Expr[]
    functions = Expr[]

    for d in doc
        if d isa TypeDefinition &&
           d.type isa ObjectTypeDefinition &&
           haskey(schema_types, jltype(d.type.name))
            jlfuncs = jlfunction(d, schema_types)
            for f in jlfuncs
                push!(functions, f)
            end
        else
            jlt = jltype(d, revisited_graph, scalar_type_map)
            if !isnothing(jlt)
                push!(types, jlt)
            end
        end
    end

    return types, functions
end
tojl(doc::Document) = tojl(doc, Dict())

jltype(x) = nothing

function jltype(
    t::TypeDefinition,
    revisited_graph::Dict{Symbol,Set{Symbol}},
    scalar_type_map::Dict,
)
    typ = t.type
    ex = if typ isa ScalarTypeDefinition
        jltype(typ, scalar_type_map)
    elseif typ isa ObjectTypeDefinition
        jltype(typ, revisited_graph)
    elseif typ isa InputObjectTypeDefinition
        jltype(typ, revisited_graph)
    else
        jltype(typ)
    end

    docstr = isnothing(something(t.description)) ? "" : jltype(t.description)
    docstr == "" && return ex
    docstr = strip(docstr)

    ex =
        if typ isa ScalarTypeDefinition ||
           typ isa UnionTypeDefinition ||
           typ isa EnumTypeDefinition
            :(Core.@doc $docstr $ex)
        else
            quote
                $("\"\"\"\n$docstr\n\"\"\"")
                $(ex.args...)
            end
        end

    return ex
end

function jltype(name::Symbol, fields::Vector{JLKwField}, graph::Dict{Symbol,Set{Symbol}})
    ex = if haskey(graph, name)
        get_property_expr = generate_custom_getproperty!(fields, name, graph[name])
        st = JLKwStruct(; ismutable = true, name = name, fields = fields)
        quote
            $(codegen_ast(st))

            StructTypes.StructType(::Type{$name}) = StructTypes.Mutable()

            StructTypes.omitempties(::Type{$name}) = true

            $(codegen_ast(get_property_expr))
        end
    else
        st = JLKwStruct(; ismutable = true, name = name, fields = fields)
        quote
            $(codegen_ast(st))

            StructTypes.StructType(::Type{$name}) = StructTypes.Mutable()
            StructTypes.omitempties(::Type{$name}) = true
        end
    end

    return ex
end

function jltype(t::ObjectTypeDefinition, graph::Dict{Symbol,Set{Symbol}})
    name = jltype(t.name)
    fields = map(t.fields_definition) do fd
        jltype(fd)
    end

    return jltype(name, fields, graph)
end

jlfunction(t::TypeDefinition, schema_types::Dict{Symbol,Symbol}) =
    jlfunction(t.type, schema_types)

function jlfunction(t::ObjectTypeDefinition, schema_types::Dict{Symbol,Symbol})
    name = jltype(t.name)
    functions = map(t.fields_definition) do fd
        jlfunction(fd, schema_types[name])
    end

    return functions
end

function jltype(t::InputObjectTypeDefinition, graph::Dict{Symbol,Set{Symbol}})
    name = jltype(t.name)
    fields = map(t.fields) do f
        jltype(f)
    end

    return jltype(name, fields, graph)
end

function jltype(t::FieldDefinition)
    name = jltype(t.name)
    typ = jltype(t.type)

    kw = if !t.type.non_null
        JLKwField(; name = name, type = typ, default = :nothing)
    else
        JLKwField(; name = name, type = typ)
    end

    if !isnothing(something(t.description))
        kw.doc = jltype(t.description)
    end

    return kw
end

function jltype(t::ArgumentsDefinition)
    map(t.input_value_definitions) do inp
        jltype(inp)
    end
end

function jlfunctionarg(t::ArgumentsDefinition)
    map(t.input_value_definitions) do inp
        jlfunctionarg(inp)
    end
end

function jlfunction(t::FieldDefinition, stype::Symbol)
    name = jltype(t.name)
    typ = jltype(t.type)
    argdefs = if isnothing(t.arguments_definition)
        Expr[]
    else
        jlfunctionarg(t.arguments_definition)::Vector{Expr}
    end

    args = filter(e -> e.head != :kw, argdefs)
    kwargs = filter(e -> e.head == :kw, argdefs)

    # Recover the original graphql syntax so that
    # the function signature can be used in the graphql response
    funcsigs = if isnothing(t.arguments_definition)
        []
    else
        gqlstr(t.arguments_definition)
    end
    # NOTE: not sure if keyword arguments need to come last
    sort!(funcsigs, by = s -> s.default)

    sig = join(map(funcsigs) do s
        v = "$(s.name): $(s.typ)"
        if s.default != ""
            v *= " = $(s.default)"
        end
        v
    end, ", ")

    input = join(map(funcsigs) do s
        "$(s.name): \$$(s.name)"
    end, ", ")

    variable_args = map(args) do a
        e = a.args[1]
        :($("$e") => $e)
    end
    variable_kw = map(kwargs) do kw
        e = kw.args[1].args[1]
        :($("$e") => $e)
    end

    body = quote
        q = inp -> begin
                s = $("""
                $(lowercase(string(stype))) $(uppercasefirst(string(name)))($sig) {
                    $(name)($input) {
                """)
                s *= inp
                s *= """
                    }
                }
                """
            end

        query = q(f.query)
        variables = Dict($(variable_args...), $(variable_kw...))
        filter!(v -> !isnothing(v[2]), variables)

        return (; query, variables)
    end

    doc::Union{Nothing,String} = if isnothing(something(t.description))
        nothing
    else
        strip(jltype(t.description))
    end
    jlf = JLFunction(; name = :(f::$name), args, kwargs, body)

    ex = if isnothing(doc)
        quote
            struct $name
                query::String
            end

            $(codegen_ast(jlf))
        end
    else
        quote
            struct $name
                query::String
            end

            $("\"\"\"\n$doc\n\"\"\"")
            $(codegen_ast(jlf))
        end
    end

    return ex
end

function jltype(t::InputValueDefinition)
    name = jltype(t.name)
    typ = jltype(t.type)

    kw = if !t.type.non_null
        JLKwField(; name = name, type = typ, default = jltype(t.default_value))
    else
        JLKwField(; name = name, type = typ)
    end

    if !isnothing(something(t.description))
        kw.doc = jltype(t.description)
    end

    return kw
end

function jlfunctionarg(t::InputValueDefinition)
    name = jltype(t.name)
    typ = jltype(t.type)

    if !t.type.non_null
        Expr(:kw, Expr(:(::), name, typ), jltype(t.default_value))
    else
        Expr(:(::), name, typ)
    end
end

function jltype(t::UnionTypeDefinition)
    name = jltype(t.name)
    types = map(something(t.types)) do tt
        jltype(tt.type)
    end
    lhs = Expr(:curly, :Union, types...)
    return Expr(:const, Expr(:(=), name, lhs))
end

function jltype(t::EnumTypeDefinition)
    name = jltype(t.name)
    enums = map(jltype, t.enums)
    ex = quote
        @enum $name begin
            $(enums...)
        end
    end
    return ex
end

function jltype(t::EnumValueDefinition)
    ex = jltype(t.value)
    if !isnothing(something(t.description))
        doc = jltype(t.description)
        ex = :(Core.@doc $doc $ex)
    end
    return ex
end

function jltype(t::ScalarTypeDefinition, scalar_type_map::Dict)
    name = jltype(t.name)
    typ = get(scalar_type_map, name, Any)
    return Expr(:const, Expr(:(=), name, typ))
end

function jltype(t::NamedType)
    typ = jltype(t.type)
    typ = get(DEFAULT_TYPE_MAP, typ, typ)
    if !t.non_null
        :(Union{$typ,Missing,Nothing})
    else
        typ
    end
end

function jltype(t::ListType)
    typ = :(Vector{$(jltype(t.type))})
    if !t.non_null
        :(Union{$typ,Missing,Nothing})
    else
        typ
    end
end
jltype(t::DefaultValue) = jltype(t.value)
jltype(t::RBNF.Token) = Symbol(t.str)
jltype(t::RBNF.Token{:int_value}) = Base.parse(Int32, t.str)
jltype(t::RBNF.Token{:float_value}) = Base.parse(Float64, t.str)
jltype(t::RBNF.Token{:single_quote_string_value}) = convert(String, t)
jltype(t::RBNF.Token{:triple_quote_string_value}) = convert(String, t)
jltype(t::Some) = jltype(something(t))
jltype(t::Nothing) = :nothing

"""
Returns the string GraphQL representation of the parsed GraphQL type.

In a sense it's a "reverse parse".
"""
function gqlstr end

function gqlstr(t::ArgumentsDefinition)
    map(t.input_value_definitions) do iv
        gqlstr(iv)
    end
end

function gqlstr(t::InputValueDefinition)
    name = gqlstr(t.name)
    typ = gqlstr(t.type)
    default = gqlstr(t.default_value)
    (; name, typ, default)
end

function gqlstr(t::NamedType)
    s = gqlstr(t.type)
    t.non_null && (s *= "!")
    s
end

function gqlstr(t::ListType)
    s = "[$(gqlstr(t.type))]"
    t.non_null && (s *= "!")
    s
end
gqlstr(t::DefaultValue) = gqlstr(t.value)
gqlstr(t::RBNF.Token) = t.str
gqlstr(::Nothing) = ""

"""
    generate_custom_getproperty!(fields::Vector{JLKwField}, typename::Symbol, cyclic_types::Set{Symbol})

Generates a custom `Base.getproperty` function for the type `typename`. This is done
in effort to preserve type inference of fields there the type has to be removed due to
defining types in a cyclic order.

    !!! The `fields` argument is mutated by removing type information for cyclic fields.
"""
function generate_custom_getproperty!(
    fields::Vector{JLKwField},
    typename::Symbol,
    cyclic_types::Set{Symbol},
)
    fnames = Symbol[]
    ftypes = Union{Expr,Symbol}[]

    for (i, f) in enumerate(fields)
        fname = f.name
        ftype = f.type
        fdefault = f.default

        ft = get_leaf_type(ftype)
        ft in cyclic_types || continue

        push!(fnames, fname)
        push!(ftypes, ftype)

        fields[i] = JLKwField(; name = fname, default = fdefault)
    end

    length(fnames) == 0 && return

    ifexpr = JLIfElse()
    for i in 1:length(fnames)
        # using symbols with :symbol syntax is tricky in exprs
        ifexpr[:(s === Symbol($("$(fnames[i])")))] =
            :(getfield(t, (Symbol($("$(fnames[i])"))))::$(ftypes[i]))
    end
    ifexpr.otherwise = :(getfield(t, s))

    ex = quote
        function Base.getproperty(t::$typename, sym::Symbol)
            $(codegen_ast(ifexpr))
        end
    end

    return ex
end

function get_leaf_type(ex::Union{Expr,Symbol})
    ex isa Symbol && return ex

    t = if ex.args[1] == :Vector
        get_leaf_type(ex.args[2])
    else
        idx = findfirst(ft -> !in(ft, (:Union, :Missing, :Nothing)), ex.args)
        get_leaf_type(ex.args[idx])
    end

    return t
end


text:
struct Toy2DContState
    x::Float64
    y::Float64
end
norm(s::Toy2DContState) = LinearAlgebra.norm(SVector{2,Float64}(s.x,s.y))

struct Toy2DContAction
    xdot::Float64
    ydot::Float64
end
norm(a::Toy2DContAction) = LinearAlgebra.norm(SVector{2,Float64}(a.xdot,a.ydot))

# Note - Hardcoding 1 x 1 grid because it can be scale-invariant
struct Toy2DParameters
    speed_limit::Float64
    speed_vals::Int64
    epsilon::Float64
    num_generative_samples::Int64
    num_bridge_samples::Int64
    horizon_limit::Int64
    axis_vals::Int64
    valid_vertices_threshold::Float64
    time_coefficient::Float64
end

# Const types
const Toy2DStateType = CMSSPState{Int64,Toy2DContState}
const Toy2DActionType = CMSSPAction{Int64,Toy2DContAction}
const Toy2DCMSSPType = CMSSP{Int64,Toy2DContState,Int64,Toy2DContAction,Toy2DParameters}
const Toy2DModalMDPType = ModalFinHorMDP{Int64,Toy2DContState,Toy2DContAction,Toy2DParameters}
const Toy2DOpenLoopVertex = OpenLoopVertex{Int64,Toy2DContState,Int64}
# Context is current collection of all valid transition points at each timestep
const Toy2DModePair = Tuple{Int64, Int64}
const Toy2DContextType = Vector{Tuple{Toy2DModePair,Toy2DContState}}

# Const values
const TOY2D_MODES = [1,2,3,4,5]
const TOY2D_GOAL_MODE = 5
const TOY2D_QUADRANTS = [(0.0,pi/2), (pi/2,pi), (pi,3*pi/2), (3*pi/2,2*pi)]

mutable struct Toy2DContextSet
    curr_time::Int64
    switch_bias_dict::Dict{Toy2DModePair,Int64}
    curr_context_set::Vector{Toy2DContextType}
end

function Toy2DContextSet()
    return Toy2DContextSet(0,Dict{Toy2DModePair,Int64}(),Vector{Toy2DContextType}(undef,0))
end

const Toy2DSolverType = HHPCSolver{Int64,Toy2DContState,Int64,Toy2DContAction,Toy2DContextSet,Toy2DParameters}

# Hardcoded - only 1 modeswitch action
# Get control actions based on parameters
function get_toy2d_actions(params::Toy2DParameters)

    actions = Vector{Toy2DActionType}(undef,0)

    idx = 1
    push!(actions, Toy2DActionType(1,idx))
    idx = 2
    push!(actions, Toy2DActionType(2,idx))

    vel_vals = polyspace_symmetric(params.speed_limit, params.speed_vals)

    # Now enter control actions based on parameters
    #vel_vals = range(-params.speed_limit, stop=params.speed_limit, step=params.speed_resolution)
    for xdot in vel_vals
        for ydot in vel_vals
            idx += 1
            push!(actions, Toy2DActionType(Toy2DContAction(xdot,ydot), idx))
        end
    end

    return actions
end

# Hardcoded tabular MDP
# 1 -- 2 -- 5
# |         |
# 3 -- 4 -- |
function get_toy2d_switch_mdp()

    # Initialize T as 5 x 1 x 5
    T = zeros(5, 2, 5)
    T[1,1,2] = 1.0
    T[1,2,3] = 1.0
    T[2,1,5] = 1.0
    T[3,1,4] = 1.0
    T[3,2,5] = 1.0
    T[4,1,5] = 1.0

    R = zeros(5,2)
    R[1,1] = -1.0 # For T[1,1,2]
    R[1,2] = -2.0 # For T[1,2,3]
    R[2,1] = -10.0
    R[3,1] = -2.0
    R[3,2] = -8.0
    R[4,1] = -1.0

    # Undiscounted
    return TabularMDP(T, R, 1.0)
end

# Create Toy2DCMSSP based on params
function create_toy2d_cmssp(params::Toy2DParameters)

    actions = get_toy2d_actions(params)
    modes = TOY2D_MODES
    switch_mdp = get_toy2d_switch_mdp()

    return Toy2DCMSSPType(actions, modes, switch_mdp, params)
end

function POMDPs.isterminal(mdp::Toy2DModalMDPType, relative_state::Toy2DContState)
    return norm(relative_state) < mdp.params.epsilon
end

function POMDPs.isterminal(cmssp::Toy2DCMSSPType, state::Toy2DStateType)
    return state.mode == TOY2D_GOAL_MODE
end

function rel_state(source::Toy2DContState, target::Toy2DContState)
    return Toy2DContState(source.x - target.x, source.y - target.y)
end

function HHPC.get_relative_state(mdp::Toy2DModalMDPType, source::Toy2DContState, target::Toy2DContState)
    return rel_state(source, target)
end

function POMDPs.convert_s(t::Type{V}, c::Toy2DContState, mdp::Toy2DModalMDPType) where V <: AbstractVector{Float64}
    v_ = SVector{2,Float64}(c.x, c.y)
    v = convert(t,v_)
    return v
end

function POMDPs.convert_s(::Type{Toy2DContState}, v::AbstractVector{Float64}, mdp::Toy2DModalMDPType)
    c = Toy2DContState(v[1], v[2])
    return c
end

function sample_toy2d(rng::RNG=Random.GLOBAL_RNG) where {RNG <: AbstractRNG}
    d = Uniform(0.0,1.0)
    return Toy2DContState(rand(rng,d), rand(rng,d))
end

function generate_start_state(cmssp::Toy2DCMSSPType, rng::RNG) where {RNG <: AbstractRNG}
    return Toy2DStateType(1, sample_toy2d(rng))
end

# Just used by generate_sr - euclidean distance + squared velocity cost
function POMDPs.reward(mdp::Toy2DModalMDPType, state::Toy2DContState,
                       cont_action::Toy2DContAction, statep::Toy2DContState)
    delta_state = Toy2DContState(statep.x - state.x, statep.y - state.y)
    reward = -1.0*(norm(delta_state)) - mdp.params.time_coefficient
    return reward
end

function HHPC.expected_reward(mdp::Toy2DModalMDPType, state::Toy2DContState,
                         cont_action::Toy2DContAction, rng::RNG=Random.GLOBAL_RNG) where {RNG <: AbstractRNG}
    params = mdp.params
    avg_rew = 0.0
    for i = 1:params.num_generative_samples
        (_, r) = next_state_reward(mdp,state, cont_action, rng)
        avg_rew += r
    end
    avg_rew /= params.num_generative_samples
    return avg_rew
end

function next_state_reward(mdp::Toy2DModalMDPType, state::Toy2DContState,
                           cont_action::Toy2DContAction, rng::RNG=Random.GLOBAL_RNG) where {RNG <: AbstractRNG}

    noise_dist = MvNormal([0.0001 + cont_action.xdot/100.0, 0.0001 + cont_action.ydot/100.0])
    xy_noise = rand(rng,noise_dist)
    new_x = clamp(state.x + cont_action.xdot + xy_noise[1], -1.0, 1.0)
    new_y = clamp(state.y + cont_action.ydot + xy_noise[2], -1.0, 1.0)

    new_state = Toy2DContState(new_x, new_y)
    reward_val = reward(mdp, state, cont_action, new_state)

    return (new_state, reward_val)
end

function next_state_reward_sim(mdp::Toy2DModalMDPType, state::Toy2DContState,
                           cont_action::Toy2DContAction, rng::RNG=Random.GLOBAL_RNG) where {RNG <: AbstractRNG}

    noise_dist = MvNormal([0.0001 + cont_action.xdot/100.0, 0.0001 + cont_action.ydot/100.0])
    xy_noise = rand(rng,noise_dist)
    new_x = clamp(state.x + cont_action.xdot + xy_noise[1], 0.0, 1.0)
    new_y = clamp(state.y + cont_action.ydot + xy_noise[2], 0.0, 1.0)

    new_state = Toy2DContState(new_x, new_y)
    reward_val = reward(mdp, state, cont_action, new_state)

    return (new_state, reward_val)
end

# In this case, independent of actual mode
function POMDPs.generate_sr(mdp::Toy2DModalMDPType, state::Toy2DContState,
                            cont_action::Toy2DContAction, rng::RNG=Random.GLOBAL_RNG) where {RNG <: AbstractRNG}
    # can be used with relative state too
    return next_state_reward(mdp, state, cont_action, rng)
end

# Needs to be bound to context set in runtime script
# We want that the state of the vertex be within the bounds
# of the transition point. If it is, set that to be the horizon
# Otherwise, set vertex to t = 0 with probability 1.0
function HHPC.update_vertices_with_context!(cmssp::Toy2DCMSSPType, vertices::Vector{Toy2DOpenLoopVertex},
                                       switch::Toy2DModePair, toy2d_context_set::Toy2DContextSet)

    # Note that this is for a specific mode, so only consider those with the same mode switch
    # Iterate over context set
    # If the mode switch tuple is the same as in the open-loop vertex
    # Then check that the pre_bridge_state is in the range of the first
    # The post bridge state is just the same as from where the mode switch is attempted
    # @assert (vertices[1].pre_bridge_state.mode, vertices[1].state.mode) == switch "Switch argument does not match OpenLoop Vertices"

    # # IMP - Only look at future context, so start with index 2
    # context_hor = length(toy2d_context_set)

    # num_verts = length(vertices)
    # num_valid_updates = 0

    # # TODO : Search more efficiently by starting from current v.tp???
    # for v in vertices

    #     tp_time = Vector{Int64}(undef,0)
    #     tp_prob = Vector{Float64}(undef,0)

    #     for (hor,context) in enumerate(toy2d_context_set[2:end])

    #         for switch_points in context
    #             if switch_points[1] != switch
    #                 continue
    #             end
    #             pt_dist = norm(rel_state(v.pre_bridge_state, switch_points[2]))
    #             if pt_dist < params.epsilon
    #                 push!(tp_time,hor)
    #                 push!(tp_prob,1.0/pt_dist)
    #             end
    #         end
    #     end

    #     # Will either be [0,1.0] or a weighted probability dist of inverse distance
    #     if isempty(tp_time)
    #         v.tp = TPDistribution([0],[1.0])
    #     else
    #         tp_prob = tp_prob/sum(tp_prob)
    #         v.tp = TPDistribution(tp_time, tp_prob)
    #         num_valid_updates += 1
    #     end
    # end

    # # If too many vertices invalidated, then return false and resample
    # return (num_valid_updates > params.valid_vertices_threshold*num_verts)

    # OR : JUST RETURN FALSE EACH TIME AND RE-SAMPLE
    return false
end

# Just generate points around centre within +- box
function HHPC.generate_goal_sample_set(cmssp::Toy2DCMSSPType, popped_cont::Toy2DContState,
                                  num_samples::Int64, rng::RNG=Random.GLOBAL_RNG) where {RNG <: AbstractRNG}

    # params = cmssp.params
    # goal_samples = Vector{Toy2DStateType}(undef,0)

    # dx = Uniform(TOY2D_GOAL_CENTRE.x - params.epsilon/2.0, TOY2D_GOAL_CENTRE.x + params.epsilon/2.0)
    # dy = Uniform(TOY2D_GOAL_CENTRE.y - params.epsilon/2.0, TOY2D_GOAL_CENTRE.y + params.epsilon/2.0)

    # for i = 1:num_samples
    #     push!(goal_samples, Toy2DStateType(TOY2D_GOAL_MODE, Toy2DContState(rand(rng,dx), rand(rng,dy))))
    # end
    goal_samples = [Toy2DStateType(TOY2D_GOAL_MODE, popped_cont)]

    return goal_samples
end

function HHPC.generate_next_valid_modes(cmssp::Toy2DCMSSPType, mode::Int64, toy2d_context_set::Toy2DContextSet)

    action_nextmodes = Vector{Tuple{Int64,Int64}}(undef,0)
    next_modes = Set{Int64}()
    curr_context_set = toy2d_context_set.curr_context_set

    for context in curr_context_set[2:end]

        for switch_points in context
            if switch_points[1][1] == mode
                next_mode = switch_points[1][2]

                if next_mode in next_modes
                    continue
                end

                mode_idx = mode_index(cmssp, mode)
                next_mode_idx = mode_index(cmssp, next_mode)

                # Need to find action to take to next mode
                for (ac_idx,action) in enumerate(cmssp.mode_actions)
                    if cmssp.modeswitch_mdp.T[mode_idx,ac_idx,next_mode_idx] > 0.0
                        push!(action_nextmodes,(ac_idx,next_mode_idx))
                        push!(next_modes, next_mode)
                        break
                    end
                end
            end
        end
    end
    return action_nextmodes
end

function HHPC.generate_bridge_sample_set(cmssp::Toy2DCMSSPType, cont_state::Toy2DContState,
                                    mode_pair::Toy2DModePair, num_samples::Int64,
                                    toy2d_context_set::Toy2DContextSet, rng::RNG=Random.GLOBAL_RNG) where {RNG <: AbstractRNG}

    params = cmssp.params
    bridge_samples = Vector{BridgeSample{Toy2DContState}}(undef,0)
    curr_context_set = toy2d_context_set.curr_context_set
    curr_time = toy2d_context_set.curr_time

    avg_samples_per_context = convert(Int64,ceil(num_samples/length(curr_context_set)))

    for (hor,context) in enumerate(curr_context_set[2:end])

        for switch_points in context

            if switch_points[1] == mode_pair
                # Generate bridge points around switch
                dx = Uniform(switch_points[2].x - params.epsilon/2.0, switch_points[2].x + params.epsilon/2.0)
                dy = Uniform(switch_points[2].y - params.epsilon/2.0, switch_points[2].y + params.epsilon/2.0)

                for _ = 1:avg_samples_per_context
                    bpx = clamp(rand(rng,dx),0.0,1.0)
                    bpy = clamp(rand(rng,dy),0.0,1.0)
                    state = Toy2DContState(bpx,bpy)
                    push!(bridge_samples, BridgeSample(state, state, TPDistribution([curr_time+hor],[1.0])))
                end
            end
        end
    end
    return bridge_samples
end

"""
Returns next state, reward, and if a mode switch failed
"""

function HHPC.simulate_cmssp(cmssp::Toy2DCMSSPType, state::Toy2DStateType, a::Toy2DActionType, t::Int64,
                        toy2d_context_set::Toy2DContextSet, rng::RNG=Random.GLOBAL_RNG) where {RNG <: AbstractRNG}

    curr_context_set = toy2d_context_set.curr_context_set
    curr_context = curr_context_set[1]

    curr_mode = state.mode
    curr_cont_state = state.continuous
    params = cmssp.params

    if typeof(a.action) <: Int64
        # action is mode switch - attempt
        mode_idx = mode_index(cmssp, curr_mode)
        mode_ac_idx = mode_actionindex(cmssp, a.action)
        next_mode_idx = findfirst(x -> x > 0.0, cmssp.modeswitch_mdp.T[mode_idx, mode_ac_idx,:])
        next_mode = cmssp.modes[next_mode_idx]

        # Action cost incurred no matter what
        reward = cmssp.modeswitch_mdp.R[mode_idx, mode_ac_idx]

        # Check if it can be done based on context
        for switch_points in curr_context
            if switch_points[1] == (curr_mode,next_mode)
                if norm(rel_state(curr_cont_state,switch_points[2])) < 2.5*params.epsilon
                    @info "Successful mode switch to mode ",next_mode
                    next_state = Toy2DStateType(next_mode, curr_cont_state)
                    return (next_state, reward, false)
                else
                    @show switch_points[2]
                    @info "Mode switch failed!"
                    next_state = curr_state
                    return (next_state,  reward, true)
                end
            end
        end

        @info "Attempted mode switch unavailable from context!"
        next_state = curr_state
        return (new_state, reward, true)
    else
        # Control action within mode
        new_cont_state, reward = next_state_reward_sim(Toy2DModalMDPType(curr_mode, cmssp.params), curr_cont_state, a.action, rng)
        next_state = Toy2DStateType(curr_mode, new_cont_state)
        return (next_state, reward, false)
    end
end

function toy2d_parse_params(filename::AbstractString)

    params_key = TOML.parsefile(filename)

    return Toy2DParameters(params_key["SPEED_LIMIT"],
                           params_key["SPEED_VALS"],
                           params_key["EPSILON"],
                           params_key["NUM_GEN_SAMPLES"],
                           params_key["NUM_BRIDGE_SAMPLES"],
                           params_key["HORIZON_LIMIT"],
                           params_key["AXIS_VALS"],
                           params_key["VALID_VERTICES_THRESHOLD"],
                           params_key["TIME_COEFFICIENT"])
end

# Context functions
function generate_start_context_set(cmssp::Toy2DCMSSPType, rng::RNG=Random.GLOBAL_RNG) where {RNG <: AbstractRNG}

    params = cmssp.params
    toy2d_context_set = Toy2DContextSet()
    init_context = Toy2DContextType(undef,0)

    # Maintains the context bias for each mode switch
    switch_bias_dict = toy2d_context_set.switch_bias_dict

    for m1 in cmssp.modes
        for m2 in cmssp.modes
            for a in cmssp.mode_actions
                if cmssp.modeswitch_mdp.T[m1,a,m2] > 0.0
                    switch_point = sample_toy2d(rng)
                    push!(init_context,(Toy2DModePair((m1,m2)), switch_point))
                    switch_bias_dict[Toy2DModePair((m1,m2))] = rand(rng, 1:4)
                end
            end
        end
    end

    push!(toy2d_context_set.curr_context_set, init_context)

    r = params.epsilon/2.0

    # Now have init_context, need to propagate forward and keep pushing
    for idx = 1:params.horizon_limit
        temp_context = toy2d_context_set.curr_context_set[idx]

        next_context = Toy2DContextType(undef,0)

        # Iterate over each and propagate with epsilon/2.0 in some radial direction
        for (switch,point) in temp_context
            theta_dist = Uniform(TOY2D_QUADRANTS[switch_bias_dict[switch]][1], TOY2D_QUADRANTS[switch_bias_dict[switch]][2])
            theta = rand(rng,theta_dist)
            new_x = clamp(point.x + r*cos(theta),0.0,1.0)
            new_y = clamp(point.y + r*sin(theta),0.0,1.0)
            new_point = Toy2DContState(new_x,new_y)

            push!(next_context,(switch,new_point))

            # If either outside range, switch to opposite TOY2D_QUADRANTS
            if new_y*new_x == 0.0 || new_x == 1.0 || new_y == 1.0
                switch_bias_dict[switch] = (switch_bias_dict[switch] + 2) % 4 + 1
            end
        end

        # Now append to context
        push!(toy2d_context_set.curr_context_set, next_context)
    end

    @assert length(toy2d_context_set.curr_context_set) == params.horizon_limit+1
    return toy2d_context_set
end

function HHPC.update_context_set!(cmssp::Toy2DCMSSPType, toy2d_context_set::Toy2DContextSet,
                            rng::RNG=Random.GLOBAL_RNG) where {RNG <: AbstractRNG}

    switch_bias_dict = toy2d_context_set.switch_bias_dict
    curr_context_set = toy2d_context_set.curr_context_set
    params = cmssp.params

    # For each switch, with 50% probability, have it either stay in place
    # Or go to the next point. In the first case, just maintain the same vector
    # Otherwise, perturb the last point and add, and left-shift
    next_context = Toy2DContextType(undef,0)
    last_context = Toy2DContextType(undef,0)

    curr_context = curr_context_set[1]
    curr_next_context = curr_context_set[2]
    curr_last_context = curr_context_set[end]

    theta_dist = Uniform(0,2*pi)
    r = params.epsilon/2.0

    for (i,switch_point) in enumerate(curr_context)
        toss = rand(rng,Uniform(0.0,1.0))
        switch = switch_point[1]
        point = switch_point[2]

        if toss < 0.35
            # Stay and most_future is the same
            push!(next_context,(switch,point))
            push!(last_context,(switch,curr_last_context[i][2]))
        else
            # Go to next one
            push!(next_context,(switch,curr_next_context[i][2]))

            # Generate new point for last context
            curr_last_context_pt = curr_last_context[i][2]
            theta_dist = Uniform(TOY2D_QUADRANTS[switch_bias_dict[switch]][1], TOY2D_QUADRANTS[switch_bias_dict[switch]][2])
            theta = rand(rng,theta_dist)
            new_x = clamp(curr_last_context_pt.x + r*cos(theta),0.0,1.0)
            new_y = clamp(curr_last_context_pt.y + r*sin(theta),0.0,1.0)
            new_point = Toy2DContState(new_x,new_y)

            push!(last_context,(switch,new_point))

            # If either outside range, switch to opposite TOY2D_QUADRANTS
            if new_y*new_x == 0.0 || new_x == 1.0 || new_y == 1.0
                switch_bias_dict[switch] = (switch_bias_dict[switch] + 2) % 4 + 1
            end
        end
    end

    # Now re-assign
    toy2d_context_set.curr_context_set[1] = next_context
    toy2d_context_set.curr_context_set[end] = last_context

    # Finally, update time
    toy2d_context_set.curr_time = toy2d_context_set.curr_time + 1
end

function HHPC.display_context_future(toy2d_context_set::Toy2DContextSet,future_time::Int64)

    # IMP - Because if 0, then first
    hor = future_time - toy2d_context_set.curr_time + 1
    println(toy2d_context_set.curr_context_set[hor])
end

## TODO : NEED TO IMPLEMENT generate_sr for MCTS


text:
using StaticArrays

# Frequently-used array types
const Mat3{T} = SMatrix{3, 3, T, 9} where T
const Vec3{T} = SVector{3, T} where T
const AbstractArray3{T} = AbstractArray{T, 3}

# Represents a symmetry (S,τ)
const SymOp = Tuple{Mat3{Int}, Vec3{Float64}}
identity_symop() = (Mat3{Int}(I), Vec3(zeros(3)))


text:
# ---
# title: 94. Binary Tree Inorder Traversal
# id: problem94
# author: zhwang
# date: 2022-02-19
# difficulty: Medium
# categories: Hash Table, Stack, Tree
# link: <https://leetcode.com/problems/binary-tree-inorder-traversal/description/>
# hidden: true
# ---
#
# Given the `root` of a binary tree, return _the inorder traversal of its nodes
# ' values_.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/09/15/inorder_1.jpg)
#
#
#
#     Input: root = [1,null,2,3]
#     Output: [1,3,2]
#
#
# **Example 2:**
#
#
#
#     Input: root = []
#     Output: []
#
#
# **Example 3:**
#
#
#
#     Input: root = [1]
#     Output: [1]
#
#
# **Example 4:**
#
# ![](https://assets.leetcode.com/uploads/2020/09/15/inorder_5.jpg)
#
#
#
#     Input: root = [1,2]
#     Output: [2,1]
#
#
# **Example 5:**
#
# ![](https://assets.leetcode.com/uploads/2020/09/15/inorder_4.jpg)
#
#
#
#     Input: root = [1,null,2]
#     Output: [1,2]
#
#
#
#
# **Constraints:**
#
#   * The number of nodes in the tree is in the range `[0, 100]`.
#   * `-100 <= Node.val <= 100`
#
#
#
# **Follow up:**
#
# Recursive solution is trivial, could you do it iteratively?
#
#
#
#
## @lc code=start
using LeetCode

inorder_traversal(::Nothing) = Int[]
function inorder_traversal(root::TreeNode)::Vector{Int}
    res, stack = Int[], []
    while !isempty(stack) || !isnothing(root)
        while !isnothing(root)
            push!(stack, root)
            root = root.left
        end
        root = pop!(stack)
        push!(res, root.val) ## inorder search
        root = root.right
    end
    res
end

##### code template for inorder traversal and preorder traversal #####
## traversal(::Nothing) = Int[]
## function traversal(root::TreeNode)::Vector{Int}
##     res, stack = Int[], []
##     while !isempty(stack) || !isnothing(root)
##         while !isnothing(root)
##             push!(stack, root)
##             ## preorder traversal: put codes here
##             root = root.left
##         end
##         root = pop!(stack)
##         ## inorder traversal: put codes here
##         root = root.right
##     end
##     res
## end

## @lc code=end
## @lc test=start
@testset "94.binary-tree-inorder-traversal.jl" begin
    @test inorder_traversal(TreeNode{Int}([1])) == [1]
    @test inorder_traversal(nothing) == Int[]
    @test inorder_traversal(TreeNode{Int}([1,nothing,2,3])) == [1,3,2]
    @test inorder_traversal(TreeNode{Int}([1, 2, 3, 4, 5, 6])) == [4,2,5,1,6,3]
    @test inorder_traversal(TreeNode{Int}([3,9,4,nothing, nothing, 5,7])) == [9,3,5,4,7]
end
## @lc test=end

## @lc info=start
# link: [solution to 94](https://leetcode-cn.com/problems/binary-tree-preorder-traversal/solution/144-er-cha-shu-de-qian-xu-bian-li-by-rex-frkw/)
## @lc info=end


text:
using Plots,JLD
function randomWalk(starting_position,L)
    mean_life_span=0
    TotData=[[0.0 for i in 1:L]]
    TotData[1][starting_position]=1
    death_p_sum=(TotData[1][1]+TotData[1][end])
    if death_p_sum>=0.9 && death_p_sum<=1
        return 0
    end
    step=2
    while true

        death_p_sum=0
        R=[0.0 for e in 1:L]
        for i in 2:(L-1)
            R[i-1:2:i+1].+=TotData[step-1][i]/2
        end
        R[1]+=TotData[step-1][1]
        R[end]+=TotData[step-1][end]
        push!(TotData, R)

        death_p_sum+=(R[1]+R[end])
        #println(R)
        # println(step)
        # println(R[1])
        # println(R[end])
        mean_life_span+=sum(R[2:19])
        if death_p_sum>=0.99999999
            break
        end
        step+=1
    end
    return mean_life_span
end

#time_=randomWalk(2,20)
MeanLifeSpan=[]
StartingPosList=[i for i in 1:20]
for StartingPos in StartingPosList
    println(StartingPos)
    push!(MeanLifeSpan,randomWalk(StartingPos,20))
end
save("C:\\Users\\Narges\\Documents\\GitHub\\computational_physics\\chapter5\\4.4\\MeanLifeSpan_startpos_deterministic_data.jld",
 "data", MeanLifeSpan)
scatter(StartingPosList, MeanLifeSpan,legends=false, color=:black, alpha=0.5,dpi=400)
plot!(StartingPosList, MeanLifeSpan,guidefontsize=9, color=:pink2,
 xlabel="the Starting Position",ylabel="Mean Life-Span",
 title="Mean Life-Span Over the Starting Position Using the Deterministic Algorithm", titlefontsize=9)
 savefig("C:\\Users\\Narges\\Documents\\GitHub\\computational_physics\\chapter5\\4.4\\Figs\\MeanLifeSpan_startpos_deterministic.png")


text:
"""
    legendre_1term_raise_lm_accuracy.jl

Simple tests to explore the numerical accuracy of the 1-term recurrence relation which
boosts (m,m) -> (m+1,m+1) in the associated Legendre polynomials.

## Tests

- `recur0`: Baseline version which is a direct translation of the math.
- `recur1`: Uses FMA instructions to compute (1-x^2). Eliminating intermediate rounding
  is expected to increase the accuracy of the entire recursion.
- `recur2`: Same as `recur1` but also iterates upward in degree/order in pairs. Stepping
  from even to odd ``(m,m) -> (m+1,m+1)`` multiplies by ``y = \\sqrt{1-x^2}`` and the
  coefficient ``μ_m``. The second half then steps from even to even
  ``(m,m) -> (m+2,m+2)`` by multiplying by ``y^2 = 1-x^2`` (without taking the square
  root) and the pair of coefficients ``μ_m`` and ``μ_{m+1}``.
- `recur3`: Same as `recur2` but replaces the FMA calculation of ``1-x^2`` with the
  product ``(1-x)(1+x)``.

## Conclusions
- Without FMA, the product ``(1-x)(1+x)`` is better for ``|x| ≈ 1`` than ``1 - x^2``.
- FMA is the best option for computing ``y = \\sqrt(1-x^2)`` and should be preferred.
- The pair-wise iteration improves numerical accuracy as ``m \\to \\finty``.
"""

using PyPlot, PyCall

function recur0(x, nstep)
    y = sqrt(one(x) - x^2)
    vals = zeros(typeof(x), 2nstep)
    z = inv(sqrt(4oftype(x, π)))
    @inbounds for ii = 1:2nstep
        μ = sqrt(one(x) + one(x)/2(ii+1))
        z = -μ * y * z
        vals[ii] = z
    end
    return vals
end

function recur1(x, nstep)
    y = sqrt(-fma(x, x, -one(x)))
    vals = zeros(typeof(x), 2nstep)
    z = inv(sqrt(4oftype(x, π)))
    @inbounds for ii = 1:2nstep
        μ = sqrt(one(x) + one(x)/2(ii+1))
        z = -μ * y * z
        vals[ii] = z
    end
    return vals
end

function recur2(x, nstep)
    y² = -fma(x, x, -one(x))
    y¹ = sqrt(y²)
    vals = zeros(typeof(x), 2nstep)
    z = inv(sqrt(4oftype(x, π)))
    @inbounds for ii = 1:2:2nstep
        μ₁ = sqrt(one(x) + one(x)/2(ii+1))
        μ₂ = sqrt(one(x) + one(x)/2(ii+2))
        vals[ii] = -μ₁ * y¹ * z
        z = μ₁ * μ₂ * y² * z
        vals[ii+1] = z
    end
    return vals
end

function recur3(x, nstep)
    y² = (one(x) - x) * (one(x) + x)
    y¹ = sqrt(y²)
    vals = zeros(typeof(x), 2nstep)
    z = inv(sqrt(4oftype(x, π)))
    @inbounds for ii = 1:2:2nstep
        μ₁ = sqrt(one(x) + one(x)/2(ii+1))
        μ₂ = sqrt(one(x) + one(x)/2(ii+2))
        vals[ii] = -μ₁ * y¹ * z
        z = μ₁ * μ₂ * y² * z
        vals[ii+1] = z
    end
    return vals
end

# Half of the lmax/mmax to iterate to.
nstep = 512
# Range of arguments, taken more densely from x ≈ 1 where catastrophic cancellation is
# expected without FMA.
x₀ = 1 .- exp2.(range(-16, -eps(), length=2560))

# Reference values calculated with extended precision.
ref = Float64.(reduce(hcat, setprecision(512) do
    recur1.(big.(x₀), nstep)
end))

v0 = reduce(hcat, recur0.(x₀, nstep))
v1 = reduce(hcat, recur1.(x₀, nstep))
v2 = reduce(hcat, recur2.(x₀, nstep))
v3 = reduce(hcat, recur3.(x₀, nstep))

colors = pyimport("matplotlib.colors")
fig, axs = subplots(2, 2, num = "legendre_1term_raise_lm_accuracy", clear = true)
function plot_compare(ax, v; title=nothing, kws...)
    lnorm = colors.LogNorm(vmin = 1, vmax=1e3)
    ulps = abs.(v .- ref) ./ eps.(ref)
    img = ax.pcolormesh(1 .- x₀, 1:2nstep, ulps,
                    norm=lnorm;
                    kws...)
    ax.set_xscale("log")
    ax.set_xlabel("argument \$(1-x)\$")
    ax.set_ylabel("degree/order")
    cbar = ax.figure.colorbar(img, ax=ax)
    cbar.set_label("ulps");
    title !== nothing && ax.set_title(title)
    ax.text(0.05, 0.9, "ulp max = $(round(maximum(ulps), digits=2))", transform=ax.transAxes)
    nothing
end
plot_compare(axs[1,1], v0, title="v0 — Baseline")
plot_compare(axs[1,2], v1, title="v1 — FMA \$(1 - x^2)\$")
plot_compare(axs[2,1], v2, title="v2 — FMA + pairwise iteration")
plot_compare(axs[2,2], v3, title="v3 — \$(1-x)(1+x)\$ + pairwise iteration")


text:
"""
    space(::SiteType"S=1/2";
          conserve_qns = false,
          conserve_sz = conserve_qns,
          conserve_szparity = false,
          qnname_sz = "Sz",
          qnname_szparity = "SzParity")

Create the Hilbert space for a site of type "S=1/2".

Optionally specify the conserved symmetries and their quantum number labels.
"""
function ITensors.space(
  ::SiteType"S=1/2";
  conserve_qns=false,
  conserve_sz=conserve_qns,
  conserve_szparity=false,
  qnname_sz="Sz",
  qnname_szparity="SzParity",
)
  if conserve_sz && conserve_szparity
    return [
      QN((qnname_sz, +1), (qnname_szparity, 1, 2)) => 1,
      QN((qnname_sz, -1), (qnname_szparity, 0, 2)) => 1,
    ]
  elseif conserve_sz
    return [QN(qnname_sz, +1) => 1, QN(qnname_sz, -1) => 1]
  elseif conserve_szparity
    return [QN(qnname_szparity, 1, 2) => 1, QN(qnname_szparity, 0, 2) => 1]
  end
  return 2
end

# Use Qubit  definition of any operator/state
# called using S=1/2 SiteType
function ITensors.val(vn::ValName, ::SiteType"S=1/2"; kwargs...)
  return val(vn, SiteType("Qubit"); kwargs...)
end

function ITensors.state(sn::StateName, ::SiteType"S=1/2"; kwargs...)
  return state(sn, SiteType("Qubit"); kwargs...)
end

ITensors.op(o::OpName, ::SiteType"S=1/2"; kwargs...) = op(o, SiteType("Qubit"); kwargs...)

# Support the tag "SpinHalf" as equivalent to "S=1/2"
space(::SiteType"SpinHalf"; kwargs...) = space(SiteType("S=1/2"); kwargs...)

val(name::ValName, ::SiteType"SpinHalf") = val(name, SiteType("S=1/2"))

ITensors.state(name::StateName, ::SiteType"SpinHalf") = state(name, SiteType("S=1/2"))

function ITensors.op(o::OpName, ::SiteType"SpinHalf"; kwargs...)
  return op(o, SiteType("S=1/2"); kwargs...)
end

# Support the tag "S=½" as equivalent to "S=1/2"

space(::SiteType"S=½"; kwargs...) = space(SiteType("S=1/2"); kwargs...)

val(name::ValName, ::SiteType"S=½") = val(name, SiteType("S=1/2"))

ITensors.state(name::StateName, ::SiteType"S=½") = state(name, SiteType("S=1/2"))

ITensors.op(o::OpName, ::SiteType"S=½"; kwargs...) = op(o, SiteType("S=1/2"); kwargs...)


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    PolyArea(outer, [inner1, inner2, ..., innerk])

A polygonal area with `outer` chain, and optional inner
chains `inner1`, `inner2`, ..., `innerk`.

Chains can be a vector of [`Point`](@ref) or a
vector of tuples with coordinates for convenience.

Most algorithms assume that the outer chain is oriented
counter-clockwise (CCW) and that all inner chains are
oriented clockwise (CW).
"""
struct PolyArea{Dim,T,C<:Chain{Dim,T}} <: Polygon{Dim,T}
  outer::C
  inners::Vector{C}

  function PolyArea{Dim,T,C}(outer, inners) where {Dim,T,C}
    @assert isclosed(outer) "invalid outer chain"
    @assert all(isclosed.(inners)) "invalid inner chains"

    # fix orientation if needed
    fix(c, o) = orientation(c) == o ? c : reverse(c)
    ofixed = fix(outer, :CCW)
    ifixed = map(c -> fix(c, :CW), inners)

    new(ofixed, ifixed)
  end
end

PolyArea(outer::C, inners=[]) where {Dim,T,C<:Chain{Dim,T}} =
  PolyArea{Dim,T,Chain{Dim,T}}(outer, inners)

PolyArea(outer::AbstractVector{P}, inners=[]) where {P<:Point} =
  PolyArea(Chain(outer), [Chain(inner) for inner in inners])

PolyArea(outer::AbstractVector{TP}, inners=[]) where {TP<:Tuple} =
  PolyArea(Point.(outer), [Point.(inner) for inner in inners])

PolyArea(outer::Vararg{P}) where {P<:Point} = PolyArea(collect(outer))

PolyArea(outer::Vararg{TP}) where {TP<:Tuple} = PolyArea(collect(Point.(outer)))

==(p1::PolyArea, p2::PolyArea) =
  p1.outer == p2.outer && p1.inners == p2.inners

function vertices(p::PolyArea{Dim,T}) where {Dim,T}
  vo = vertices(p.outer)
  vi = reduce(vcat, vertices(inner) for inner in p.inners; init=Point{Dim,T}[])
  [vo; vi]
end

nvertices(p::PolyArea) = nvertices(p.outer) + mapreduce(nvertices, +, p.inners, init=0)

centroid(p::PolyArea) = centroid(p.outer)

"""
    chains(polyarea)

Return the outer and inner chains of the polygon.
"""
chains(p::PolyArea) = [p.outer; p.inners]

"""
    hasholes(polyarea)

Tells whether or not the `polyarea` contains holes.
"""
hasholes(p::PolyArea) = !isempty(p.inners)

"""
    issimple(polyarea)

Tells whether or not the `polyarea` is a simple polygon.
See https://en.wikipedia.org/wiki/Simple_polygon.
"""
issimple(p::PolyArea) = !hasholes(p) && issimple(p.outer)

"""
    windingnumber(point, polyarea)

Winding number of `point` with respect to the `polyarea`.
"""
windingnumber(point::Point, p::PolyArea) =
  windingnumber(point, p.outer)

"""
    orientation(polyarea)

Returns the orientation of the `polyarea` as either
counter-clockwise (CCW) or clockwise (CW).

For polygons with holes, returns a list of orientations.
"""
orientation(p::PolyArea, algo=WindingOrientation()) =
  orientation.([p.outer; p.inners], Ref(algo))

"""
    unique(polyarea)

Return a new `polyarea` without duplicate vertices.
"""
Base.unique(p::PolyArea) = unique!(deepcopy(p))

"""
    unique!(polyarea)

Remove duplicate vertices in `polyarea`.
"""
function Base.unique!(p::PolyArea)
  close!(unique!(open!(p.outer)))
  hasholes(p) && foreach(c->close!(unique!(open!(c))), p.inners)
  p
end

"""
    bridge(polyarea)

Transform `polyarea` with holes into a single
outer chain via bridges.

## References

* Held. 1998. [FIST: Fast Industrial-Strength Triangulation of Polygons]
  (https://link.springer.com/article/10.1007/s00453-001-0028-4)
"""
function bridge(p::PolyArea)
  hasholes(p) || return first(chains(p))

  # retrieve chains with coordinates
  pchains = [coordinates.(vertices(c)) for c in chains(p)]

  # sort vertices lexicographically
  coords  = [coord for pchain in pchains for coord in pchain]
  indices = sortperm(sortperm(coords))

  # each chain has its own set of indices
  pinds = Vector{Int}[]; offset = 0
  for nvertex in length.(pchains)
    push!(pinds, indices[offset+1:offset+nvertex])
    offset += nvertex
  end

  # sort chains based on leftmost vertex
  leftmost = argmin.(pinds)
  minimums = getindex.(pinds, leftmost)
  reorder  = sortperm(minimums)
  leftmost = leftmost[reorder]
  minimums = minimums[reorder]
  pchains  = pchains[reorder]
  pinds    = pinds[reorder]

  # initialize outer boundary
  outer = first(pchains)
  oinds = first(pinds)

  # merge islands into outer boundary
  for i in 2:length(pchains)
    l = leftmost[i]
    m = minimums[i]
    c = pchains[i]
    o = pinds[i]

    # find closest vertex in boundary
    dmin, jmin = Inf, 0
    for j in findall(oinds .≤ m)
      d = sum(abs, outer[j] - c[l])
      if d < dmin
        dmin, jmin = d, j
      end
    end

    # insert island at closest vertex
    island = push!(circshift(c, -l+1), c[l])
    iinds  = push!(circshift(o, -l+1), o[l])
    outer = [outer[1:jmin]; island; outer[jmin:end]]
    oinds = [oinds[1:jmin]; iinds;  oinds[jmin:end]]
  end

  # close boundary
  push!(outer, first(outer))

  Chain(Point.(outer))
end

function Base.show(io::IO, p::PolyArea)
  outer = p.outer
  inner = isempty(p.inners) ? "" : ", "*join(p.inners, ", ")
  print(io, "PolyArea($outer$inner)")
end

function Base.show(io::IO, ::MIME"text/plain", p::PolyArea{Dim,T}) where {Dim,T}
  outer = "    └─$(p.outer)"
  inner = ["    └─$v" for v in p.inners]
  println(io, "PolyArea{$Dim,$T}")
  println(io, "  outer")
  if isempty(inner)
    print(io, outer)
  else
    println(io, outer)
    println(io, "  inner")
    print(io, join(inner, "\n"))
  end
end


text:
using Test, CubicEoS, CubicEoSDatabase

"Non-recursively compare field values of `x` and `y`."
function compare_structs(x::T, y::T) where {T}
    for field in fieldnames(T)
        @test getfield(x, field) == getfield(y, field)
    end
end

@testset "src/dbload.jl" begin
    @testset "Default databases" begin
        comp_explicit = load(
            BrusilovskyEoSComponent;
            name="methane",
            component_dbs=(Data.martinez(), Data.brusilovsky_comp())
        )
        comp_default = load(BrusilovskyEoSComponent,
            name="methane"
        )
        compare_structs(comp_explicit, comp_default)

        mix_explicit = load(
            BrusilovskyEoSMixture;
            names = ("methane", "n-pentane"),
            component_dbs=(Data.martinez(), Data.brusilovsky_comp()),
            mix_eos_db = Data.brusilovsky_mix()
        )
        mix_default = load(
            BrusilovskyEoSMixture;
            names = ("methane", "n-pentane"),
        )
        compare_structs(mix_explicit, mix_default)
    end
end

max_stars_repo_name: vvpisarev/CubicEoS

max_issues_repo_name: vvpisarev/CubicEoS

max_forks_repo_name: vvpisarev/CubicEoS


text:
using Random, StatsBase
using DataFrames
using DocStringExtensions
using CSV
include("Input.jl")
include("DEB.jl")
include("Inference.jl")

"""
Function to gather data from model object. Returns tidy DataFrame with a single row (= 1 time-step).
$(TYPEDSIGNATURES)
"""
function record_data(
        m::Model
    )

    # If there area any Fleas left, record individual-level data
    if length(m.fleas)>0
        cum_repro_test_animal = [if f.unique_id==0 f.cum_repro end for f in m.fleas]
        carapace_length_test_animal = [if f.unique_id==0 f.L/f.shape_factor end for f in m.fleas]

        ind_statevar_mat = hcat([[
                f.food_taken .* m.timestep,
                f.S_A,
                f.S_C,
                f.cum_repro,
                f.L / f.shape_factor,
                f.functresp,
                f.Xi,
                f.D[1],
                f.D[2],
                f.D[3],
                f.f_stress,
                f.g_stress,
                f.k_M_stress_1,
                f.k_M_stress_2,
                f.k_J_stress,
                f.repro_stress_1,
                f.repro_stress_2
                ] for f in m.fleas]...)

        ingestion_rate_mean,
        S_A_mean,
        S_C_mean,
        cum_repro_mean,
        carapace_length_mean,
        functresp_mean,
        Xi_mean,
        D1_mean,
        D2_mean,
        D3_mean,
        f_stress_mean,
        g_stress_mean,
        k_M_stress_1_mean,
        k_M_stress_2_mean,
        k_J_stress_mean,
        repro_stress_1_mean,
        repro_stress_2_mean =
        [robust_mean(x) for x in eachrow(ind_statevar_mat)]

        # record life stage structure
        X_juv = []
        X_ad = []
        try
            push!(X_juv, sum([f.U_H < f.U_Hp for f in m.fleas]))
        catch
            push!(X_juv, 0)
        end
        try
            push!(X_ad, sum([f.U_H >= f.U_Hp for f in m.fleas]))
        catch
            push!(X_ad, 0)
        end

        X_matrix = Array([])

        for j in 1:m.n_spec
            if j==1
                X_matrix = [get_species_abundance(m.fleas, j)]
            else
                X_matrix = hcat(X_matrix, get_species_abundance(m.fleas, j))
            end
        end

        if ndims(X_matrix)==1
            X_matrix = add_dim(X_matrix)
        end

        model_out = hcat(
                    DataFrame(
                        tday            = [m.tday],
                        F               = [m.food_density],
                        cum_repro       = robust_mean(cum_repro_test_animal),
                        carapace_length = robust_mean(carapace_length_test_animal),
                        stressor_1_conc = [m.stressor_conc[1]],
                        stressor_2_conc = [m.stressor_conc[2]],
                        stressor_3_conc = [m.stressor_conc[3]],
                        ingestion_rate_mean = ingestion_rate_mean,
                        S_A_mean = S_A_mean,
                        S_C_mean = S_C_mean,
                        cum_repro_mean       = cum_repro_mean,
                        carapace_length_mean = carapace_length_mean,
                        functresp_mean = functresp_mean,
                        Xi_mean = Xi_mean,
                        D1_mean = D1_mean,
                        D2_mean = D2_mean,
                        D3_mean = D3_mean,
                        f_stress_mean = f_stress_mean,
                        g_stress_mean = g_stress_mean,
                        k_M_stress_1  = k_M_stress_1_mean,
                        k_M_stress_2  = k_M_stress_2_mean,
                        k_J_stress    = k_J_stress_mean,
                        repro_stress_1 = repro_stress_1_mean,
                        repro_stress_2 = repro_stress_2_mean,
                        X_juv = X_juv,
                        X_ad  = X_ad,
                        X     = [sum(X_matrix)]
                    ),
                    DataFrame(X_matrix)
                )

        return model_out
    else
        model_out = hcat(
              DataFrame(
                  tday            = [m.tday],
                  F               = [m.food_density],
                  cum_repro       = Array{Union{Float64,Missing}}([missing]),
                  carapace_length = Array{Union{Float64,Missing}}([missing]),
                  stressor_1_conc = [m.stressor_conc[1]],
                  stressor_2_conc = [m.stressor_conc[2]],
                  stressor_3_conc = [m.stressor_conc[3]],
                  ingestion_rate_mean = Array{Union{Float64,Missing}}([missing]),
                  S_A_mean = Array{Union{Float64,Missing}}([missing]),
                  S_C_mean = Array{Union{Float64,Missing}}([missing]),
                  cum_repro_mean       = Array{Union{Float64,Missing}}([missing]),
                  carapace_length_mean = Array{Union{Float64,Missing}}([missing]),
                  functresp_mean = Array{Union{Float64,Missing}}([missing]),
                  Xi_mean = Array{Union{Float64,Missing}}([missing]),
                  D1_mean = Array{Union{Float64,Missing}}([missing]),
                  D2_mean = Array{Union{Float64,Missing}}([missing]),
                  D3_mean = Array{Union{Float64,Missing}}([missing]),
                  f_stress_mean = Array{Union{Float64,Missing}}([missing]),
                  g_stress_mean = Array{Union{Float64,Missing}}([missing]),
                  k_M_stress_1  = Array{Union{Float64,Missing}}([missing]),
                  k_M_stress_2  = Array{Union{Float64,Missing}}([missing]),
                  k_J_stress    = Array{Union{Float64,Missing}}([missing]),
                  repro_stress_1 = Array{Union{Float64,Missing}}([missing]),
                  repro_stress_2 = Array{Union{Float64,Missing}}([missing]),
                  X_juv = [0.],
                  X_ad  = [0.],
                  X     = [0.]
              ),
              DataFrame(hcat(repeat([0.], length(m.fleas_dict_list))'...))
          )
    end
end

"""
Function to run the DEB-IBM, given global parameters in global_dict and individual-level parameters in fleas_dict.
$(SIGNATURES)
"""
function Run(
        global_dict::Dict{String,Any},
        fleas_dict::Dict{String,Any};
        cu_ts=nothing,
        ni_ts=nothing,
        zn_ts=nothing
    )

    #### ---- Initialize arrays to record data ---- ####
    tt = []
    Nt_juv = []
    Nt_ad = []
    Ft = []
    stressor_1_conc = []
    stressor_2_conc = []
    stressor_3_conc = []
    mean_damage_1 = []

    cum_repro = []
    carapace_length = []
    ingestion_rate = []
    functresp = []
    assim_rate = []

    max_steps = (global_dict["tmax"]*global_dict["timestep"])+1

    # if no exposure data has been passed on, use constant values given in global dictionary
    if cu_ts == nothing
        cu_ts = repeat([global_dict["stressor_conc"][1]], Integer(max_steps))
    end
    if ni_ts == nothing
        ni_ts = repeat([global_dict["stressor_conc"][2]], Integer(max_steps))
    end
    if zn_ts == nothing
        zn_ts = repeat([global_dict["stressor_conc"][3]], Integer(max_steps))
    end

    # make a new model object
    m = Model(
            global_dict["remove_juveniles"],
            global_dict["renew_medium"],
            global_dict["frct_medium_renewal"],
            global_dict["timestep"],
            0, # tstep
            0, # tday
            global_dict["food_level"], # food_level
            global_dict["food_level"], # initial food_density
            global_dict["age_range_init"],
            0,
            [],
            fleas_dict,
            global_dict["stressor_conc"]
    )

    # generate initial juveniles
    for i in 1:global_dict["N0_juvenile"]
        NewFlea = make_flea(m, true, false)
        m.unique_id_count += 1
        push!(m.fleas, NewFlea)
    end
    # generate initial adults
    for i in 1:global_dict["N0_adult"]
        NewFlea = make_flea(m, false, true)
        m.unique_id_count +=1
        push!(m.fleas, NewFlea)
    end

    popsize=0

    # for every timestep
    for t in 1:max_steps
        # update time in days
        tday = (t-1)/m.timestep
        # read exposure concentrations
        m.stressor_conc[1] = cu_ts[t]
        m.stressor_conc[2] = ni_ts[t]
        m.stressor_conc[3] = zn_ts[t]
        # execute model step
        m = Step(m)
        # record population size
        popsize = length([f for f in m.fleas])

        if popsize > 5000
            @warn "Aborted simulation due to large population size"
            break
        end

        # collect time-series data
        try
            push!(Nt_juv, sum([f.U_H < f.U_Hp for f in m.fleas]))
        catch
            push!(Nt_juv, 0)
        end
        try
            push!(Nt_ad, sum([f.U_H >= f.U_Hp for f in m.fleas]))
        catch
            push!(Nt_ad, 0)
        end

        push!(Ft, m.food_density)
        push!(tt, tday)
        push!(stressor_1_conc, m.stressor_conc[1])
        push!(stressor_2_conc, m.stressor_conc[2])
        push!(stressor_3_conc, m.stressor_conc[3])
        try
            push!(cum_repro, mean([f.cum_repro for f in m.fleas]))
            push!(carapace_length, [if f.unique_id==0 f.L/f.shape_factor end for f in m.fleas][1])
            push!(ingestion_rate, mean([f.food_taken for f in m.fleas])*m.timestep)
            push!(functresp, mean(skipmissing([f.functresp for f in m.fleas])))
            push!(assim_rate, mean(skipmissing([f.S_A*m.timestep for f in m.fleas])))
        catch
            push!(cum_repro, missing)
            push!(carapace_length, missing)
            push!(ingestion_rate, missing)
            push!(functresp, missing)
            push!(assim_rate, missing)
        end
    end

    model_out = DataFrame(
            tday = tt, # time in days
            N_juv = Nt_juv, # juveniles
            N_ad  = Nt_ad, # adults
            N = Nt_juv .+ Nt_ad, # population
            F = Ft, # food density
            stressor_1_conc = stressor_1_conc,
            stressor_2_conc = stressor_2_conc,
            stressor_3_conc = stressor_3_conc,
            cum_repro = cum_repro,
            carapace_length = carapace_length,
            ingestion_rate = ingestion_rate,
            functresp = functresp,
            assim_rate = assim_rate
        )

    model_out[isnothing.(model_out.carapace_length),:carapace_length] .= missing
    model_out[isnothing.(model_out.cum_repro),:cum_repro] .= missing

    # return time-series data as dataframe
    return model_out
end

"""
Function to simulate single-stressor toxicity tests. \n
$(SIGNATURES)
"""
function SingleStressorTest(
        global_dict::Dict{String, Any},
        fleas_dict::Dict{String, Any},
        stressor::Symbol,
        exposure;
        accepted_DEB = nothing,
        params_DEB = ["p_Am", "F_m", "y_EX", "mortality_constant"],
        accepted_stressor = nothing,
        replicates = 10
    )

    df = DataFrame()

    for rep in 1:replicates
        #### ---- Obtain parameter values ---- ####

        # make sure to always start with default parameters
        fleas_dict = copy(fleas_dict)

        # IF KDEs havee been passed on, sample parameter values from KDE
        # otherwise, use point estimates from fleas_dict

        # get DEB parameters
        if accepted_DEB != nothing
            sample = posterior_sample(accepted_DEB)
            for (i,par) in enumerate(params_DEB)
                    fleas_dict[par] = sample[i]
            end
            fleas_dict = fleas_dict |> convert_parameters |> estimate_egg_weight
        end

        # get toxicity parameters
        if accepted_stressor != nothing
            par_sample = posterior_sample(accepted_stressor)
            if stressor in [:Cu_ts, :Cu]
                fleas_dict = assign_params(fleas_dict, par_sample, string(stressor[1:2]))
            elseif stressor in [:Ni_ts, :Ni]
                fleas_dict = assign_params(fleas_dict, par_sample, string(stressor[1:2]))
            elseif stressor in [:Zn_ts, :Zn]
                fleas_dict = assign_params(fleas_dict, par_sample, string(stressor[1:2]))
            else
                error("Unknown stressor")
            end
        end

        # update compound parameters
        fleas_dict = fleas_dict |> convert_parameters |> estimate_egg_weight

        #### ---- Run the simulations ---- ####
        m = DataFrame()
        if typeof(exposure)==DataFrame
            for trt in unique(exposure.treatment)
                # get exposure time series
                stressor_ts = exposure[exposure.treatment.==trt,stressor]
                # run the model
                if stressor == :Cu_ts
                    m_trt = Run(global_dict, fleas_dict; cu_ts=stressor_ts)
                    m_trt.treatment = trt
                    m_trt.replicate = rep
                    append!(m, m_trt)
                elseif stressor == :Ni_ts
                    m_trt = Run(global_dict, fleas_dict; ni_ts=stressor_ts)
                    m_trt.treatment = trt
                    m_trt.replicate = rep
                    append!(m, m_trt)
                elseif stressor == :Zn_ts
                    m_trt = Run(global_dict, fleas_dict; zn_ts=stressor_ts)
                    m_trt.treatment = trt
                    m_trt.replicate = rep
                    append!(m, m_trt)
                else
                    error("Stressor has to be one of :Cu_ts, :Ni_ts, Zn_ts")
                end
            end

        elseif typeof(exposure)==Array{Float64,1}
            treatment_names = string(stressor) .* string.([x[1]-1 for x in collect(enumerate(exposure))])
            treatment_names[1] = "Co"
            # for every treatment
            for (i,trt) in enumerate(exposure)
                # run simulation with constant exposure
                if stressor == :Cu
                    global_dict["stressor_conc"][1] = trt
                    m_trt = Run(global_dict, fleas_dict)
                    m_trt[!,:treatment] .= treatment_names[i]
                    m_trt[!,:replicate] .= rep
                    append!(m, m_trt)
                elseif stressor == :Ni
                    global_dict["stressor_conc"][2] = trt
                    m_trt = Run(global_dict, fleas_dict)
                    m_trt[!,:treatment] .= treatment_names[i]
                    m_trt[!,:replicate] .= rep
                    append!(m, m_trt)
                elseif stressor == :Zn
                    global_dict["stressor_conc"][3] = trt
                    m_trt = Run(global_dict, fleas_dict)
                    m_trt[!,:treatment] .= treatment_names[i]
                    m_trt[!,:replicate] .= rep
                    append!(m, m_trt)
                else
                    error("Stressor has to be one of :Cu, :Ni, :Zn")
                end
            end
        else
            error("Exposure has to be of type DataFrame or Array{Float64, 1}, not "*string(typeof(exposure)))
        end
        # add a column for relative responses
        ref = by(m[m.treatment.=="Co",:], [:tday]) do df
            DataFrame(ref=mean(df.N))
        end

        m = join(m, ref, on=:tday)
        m.response = m.N ./ m.ref
        select!(m, Not(:ref))

        append!(df, m)

    end
    return df
end

"""
Function to simulate mixture treatments.
$(TYPEDSIGNATURES)
"""
function MixtureTest(
        global_dict,
        fleas_dict,
        exposure;
        accepted_DEB = nothing,
        params_DEB = ["p_Am", "F_m", "y_EX", "mortality_constant"],
        accepted_stressor_array = nothing,
        replicates = 10,
        treatments=["Co", "Mix1", "Mix2", "Mix3", "Mix4", "Mix5"]
    )

    df = DataFrame()
    for rep in 1:replicates
        m = DataFrame()

        # sample DEB parameters
        if accepted_DEB != nothing
            sample = posterior_sample(accepted_DEB)
            for (i,par) in enumerate(params_DEB)
                    fleas_dict[par] = sample[i]
            end
            fleas_dict = fleas_dict |> convert_parameters |> estimate_egg_weight
        end

        # sample toxicity parameters
        if accepted_stressor_array!=nothing
            for (i,accepted_stressor) in enumerate(accepted_stressor_array)
                particle = posterior_sample(accepted_stressor)
                fleas_dict = assign_params(fleas_dict, particle, i)
            end
        end

        if typeof(exposure)==DataFrame
            for trt in unique(exposure.treatment)
                exp_trt = exposure[exposure.treatment.==trt,:]
                m_trt = Run(global_dict, fleas_dict; cu_ts=exp_trt[:,:Cu_ts], ni_ts=exp_trt[:,:Ni_ts], zn_ts=exp_trt[:,:Zn_ts])
                m_trt.treatment = trt
                append!(m, m_trt)
            end
        else
            for (i,trt) in enumerate(eachrow(exposure))
                global_dict["stressor_conc"][1] = trt[1]
                global_dict["stressor_conc"][2] = trt[2]
                global_dict["stressor_conc"][3] = trt[3]
                m_trt = Run(global_dict, fleas_dict)
                m_trt[!,:treatment] .= treatments[i]
                append!(m, m_trt)
            end
        end

        # Finally, add a column for relative responses
        ref = by(m[m.treatment.=="Co",:], [:tday]) do df
            DataFrame(ref=mean(df.N))
        end

        m = join(m, ref, on=:tday)
        m.response = m.N ./ m.ref
        select!(m, Not(:ref))
        m[!,:replicate] .= rep
        append!(df, m)
    end

    return df
end


text:
# Quick notes for tests. None is actually implemented. This file is not even supposed to work.

using FactCheck

facts("Module: plot") do
  context() do
    # a: no sc, i.e. only one scenario per reservoir.     => with sc, must throw an error.
    # d: no sc, check the second output argument is empty.
    a = getDataFrameForRiver([VesdreReservoir_adam2015_95_avg, VesdreReservoir_adam2015_95_low, VesdreReservoir_adam2015_95_high], 1:ruleLength, "Vesdre", multipleScenarios=:None)
    b = getDataFrameForRiver([sc, VesdreReservoir_adam2015_95_avg, VesdreReservoir_adam2015_95_low, VesdreReservoir_adam2015_95_high], 1:ruleLength, "Vesdre", multipleScenarios=:Merge)
    c = getDataFrameForRiver([sc, VesdreReservoir_adam2015_95_avg, VesdreReservoir_adam2015_95_low, VesdreReservoir_adam2015_95_high], 1:ruleLength, "Vesdre", multipleScenarios=:Split)
    d = getDataFrameForRiver([VesdreReservoir_adam2015_95_avg, VesdreReservoir_adam2015_95_low, VesdreReservoir_adam2015_95_high], 1:ruleLength, "Vesdre", multipleScenarios=:Split)
  end
end

facts("Module: read") do
  context() do
    b = makeScenarios(Vesdre_ts, sampleDuration=Week(1))
    c = makeScenarios(Vesdre_ts, sampleDuration=Day(1), resamplingRate=28) ---> throws()
  end
end

facts("Module: reservoir") do
  context("accessors: statistical") do
    splitIntoSets(Vesdre.scenarios, [1])
    splitIntoSets(Vesdre.scenarios, [.8, .2])
  end

  context("accessors") do
    totalInflow(VesdreReservoir, 1)
    totalInflow(VesdreReservoir, 1:22)
    totalInflow(VesdreReservoir, 1, asTimeSeries=true)
    totalInflow(VesdreReservoir, 1:22, asTimeSeries=true)
  end
end


text:
using CSV, DataFrames, DSP, HDF5, Random, Statistics

function fp32int8(A::AbstractArray)
    A .-= mean(A,dims=1)
    maxs = maximum(abs,A,dims=1)
    factors = 127 ./ maxs
    A .*= factors
    return round.(Int8,A)
end

function fp32uint8(A::AbstractArray)
    m = min(minimum(A),0)
    factor = 255 / (maximum(A) - m)
    return round.(UInt8,(A .- m) .* factor)
end

function vel2acc(A::AbstractArray,freqs::AbstractArray)
    Xw = fft(A,1)
    Xw .*= freqs .* 1im .* 2π
    return real.(ifft(Xw,1))
end

function taper!(A::AbstractArray;α::Real=0.05)
    window = tukey(size(A,1),α)
    A .*= window
    return nothing
end

function randamp!(A::AbstractArray)

    return nothing
end

function wav2spec(A::AbstractArray,n,overlap,fs)
    spec = spectrogram(A[:,1],n,overlap,fs=fs)
    Nrows, Ncols = size(spec.power)
    out = zeros(eltype(A),Nrows,Ncols,3)
    out[:,:,1] .= spec.power
    for ii = 2:3
        out[:,:,ii] .= spectrogram(A[:,ii],n,overlap,fs=fs).power
    end

    # make smallest value across channels equal to 1
    out .+= (1 - minimum(out))\

    # return log10 of power
    return log10.(out)
end

function ricker(f,n,dt,t0)
    # Create the wavelet and shift in time if needed
    T = dt*(n-1)
    t = 0:dt:T
    tau = t .- t0
    s = (1 .-tau .* tau .* f^2 * pi^2) .* exp.(-tau .^ 2 .* pi^2 .*f^2)
    return s
end

function specplot(A::AbstractArray;n=50,overlap=45,fs=100.)
    spec = spectrogram(A,n,overlap,fs=fs)
    t = 0 : 1 / fs : size(A,1) / fs - 1 / fs
    p1 = plot(t,A)
    xlims!(p1,minimum(spec.time),maximum(spec.time))
    p2 = heatmap(
        spec.time,
        spec.freq,
        10 .* log10.(spec.power),
        legend=:none,
        xlabel="Time [s]",
        ylabel="Frequency [Hz]",
    )
    l = @layout [a;b]
    p = plot(p1,p2, layout = l)
    display(p)
    return nothing
end

function highpass!(A::AbstractArray{<:AbstractFloat},freq::Real,fs::Real)
    T = eltype(A)
    # create filter
    responsetype = Highpass(T(freq); fs=fs)
    designmethod = Butterworth(T,4)
    A[:,:] .= filt(digitalfilter(responsetype, designmethod), @view(A[:,:]))
    return nothing
end

function detrend!(X::AbstractArray{<:AbstractFloat})
    T = eltype(X)
    N = size(X,1)

    # create linear trend matrix
    A = similar(X,T,N,2)
    A[:,2] .= T(1)
    A[:,1] .= range(T(0),T(1),length=N)
    # create linear trend matrix
    R = transpose(A) * A

    # do the matrix inverse for 2x2 matrix
    # this is really slow on GPU
    Rinv = inv(R)
    factor = Rinv * transpose(A)

    # remove trend
    X .-= A * (factor * X)
    return nothing
end

function builddirs(DATADIR)
    WAVES = joinpath(DATADIR,"waves")
    WAVETRAIN = joinpath(WAVES,"train")
    WAVETEST = joinpath(WAVES,"test")
    SPECS = joinpath(DATADIR,"specs")
    SPECTRAIN = joinpath(SPECS,"train")
    SPECTEST = joinpath(SPECS,"test")
    for dir in [WAVES,WAVETRAIN,WAVETEST,SPECS,SPECTRAIN,SPECTEST]
        if !isdir(dir)
            mkpath(dir)
        end
    end

    # make directories for each class
    for dir in [WAVETRAIN,WAVETEST,SPECTRAIN,SPECTEST]
        for class in ["noise","p-wave","s-wave"]
            tmpdir = joinpath(dir,class)
            if !isdir(tmpdir)
                mkpath(tmpdir)
            end
        end
    end
    return WAVES, WAVETRAIN, WAVETEST, SPECS, SPECTRAIN, SPECTEST
end

function test_train_split(df::DataFrame)
    ind = findall(df[:,:receiver_type] .∈ Ref(["HH","EH"]))
    Hdf = df[ind,:]

    # split into test and train datasets
    trainind = findall((Hdf[:,:receiver_latitude] .> 30) .& (Hdf[:,:receiver_latitude] .< 50) .& (Hdf[:,:receiver_longitude] .> -130) .& (Hdf[:,:receiver_longitude] .< -108))
    testind = collect(1:size(Hdf,1))
    deleteat!(testind,trainind)

    # filter out data in W Canada and Texas
    CDind = findall((Hdf[:,:receiver_latitude] .> 50) .& (Hdf[:,:receiver_latitude] .< 55.5) .& (Hdf[:,:receiver_longitude] .> -130.5) .& (Hdf[:,:receiver_longitude] .< -112))
    TXind = findall((Hdf[:,:receiver_latitude] .> 30) .& (Hdf[:,:receiver_latitude] .< 50) .& (Hdf[:,:receiver_longitude] .> -109) .& (Hdf[:,:receiver_longitude] .< -101.5))
    filter!(x-> x ∉ CDind,testind)
    filter!(x-> x ∉ TXind,testind)

    # separate into NOISE and EQ datasets for test / train
    EQtrainind = findall(.!ismissing.(Hdf[trainind,:p_status]))
    EQtestind = findall(.!ismissing.(Hdf[testind,:p_status]))
    noisetrainind = findall(ismissing.(Hdf[trainind,:p_status]))
    noisetestind = findall(ismissing.(Hdf[testind,:p_status]))
    EQtrain = Hdf[trainind[EQtrainind],:]
    EQtest = Hdf[testind[EQtestind],:]
    NOISEtrain = Hdf[trainind[noisetrainind],:]
    NOISEtest = Hdf[testind[noisetestind],:]

    # take 80 / 20 split for test / train EQ
    NEQtest = size(EQtest,1)
    NEQtrain = NEQtest * 4
    train80 = randperm(size(EQtrain,1))[1:NEQtrain]
    EQtrain = EQtrain[train80,:]
    EQtest = EQtest[randperm(NEQtest)[1:NEQtest],:]

    # even out number of stations for NOISE datasets
    NNtest = size(NOISEtest,1)
    NNtrain = size(NOISEtrain,1)
    NN = min(NNtest,NNtrain)
    NOISEtest = NOISEtest[randperm(NNtest)[1:NN],:]
    NOISEtrain = NOISEtrain[randperm(NNtrain)[1:NN],:]
    return EQtrain, EQtest, NOISEtrain, NOISEtest
end

function augment_noise!(A::AbstractArray,spikes,ramp,)
end

function preprocess_noise(
    df::DataFrame,
    datafile,
    DATADIR,
    nwindows;
    train=true,
    α=0.065,
    freq=2.,
    window_length=2.,
    n = 50,
    overlap=45,
    fs=100.,
)

    # load (meta)data
    h5file = h5open(datafile,"r")
    N = size(df,1)
    if train
        SPEC = joinpath(DATADIR,"specs","train","noise")
        WAVE = joinpath(DATADIR,"waves","train","noise")
    else
        SPEC = joinpath(DATADIR,"specs","test","noise")
        WAVE = joinpath(DATADIR,"waves","test","noise")
    end

    # start after taper
    window_samples = round(Int,window_length * 100)
    total_windows = round(Int,6000 / window_samples)
    startwindow = round(Int,6000 * α / window_samples + 1)
    for ii = 1:N
        path = df[ii,:trace_name]
        println("Reading file $path $ii of $N")
        data = read(h5file,"data/" * path) |> transpose |> Array
        detrend!(data)
        taper!(data,α=α)
        highpass!(data,freq,100.)
        for ii = startwindow:min(total_windows - startwindow,startwindow+nwindows-1)
            wave = data[ii*window_samples:(ii+1)*window_samples-1,:]
            qwave = fp32uint8(wave)
            spec = wav2spec(wave,n,overlap,fs)
            qspec = fp32uint8(spec)
            waveout = joinpath(WAVE,"$(path)_$ii.jpg")
            specout = joinpath(SPEC,"$(path)_$ii.jpg")
            save(waveout,qwave)
            save(specout,qspec)
        end
    end
    close(h5file)
    return nothing
end

function preprocess_waves(df::DataFrame,datafile,DIR)
    # load (meta)data
    h5file = h5open(datafile,"r")
    N = size(df,1)
    PWAVE = joinpath(DIR,"pwave")
    SWAVE = joinpath(DIR,"swave")
    for ii = 1:N
        path = df[ii,:trace_name]
        data = read(file,"data/" * path) |> transpose |> Array
        detrend!(data)
        taper!(data,α=0.1)
        highpass!(data,2.,100.)
    end
end

### training workflow
# 1. read from HDF5
# 2. detrend
# 3. taper
# 4. highpass filter above 1 Hz
# 4. convert to acceleration
# 5. window earthquakes around arrival randomly
# 5a. add gaussian normal noise
# 5b. clip amplitudes randomly
# 6. convert to spectrogram or not
# 6a. perform spectral averaging
# 7. quantization
# 8. save into train / test

# directory structure
BASEDIR = expanduser("~/CS249FINAL")
DATADIR = joinpath(BASEDIR,"data")
WAVES, WAVETRAIN, WAVETEST, SPECS, SPECTRAIN, SPECTEST = builddirs(DATADIR)

# load metadata
df = CSV.File(joinpath(DATADIR,"merge.csv"),"r") |> DataFrame
EQtrain, EQtest, NOISEtrain, NOISEtest = test_train_split(df)

# need 22 two-second windows per training noise sample
# need 6 two-second windows per test noise sample


text:
abstract type AbstractRoute end
struct StaticRoute <: AbstractRoute
    url ::String
    StaticRoute(url) = new(url)
end
struct DynamicRoute{ST,VT} <: AbstractRoute
    segments_length ::Int
    static_segments ::ST
    variables ::VT
    tailed::Bool

    DynamicRoute(segments_length, static_segments::ST, variables::VT, tailed) where {ST, VT} = new{ST, VT}(segments_length, static_segments, variables, tailed)
end

function DynamicRoute(url::AbstractString)
    parts = split(url, '/', keepempty = false)
    segments_vector = Tuple{Int, String}[]
    variables = Dict{Symbol, Int}()
    segments_count = 0
    for  part in parts
        segments_count += 1
        if length(part) > 2 && part[1] == '<' && part[end] == '>'
            var_name = Symbol(part[2:end-1])
            var_name in keys(variables) && ArgumentError("duplicated variable name $(var_name)")
            variables[var_name] = segments_count
        elseif part != "*"
            push!(segments_vector,(segments_count, String(part)))
        end
    end
    #Route is tailed if it url hasn't / at end and the last segment is dynamic
    tailed = false
    if url[end] != '/' && (isempty(segments_vector) || segments_vector[end][1] != length(parts))
        tailed = true
    end
    return DynamicRoute(
        segments_count,
        (segments_vector...,),
        (;variables...),
        tailed
    )
end

struct RouteHandler{RT <: AbstractRoute, FT <: Function}
    route::RT
    handler::FT
    RouteHandler(route::RT, handler::FT) where {RT <: AbstractRoute, FT <: Function} = new{RT, FT}(route, handler)
end

function try_handle(route_handler::RouteHandler{StaticRoute, FT}, path::AbstractString, request::HTTP.Request, args...) where {FT}
    route_handler.route.url != path && return nothing
    return route_handler.handler(request, args...)
end

function args_tuple(route::DynamicRoute{ST, NamedTuple{Names, T}}, parts) where {ST, Names, T}
    return NamedTuple{Names}(getindex.(Ref(parts), values(route.variables)))
end
function try_handle(route_handler::RouteHandler{<:DynamicRoute, FT}, path::AbstractString, request::HTTP.Request, args...) where {FT}
    route = route_handler.route
    parts_length = route.segments_length
    limit = route.tailed ? parts_length : 0
    parts = split(lstrip(path, '/'), '/', limit = limit, keepempty = false)
    length(parts) != parts_length && return nothing
    for segment in route.static_segments
        parts[segment[1]] != segment[2] && return nothing
    end
    return route_handler.handler(request, args...;args_tuple(route, parts)...)
end

function _make_route(url::AbstractString)
    parts = split(url, '/', keepempty = false)
    for part in parts
        if part == "*" || (length(part) > 2 && part[1] == '<' && part[end] == '>')
            return DynamicRoute(url)
        end
    end
    return StaticRoute(url)
end

struct Route{RH}
    method ::Union{Nothing,String}
    route_handler ::RH
    Route(method, route_handler::RH) where {RH} = new{RH}(method, route_handler)
end
function Route(handler::Function, method, path::AbstractString)
    return Route(
        method,
        RouteHandler(
            _make_route(path),
            handler
        )
    )
end
Route(handler::Function, path::AbstractString) = Route(handler, nothing, path)

function try_handle(route::Route, path::AbstractString, request::HTTP.Request, args...)
    (!isnothing(route.method) && route.method != request.method) && return nothing
    return try_handle(route.route_handler, path, request, args...)
end

@inline function _handle(route_tuple::Tuple, path::AbstractString, request::HTTP.Request, args...)
    res = try_handle(route_tuple[1], path, request, args...)
    return !isnothing(res) ?
            res :
            _handle(Base.tail(route_tuple), path, request, args...)
end

@inline function _handle(route_tuple::Tuple{}, path::AbstractString, request::HTTP.Request, args...)
    return HTTP.Response(404)
end

function handle(route_tuple::Tuple, path::AbstractString, request::HTTP.Request, args...)::HTTP.Response
    return _handle(route_tuple, path, request, args...)
end

mutable struct Router
    routes::Tuple
    Router() = new(())
    Router(routes::Tuple{Vararg{Route}}) = new(routes)
end

function add_route!(router::Router, route::Route)
    router.routes = (router.routes..., route)
end
add_route!(handler::Function, router::Router, method, url) = add_route!(router, Route(handler, method, url))
add_route!(handler::Function, router::Router, url) = add_route!(handler, router, nothing, url)

function HTTP.handle(router::Router, request::HTTP.Request, args...)
    path = HTTP.URI(request.target).path
    return handle(router.routes, path, request, args...)
end


text:
using Revise
Revise.revise()

using Term
import Term: install_stacktrace

install_stacktrace()
# -------------------------------- nested tags ------------------------------- #
println(
    RenderableText(
        "You [green]can nest [blue underline] one [on_gold3 black bold]style[/on_gold3 black bold] inside  [/blue underline] another [/green]too!",
    ),
)

# ----------------------------- >2 tags per line ----------------------------- #
println(
    RenderableText(
        "[red]reeed [blue]blueeee[/blue] reeeed [green] greeen [/green] red red red[/red]",
    ),
)

# ------------------------------ text reshaping ------------------------------ #
text =
    "This is a [blue][gold3]very[/gold3] long[/blue] piece of [green]text with [red]nested[/red] styles[/green]!! "^6

for width in [9, 29, 37, 50]
    print("\n\n")
    println("Width: $width")
    println("."^width)

    print(RenderableText(text; width = width))
end

text = "[red]TEST[/red]aaa"^20
# text = "TESTaaa"^10

for width in [21, 44, 81]
    print("\n\n")
    println("Width: $width")
    println("."^width)

    print(RenderableText(text; width = width))
end


text:
using Sockets

function readview(sock)
	view = readline(sock)
	view == "" && return
	size = length(view) - 1
	view *= "\n"
	for i = 1:size
		line = readline(sock) * "\n"
		view = view * line
	end
	return view
end

sock = connect(get(ARGS, 1, "localhost"), get(ARGS, 2, 63187))
while isopen(sock)
	view = readview(sock)
	view == "" && break
	print(view, "Command (q<>^v): ")
	cmd = read(stdin, Char)
	print("\n")
	if cmd == 'q'
		break
	elseif cmd == '\n'
		cmd = '^'
	end
	write(sock, cmd)
end
close(sock)

max_stars_repo_name: markusfisch/terrain-parsing-game

max_issues_repo_name: markusfisch/terrain-parsing-game

max_forks_repo_name: markusfisch/terrain-parsing-game


text:
x = 3
x = x+1
println("hi")

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
using ProgressiveAligner
using Base.Test

# write your own tests here
@test 1 == 1
tests = [
    "DataReaderTests",
    #"DataWriterTests",
    "ClusteringTests",
    "ProfileAlignerTests"
    ]

println("Running tests:")
for t in tests
    println(" * $(t)")
    include("$(t).jl")
end