text:
"UnorderedOrdering - Indicates that a collection's is not forward or reverse ordered."
struct UnorderedOrdering <: Ordering end
const Unordered = UnorderedOrdering()

"UnkownOrdering - Indicates that a collection's ordering is not known."
struct UnkownOrdering <: Ordering end
const UnkownOrder = UnkownOrdering()

"""
    is_forward(x) -> Bool

Returns `true` if `x` is sorted forward.
"""
is_forward(x) = issorted(x)
is_forward(::ForwardOrdering) = true
is_forward(::Ordering) = false
is_forward(::AbstractUnitRange) = true
is_forward(x::AbstractRange) = step(x) > 0

"""
    is_reverse(x) -> Bool

Returns `true` if `x` is sorted in reverse.
"""
is_reverse(x) = issorted(x, order=Reverse)
is_reverse(::ReverseOrdering) = true
is_reverse(::Ordering) = false
is_reverse(::AbstractUnitRange) = false
is_reverse(x::AbstractRange) = step(x) < 0

"""
    order(x) -> Ordering

Returns the ordering of `x`.
"""
order(x::T) where {T} = _order(order(T), x)
function _order(::UnkownOrdering, x)
    if is_reverse(x)
        return Reverse
    elseif is_forward(x)
        return Forward
    else
        return Unordered
    end
end
_order(xo::Ordering, x) = xo
order(::Type{T}) where {T} = UnkownOrder
order(::Type{T}) where {T<:AbstractUnitRange} = Forward

"""
    is_ordered(x) -> Bool

Returns `true` if `x` is ordered. `is_ordered` should return the same value that
`issorted` would on `x` except it doesn't specify how it's sorted (e.g.,
forward, reverse, etc).
"""
is_ordered(::Type{T}) where {T} = false
is_ordered(::Type{T}) where {T<:AbstractRange} = true
is_ordered(x::X) where {X} = is_ordered(X) ? true : is_forward(x) || is_reverse(x)

"""
    ordmax(x) = ordmax(x, order(x))
    ordmax(x::T, ::Ordering) -> T

Finds the maximum of `x` using information about its ordering.
"""
ordmax(x) = ordmax(x, order(x))
ordmax(x, ::ForwardOrdering) = last(x)
ordmax(x, ::ReverseOrdering) = first(x)
ordmax(x, ::UnorderedOrdering) = maximum(x)

"""
    ordmin(x) = ordmin(x, order(x))
    ordmin(x::T, ::Ordering) -> T

Finds the minimum of `x` using information about its ordering.
"""
ordmin(x) = ordmin(x, order(x))
ordmin(x, ::ForwardOrdering) = first(x)
ordmin(x, ::ReverseOrdering) = last(x)
ordmin(x, ::UnorderedOrdering) = minimum(x)

"""
    find_max(x)

Returns the index of the maximum value for `x`. Differes from `findmax` by
accounting for any sorting.
"""
find_max(x) = find_max(x, order(x))
find_max(x, ::ForwardOrdering) = (last(x), lastindex(x))
find_max(x, ::ReverseOrdering) = (first(x), firstindex(x))
find_max(x, ::UnorderedOrdering) = findmax(x)

"""
    find_min(x)

Returns the index of the minimum value for `x`. Differes from `findmin` by
accounting for any sorting.
"""
find_min(x) = find_min(x, order(x))
find_min(x, ::ForwardOrdering) = (first(x), firstindex(x))
find_min(x, ::ReverseOrdering) = (last(x), lastindex(x))
find_min(x, ::UnorderedOrdering) = findmin(x)

"""
    is_within(x, y) -> Bool

Returns `true` if all of `x` is found within `y`.
"""
is_within(x, y) = is_within(x, order(x), y, order(y))
is_within(x, xo, y, yo) = (ordmin(x, xo) >= ordmin(y, yo)) && (ordmax(x, xo) <= ordmax(y, yo))

"""
    gtmax(x, y) -> Bool

Returns `true` if the maximum of `x` is greater than that of `y`.
"""
gtmax(x, y) = gtmax(x, order(x), y, order(y))
gtmax(x, xo, y, yo) = ordmax(x, xo) > ordmax(y, yo)

"""
    ltmax(x, y) -> Bool

Returns `true` if the maximum of `x` is less than that of `y`.
"""
ltmax(x, y) = ltmax(x, order(x), y, order(y))
ltmax(x, xo, y, yo) = ordmax(x, xo) < ordmax(y, yo)

"""
    eqmax(x, y) -> Bool

Returns `true` if the maximum of `x` and `y` are equal.
"""
eqmax(x, y) = eqmax(x, order(x), y, order(y))
eqmax(x, xo, y, yo) = ordmax(x, xo) == ordmax(y, yo)

"""
    gtmin(x, y) -> Bool

Returns `true` if the minimum of `x` is greater than that of `y`.
"""
gtmin(x, y) = gtmin(x, order(x), y, order(y))
gtmin(x, xo, y, yo) = ordmin(x, xo) > ordmin(y, yo)

"""
    ltmin(x, y) -> Bool

Returns `true` if the minimum of `x` is less than that of `y`.
"""
ltmin(x, y) = ltmin(x, order(x), y, order(y))
ltmin(x, xo, y, yo) = ordmin(x, xo) < ordmin(y, yo)

"""
    eqmin(x, y) -> Bool

Returns `true` if the minimum of `x` and `y` are equal.
"""
eqmin(x, y) = eqmin(x, order(x), y, order(y))
eqmin(x, xo, y, yo) = ordmin(x, xo) == ordmin(y, yo)

"""
    group_max(x, y[, z...])

Returns the maximum value of all collctions.
"""
group_max(x, y, z...) = max(group_max(x, y), group_max(z...))
group_max(x) = ordmax(x)
group_max(x, y) = _group_max(x, order(x), y, order(y))
_group_max(x, xo, y, yo) = max(ordmax(x, xo), ordmax(y, yo))

"""
    group_min(x, y[, z...])

Returns the minimum value of all collctions.
"""
group_min(x, y, z...) = min(group_min(x, y), group_min(z...))
group_min(x) = ordmin(x)
group_min(x, y) = _group_min(x, order(x), y, order(y))
_group_min(x, xo, y, yo) = min(ordmin(x, xo), ordmin(y, yo))

"""
    cmpmax(x, y)
"""
cmpmax(x, y) = cmpmax(x, order(x), y, order(y))
cmpmax(x, xo, y, yo) = ltmax(x, xo, y, yo) ? -1 : (gtmax(x, xo, y, yo) ? 1 : 0)

"""
    cmpmin(x, y)
"""
cmpmin(x, y) = cmpmin(x, order(x), y, order(y))
cmpmin(x, xo, y, yo) = ltmin(x, xo, y, yo) ? -1 : (gtmin(x, xo, y, yo) ? 1 : 0)

"""
    min_of_group_max(x, y)

Returns the minimum of maximum of `x` and `y`. Functionally equivalent to
`min(maximum(x), maximum(y))` but uses trait information about ordering for
improved performance.
"""
min_of_group_max(x, y) = min_of_group_max(x, order(x), y, order(y))
min_of_group_max(x, xo, y, yo) = min(ordmax(x, xo), ordmax(y, yo))

"""
    max_of_group_min(x, y)

Returns the maximum of minimum of `x` and `y`. Functionally equivalent to
`max(minimum(x), minimum(y))` but uses trait information about ordering for
improved performance.
"""
max_of_group_min(x, y) = max_of_group_min(x, order(x), y, order(y))
max_of_group_min(x, xo, y, yo) = max(ordmin(x, xo), ordmin(y, yo))

"""
    is_before(x::T, y::T, collection::AbstractVector{T}) -> Bool

Returns `true` if `x` is before `y` in `collection`.
"""
function is_before(x, y, collection)
    findfirst(isequal(x), collection) < findfirst(isequal(y), collection)
end

"""
    is_before(x::AbstractVector{T}, y::AbstractVector{T}) -> is_before(order(x), order(y), x, y)
    is_before(::Ordering, ::Ordering, x, y) -> Bool

Returns `true` if all elements in `x` are before all elements in `y`. Functionally
equivalent to `all(x .< y)`.
"""
is_before(x, y) = is_before(x, order(x), y, order(y))
is_before(x, xo, y, yo) = ordmax(x, xo) < ordmin(y, yo)

"""
    is_after(x::T, y::T, collection::AbstractVector{T}) -> Bool

Returns `true` if `x` is after `y` in `collection`.
"""
function is_after(x, y, collection)
    findfirst(isequal(x), collection) > findfirst(isequal(y), collection)
end

"""
    is_after(x::AbstractVector{T}, y::AbstractVector{T}) -> is_after(order(x), order(y), x, y)
    is_after(::Ordering, ::Ordering, x, y) -> Bool

Returns `true` if all elements in `x` are after all elements in `y`. Functionally
equivalent to `all(x .> y)`.
"""
is_after(x, y) = is_after(x, order(x), y, order(y))
is_after(x, xo, y, yo) = ordmin(x, xo) > ordmax(y, yo)

"""
    is_contiguous(x, y) = is_contiguous(order(x), order(y), x, y)
    is_contiguous(::Ordering, ::Ordering, x, y) -> Bool

Returns `true` if one of the ends of `x` may be extended by a single overlapping
end of `y`.

# Example
```
julia> is_contiguous(1:3, 3:4) == true

julia> is_contiguous(3:-1:1, 3:4) == true

julia> is_contiguous(3:-1:1, 4:-1:3) == true

julia> is_contiguous(1:3, 4:-1:3) == true

julia> is_contiguous(1:3, 2:4) == false
```
"""
is_contiguous(x, y) = is_contiguous(x, order(x), y, order(y))
function is_contiguous(x, ::ForwardOrdering, y, yo)
    return last(x) == ordmin(y, yo) || first(x) == ordmax(y, yo)
end
function is_contiguous(x, ::ReverseOrdering, y, yo)
    return last(x) == ordmax(y, yo) || first(x) == ordmin(y, yo)
end

"""
    next_type(x::T)

Returns the immediately greater value of type `T`.
"""
function next_type(x::AbstractString)
    isempty(x) && return ""
    return x[1:prevind(x, lastindex(x))] * (last(x) + 1)
end
next_type(x::Symbol) = Symbol(next_type(string(x)))
next_type(x::AbstractChar) = x + 1
next_type(x::T) where {T<:AbstractFloat} = nextfloat(x)
next_type(x::T) where {T} = x + one(T)

"""
    prev_type(x::T)

Returns the immediately lesser value of type `T`.
"""
function prev_type(x::AbstractString)
    isempty(x) && return ""
    return x[1:prevind(x, lastindex(x))] * (last(x) - 1)
end
prev_type(x::Symbol) = Symbol(prev_type(string(x)))
prev_type(x::AbstractChar) = x - 1
prev_type(x::T) where {T<:AbstractFloat} = prevfloat(x)
prev_type(x::T) where {T} = x - one(T)


text:
@testset "Noop Codec" begin
    source = IOBuffer("")
    stream = TranscodingStream(Noop(), source)
    @test eof(stream)
    @inferred eof(stream)
    @test read(stream) == UInt8[]
    @test occursin("mode=read", repr(stream))

    source = IOBuffer("foo")
    stream = TranscodingStream(Noop(), source)
    @test !eof(stream)
    @test read(stream) == b"foo"
    close(stream)

    data = rand(UInt8, 100_000)
    source = IOBuffer(data)
    stream = TranscodingStream(Noop(), source)
    @test !eof(stream)
    @test read(stream) == data
    close(stream)

    stream = TranscodingStream(Noop(), IOBuffer())
    @test_throws EOFError read(stream, UInt8)
    @test_throws EOFError unsafe_read(stream, pointer(Vector{UInt8}(undef, 3)), 3)
    close(stream)

    stream = TranscodingStream(Noop(), IOBuffer("foobar"), bufsize=1)
    @test read(stream, UInt8) === UInt8('f')
    data = Vector{UInt8}(undef, 5)
    unsafe_read(stream, pointer(data), 5) === nothing
    @test data == b"oobar"
    close(stream)

    sink = IOBuffer()
    stream = TranscodingStream(Noop(), sink)
    @test write(stream, "foo") === 3
    @test occursin("mode=write", repr(stream))
    flush(stream)
    @test take!(sink) == b"foo"
    close(stream)

    data = rand(UInt8, 100_000)
    sink = IOBuffer()
    stream = TranscodingStream(Noop(), sink)
    for i in 1:10_000
        @assert write(stream, data[10(i-1)+1:10i]) == 10
    end
    flush(stream)
    @test take!(sink) == data
    close(stream)

    stream = TranscodingStream(Noop(), IOBuffer(b"foobarbaz"))
    @test position(stream) === 0
    read(stream, UInt8)
    @test position(stream) === 1
    read(stream)
    @test position(stream) === 9

    data = collect(0x00:0x0f)
    stream = TranscodingStream(Noop(), IOBuffer(data))
    @test !ismarked(stream)
    mark(stream)
    @test ismarked(stream)
    @test [read(stream, UInt8) for _ in 1:3] == data[1:3]
    reset(stream)
    @test !ismarked(stream)
    @test [read(stream, UInt8) for _ in 1:3] == data[1:3]
    mark(stream)
    @test ismarked(stream)
    unmark(stream)
    @test !ismarked(stream)
    close(stream)

    stream = TranscodingStream(Noop(), IOBuffer(b"foobarbaz"))
    seek(stream, 2)
    @test read(stream, 3) == b"oba"
    seek(stream, 0)
    @test read(stream, 3) == b"foo"
    seekstart(stream)
    @test read(stream, 3) == b"foo"
    seekend(stream)
    @test eof(stream)
    close(stream)

    data = collect(0x00:0x0f)
    stream = TranscodingStream(Noop(), IOBuffer(data))
    @test read(stream, UInt8) == data[1]
    skip(stream, 1)
    @test read(stream, UInt8) == data[3]
    skip(stream, 5)
    @test read(stream, UInt8) == data[9]
    skip(stream, 7)
    @test eof(stream)
    close(stream)

    # skip offset > bufsize
    data = collect(0x00:0x0f)
    stream = TranscodingStream(Noop(), IOBuffer(data), bufsize=2)
    @test read(stream, UInt8) == data[1]
    skip(stream, 4)
    @test read(stream, UInt8) == data[6]
    skip(stream, 3)
    @test read(stream, UInt8) == data[10]
    skip(stream, 6)
    @test eof(stream)
    close(stream)

    stream = TranscodingStream(Noop(), IOBuffer("foo"))
    read(stream, UInt8)
    @test_throws ArgumentError skip(stream, -1)
    skip(stream, 100)
    @test eof(stream)
    close(stream)

    stream = TranscodingStream(Noop(), IOBuffer("foo"))
    out = zeros(UInt8, 3)
    @test bytesavailable(stream) == 0
    @test TranscodingStreams.unsafe_read(stream, pointer(out), 10) == 3
    @test out == b"foo"
    close(stream)

    data = rand(UInt8, 1999)
    # unmarked
    stream = TranscodingStream(Noop(), IOBuffer(data), bufsize=7)
    @test hash(read(stream)) == hash(data)
    @test length(stream.state.buffer1.data) == 7
    # marked
    stream = TranscodingStream(Noop(), IOBuffer(data), bufsize=7)
    mark(stream)
    @test hash(read(stream)) == hash(data)
    @test hash(stream.state.buffer1.data[1:length(data)]) == hash(data)
    close(stream)

    stream = NoopStream(NoopStream(IOBuffer("foobar")))
    @test read(stream) == b"foobar"
    close(stream)

    stream = NoopStream(NoopStream(NoopStream(IOBuffer("foobar"))))
    @test read(stream) == b"foobar"
    close(stream)

    # Two buffers are the same object.
    stream = NoopStream(IOBuffer("foo"))
    @test stream.state.buffer1 === stream.state.buffer2

    # Nested NoopStreams share the same buffer.
    s0 = IOBuffer("foo")
    s1 = NoopStream(s0)
    s2 = NoopStream(s1)
    s3 = NoopStream(s2)
    @test s1.state.buffer1 === s2.state.buffer1 === s3.state.buffer1 ===
          s1.state.buffer2 === s2.state.buffer2 === s3.state.buffer2

    stream = TranscodingStream(Noop(), IOBuffer(b"foobar"))
    @test TranscodingStreams.stats(stream).in === Int64(0)
    @test TranscodingStreams.stats(stream).out === Int64(0)
    read(stream)
    @test TranscodingStreams.stats(stream).in === Int64(6)
    @test TranscodingStreams.stats(stream).out === Int64(6)
    close(stream)

    stream = TranscodingStream(Noop(), IOBuffer())
    @test TranscodingStreams.stats(stream).in === Int64(0)
    @test TranscodingStreams.stats(stream).out === Int64(0)
    write(stream, b"foobar")
    flush(stream)
    @test TranscodingStreams.stats(stream).in === Int64(6)
    @test TranscodingStreams.stats(stream).out === Int64(6)
    close(stream)

    stream = TranscodingStream(Noop(), IOBuffer())
    @test stream.state.mode == :idle
    @test write(stream) == 0
    @test stream.state.mode == :write
    close(stream)

    stream = NoopStream(IOBuffer("foobar"))
    @test bytesavailable(stream) === 0
    @test readavailable(stream) == b""
    @test read(stream, UInt8) === UInt8('f')
    @test bytesavailable(stream) === 5
    @test readavailable(stream) == b"oobar"
    close(stream)

    data = b""
    @test transcode(Noop, data)  == data
    @test transcode(Noop, data) !== data
    data = b"foo"
    @test transcode(Noop, data)  == data
    @test transcode(Noop, data) !== data

    data = b""
    @test transcode(Noop(), data)  == data
    @test transcode(Noop(), data) !== data
    data = b"foo"
    @test transcode(Noop(), data)  == data
    @test transcode(Noop(), data) !== data

    TranscodingStreams.test_roundtrip_transcode(Noop, Noop)
    TranscodingStreams.test_roundtrip_read(NoopStream, NoopStream)
    TranscodingStreams.test_roundtrip_write(NoopStream, NoopStream)
    TranscodingStreams.test_roundtrip_lines(NoopStream, NoopStream)

    # switch write => read
    stream = NoopStream(IOBuffer(b"foobar", read=true, write=true))
    @test_throws ArgumentError begin
        write(stream, b"xyz")
        read(stream, 3)
    end

    # switch read => write
    stream = NoopStream(IOBuffer(b"foobar", read=true, write=true))
    @test_throws ArgumentError begin
        read(stream, 3)
        write(stream, b"xyz")
    end

    stream = NoopStream(IOBuffer(""))
    @test TranscodingStreams.unread(stream, b"foo") === nothing
    @test read(stream, 3) == b"foo"
    close(stream)

    stream = NoopStream(IOBuffer("foo"))
    @test read(stream, 3) == b"foo"
    @test TranscodingStreams.unread(stream, b"bar") === nothing
    @test read(stream, 3) == b"bar"
    close(stream)

    stream = NoopStream(IOBuffer("foobar"))
    @test TranscodingStreams.unread(stream, b"baz") === nothing
    @test read(stream, 3) == b"baz"
    @test read(stream, 3) == b"foo"
    @test read(stream, 3) == b"bar"
    @test eof(stream)
    close(stream)

    stream = NoopStream(IOBuffer("foobar"))
    @test read(stream, 3) == b"foo"
    @test TranscodingStreams.unread(stream, b"baz") === nothing
    @test read(stream, 3) == b"baz"
    @test read(stream, 3) == b"bar"
    @test eof(stream)
    close(stream)

    stream = NoopStream(IOBuffer("foobar"))
    @test read(stream, 3) == b"foo"
    @test read(stream, 3) == b"bar"
    @test TranscodingStreams.unread(stream, b"baz") === nothing
    @test read(stream, 3) == b"baz"
    @test eof(stream)
    close(stream)

    stream = NoopStream(IOBuffer("foobar"))
    @test_throws ArgumentError TranscodingStreams.unsafe_unread(stream, pointer(b"foo"), -1)
    close(stream)

    stream = NoopStream(IOBuffer(""))
    unsafe_write(stream, C_NULL, 0)
    @test eof(stream)  # write
    close(stream)
    @test eof(stream)  # close

    @testset "readuntil" begin
        stream = NoopStream(IOBuffer(""))
        data = readuntil(stream, 0x00)
        @test data isa Vector{UInt8}
        @test isempty(data)

        stream = NoopStream(IOBuffer("foo,bar"))
        @test readuntil(stream, UInt8(',')) == b"foo"
        @test read(stream) == b"bar"

        stream = NoopStream(IOBuffer("foo,bar"))
        @test readuntil(stream, UInt8(','), keep = false) == b"foo"
        @test read(stream) == b"bar"

        stream = NoopStream(IOBuffer("foo,bar"))
        @test readuntil(stream, UInt8(','), keep = true) == b"foo,"
        @test read(stream) == b"bar"
    end

    @test_throws ArgumentError NoopStream(IOBuffer(), bufsize=0)
    @test_throws ArgumentError NoopStream(let s = IOBuffer(); close(s); s; end)
    @test_throws ArgumentError TranscodingStream(Noop(), IOBuffer(), bufsize=0)
    @test_throws ArgumentError TranscodingStream(Noop(), IOBuffer(), sharedbuf=true)
end

ext: jl


text:
## Define DREAMR types (based on US, UA)
@enum DREAMR_MODETYPE FLIGHT=1 RIDE=2
@enum HOP_ACTION HOPON=1 HOPOFF=2 STAY=3

## Mode Info
const DREAMR_MODES = [FLIGHT, RIDE]
const DREAMR_GOAL_MODE = FLIGHT
const DREAMR_MODE_ACTIONS = [HOPON, HOPOFF, STAY]

const HOP_INFLATION_PARAMETER = 2.0

mutable struct DREAMRModeAction
    mode_action::HOP_ACTION
    car_id::String
end

Base.isequal(dma1::DREAMRModeAction, dma2::DREAMRModeAction) = isequal(dma1.mode_action, dma2.mode_action) && isequal(dma1.car_id, dma2.car_id)

DREAMRModeAction(mode_action::HOP_ACTION) = DREAMRModeAction(mode_action, "")

## Car route stuff
struct DREAMRVertexMetadata
    car_id::String
    vertex_id::String
end

Base.zero(::Type{DREAMRVertexMetadata}) = DREAMRVertexMetadata("", "")

struct DREAMRBookkeeping
    active_car_to_idx_ranges::Dict{String, Vector{Vector{Int64}}} # Maps car name to range of vertices that represent hopon
    route_vert_id_to_idx_pair::Dict{String, Vector{Int64}} # Maps a specific route ID (car-id+"-"+vertex_id) to the indices of hopon and hopoff versions
    goal_vertex_idx::Int64
end
DREAMRBookkeeping() = DREAMRBookkeeping(Dict{String, Vector{Vector{Int64}}}(), Dict{String, Vector{Int64}}(), 0)

## Context
mutable struct DREAMRContextSet
    curr_epoch::Int64
    epochs_dict::Dict
    car_id::String
end

function get_dreamr_episode_context(ep_dict::Dict)
    return DREAMRContextSet(0, ep_dict, "")
end

const DREAMRStateType{US} = CMSSPState{DREAMR_MODETYPE, US} where {US <: UAVState}
const DREAMRActionType{UA} = CMSSPAction{DREAMRModeAction, UA} where {UA <: UAVAction}
const DREAMRCMSSPType{US,UA} = CMSSP{DREAMR_MODETYPE, US, DREAMRModeAction, UA, Parameters, DREAMRContextSet} where {US <: UAVState, UA <: UAVAction}
const DREAMRCFMDPType{US,UA} = ModalFinHorMDP{DREAMR_MODETYPE, US, UA, Parameters} where {US <: UAVState, UA <: UAVAction}
const DREAMRUFMDPType{US,UA} = ModalInfHorMDP{DREAMR_MODETYPE, US, UA, Parameters} where {US <: UAVState, UA <: UAVAction}
const DREAMRHopoffMDPType = ModalFinHorMDP{DREAMR_MODETYPE, Nothing, Nothing, Parameters}
const DREAMROpenLoopVertex{US} = OpenLoopVertex{DREAMR_MODETYPE, US, DREAMRModeAction, DREAMRVertexMetadata} where {US <: UAVState}
const DREAMRGraphTracker{US,UA,RNG} = GraphTracker{DREAMR_MODETYPE, US, DREAMRModeAction, DREAMRVertexMetadata, DREAMRBookkeeping, RNG} where {US <: UAVState, UA <: UAVAction, RNG <: AbstractRNG}
const DREAMRModePair = Tuple{DREAMR_MODETYPE, DREAMR_MODETYPE}
const DREAMRSolverType{US,UA,RNG} = HHPCSolver{DREAMR_MODETYPE, US, DREAMRModeAction, DREAMRVertexMetadata, DREAMRBookkeeping, RNG} where {US <: UAVState, UA <: UAVAction, RNG <: AbstractRNG}

# This is only really needed for baseline
function get_dreamr_actions(::Type{UA}, params::Parameters) where {UA <: UAVAction}

    actions = Vector{DREAMRActionType{UA}}(undef, 0)
    idx = 1
    for dma in DREAMR_MODE_ACTIONS
        push!(actions, DREAMRActionType{UA}(DREAMRModeAction(dma), idx))
        idx += 1
    end

    uav_dynamics_actions = get_uav_dynamics_actions(UA, params)
    for uavda in uav_dynamics_actions
        push!(actions, DREAMRActionType{UA}(uavda, idx))
        idx += 1
    end

    return actions
end

function get_dreamr_mdp()

    T = zeros(2, 3, 2)
    # Flight + hopon = ride
    T[1, 1, 2] = 1.0
    # Ride + hopoff = flight
    T[2, 2, 1] = 1.0
    # Ride + STAY = Ride
    T[2, 3, 2] = 1.0

    R = zeros(2, 3)

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

function create_dreamr_cmssp(::Type{US}, ::Type{UA}, context_set::DREAMRContextSet,
                             goal_point::Point, params::Parameters) where {US <: UAVState, UA <: UAVAction}

    actions = get_dreamr_actions(UA, params)
    switch_mdp = get_dreamr_mdp()
    goal_state = DREAMRStateType{US}(FLIGHT, get_state_at_rest(US, goal_point))

    return DREAMRCMSSPType{US,UA}(actions, DREAMR_MODES, switch_mdp, goal_state, params, context_set)
end

function POMDPs.actions(cmssp::DREAMRCMSSPType, state::DREAMRStateType)

    if cmssp.curr_context_set.curr_epoch >= length(collect(keys(cmssp.curr_context_set.epochs_dict)))
        return cmssp.actions[4:end]
    end

    if state.mode == FLIGHT
        return [cmssp.actions[1]; cmssp.actions[4:end]]
    else
        return cmssp.actions[2:3]
    end
end

function set_dreamr_goal!(cmssp::DREAMRCMSSPType{US,UA}, goal_point::Point) where {US <: UAVState, UA <: UAVAction}
    cmssp.goal_state = DREAMRStateType{US}(FLIGHT, get_state_at_rest(US, goal_point))
end

function POMDPs.isterminal(cmssp::DREAMRCMSSPType, state::DREAMRStateType)

    curr_point = get_position(state.continuous)
    curr_speed = get_speed(state.continuous)

    goal_point = get_position(cmssp.goal_state.continuous)

    dist = point_dist(curr_point, goal_point)

    params = cmssp.params

    return (state.mode == FLIGHT) && (dist < 2.0*params.time_params.MDP_TIMESTEP*params.scale_params.HOP_DISTANCE_THRESHOLD) &&
            (curr_speed < params.scale_params.XYDOT_HOP_THRESH)
end

function POMDPs.isterminal(mdp::Union{DREAMRCFMDPType{US,UA},DREAMRUFMDPType{US,UA}}, relative_state::US) where {US <: UAVState, UA <: UAVAction}

    curr_pos_norm = point_norm(get_position(relative_state))
    curr_speed = get_speed(relative_state)

    return (curr_pos_norm < mdp.params.time_params.MDP_TIMESTEP*mdp.params.scale_params.HOP_DISTANCE_THRESHOLD) &&
        (curr_speed < mdp.params.scale_params.XYDOT_HOP_THRESH)
end

function HHPC.get_relative_state(mdp::Union{DREAMRCFMDPType{US,UA},DREAMRUFMDPType{US,UA}}, source::US, target::US) where {US <: UAVState, UA <: UAVAction}
    return rel_state(source, target)
end

## Conversion for multirotor - can add others as needed
function POMDPs.convert_s(::Type{Vector{Float64}}, s::MultiRotorUAVState,
                          mdp::DREAMRCFMDPType{MultiRotorUAVState, MultiRotorUAVAction})
    v = [s.x, s.y, s.xdot, s.ydot]
end

function POMDPs.convert_s(::Type{MultiRotorUAVState}, v::AbstractVector{Float64},
                          mdp::DREAMRCFMDPType{MultiRotorUAVState, MultiRotorUAVAction})
    s = MultiRotorUAVState(v[1], v[2], v[3], v[4])
end

function POMDPs.convert_s(::Type{Vector{Float64}}, s::MultiRotorUAVState,
                          mdp::DREAMRUFMDPType{MultiRotorUAVState,MultiRotorUAVAction})
    v = [s.x, s.y, s.xdot, s.ydot]
end

function POMDPs.convert_s(::Type{MultiRotorUAVState}, v::AbstractVector{Float64},
                          mdp::DREAMRUFMDPType{MultiRotorUAVState,MultiRotorUAVAction})
    s = MultiRotorUAVState(v[1], v[2], v[3], v[4])
end

function POMDPs.reward(mdp::DREAMRCFMDPType{US,UA}, state::US, cont_action::UA, statep::US) where {US <: UAVState, UA <: UAVAction}

    cost = mdp.params.cost_params.TIME_COEFFICIENT*mdp.params.time_params.MDP_TIMESTEP
    cost += dynamics_cost(mdp.params, state, statep)

    return -cost
end

function POMDPs.reward(mdp::DREAMRUFMDPType{US,UA}, state::US, cont_action::UA, statep::US) where {US <: UAVState, UA <: UAVAction}

    cost = mdp.params.cost_params.TIME_COEFFICIENT*mdp.params.time_params.MDP_TIMESTEP
    cost += dynamics_cost(mdp.params, state, statep)

    if isterminal(mdp, statep)
        cost -= mdp.terminal_reward
    end

    return -cost
end

# FH version of generate_sr
function POMDPs.generate_sr(mdp::Union{DREAMRCFMDPType{US,UA},DREAMRUFMDPType{US,UA}}, state::US, cont_action::UA,
                            rng::RNG=Random.GLOBAL_RNG) where {US <: UAVState, UA <: UAVAction, RNG <: AbstractRNG}

    statep = next_state(mdp.params, state, cont_action, rng)
    reward = POMDPs.reward(mdp, state, cont_action, statep)

    return (statep, reward)
end

function get_cf_mdp(::Type{US}, ::Type{UA},
                    params::Parameters, beta::Float64) where {US <:UAVState, UA <: UAVAction}

    uav_dynamics_actions = get_uav_dynamics_actions(UA, params)

    return DREAMRCFMDPType{US,UA}(FLIGHT, params, uav_dynamics_actions, beta, params.time_params.HORIZON_LIM)
end

function get_uf_mdp(::Type{US}, ::Type{UA}, params::Parameters) where {US <:UAVState, UA <: UAVAction}

    uav_dynamics_actions = get_uav_dynamics_actions(UA, params)

    return DREAMRUFMDPType{US,UA}(FLIGHT, uav_dynamics_actions, params, params.cost_params.FLIGHT_REACH_REWARD) # beta = 1 and horizon = 0
end

function HHPC.expected_reward(mdp::Union{DREAMRCFMDPType{US,UA},DREAMRUFMDPType{US,UA}},
                              state::US, cont_action::UA, rng::RNG=Random.GLOBAL_RNG) where {US <: UAVState, UA <: UAVAction, RNG <: AbstractRNG}

    params = mdp.params
    avg_rew = 0.0

    for i = 1:params.scale_params.MC_GENERATIVE_NUMSAMPLES
        (_, r) = generate_sr(mdp, state, cont_action, rng)
        avg_rew += r
    end

    avg_rew /= params.scale_params.MC_GENERATIVE_NUMSAMPLES
    return avg_rew
end

function HHPC.generate_goal_vertex_set!(cmssp::DREAMRCMSSPType{US,UA}, vertex::DREAMROpenLoopVertex{US},
                                        graph_tracker::DREAMRGraphTracker,
                                        rng::RNG=Random.GLOBAL_RNG) where {US <: UAVState, UA <: UAVAction, RNG <: AbstractRNG}

    # If current goal vertex idx is 0, add goal and send it as nbr
    if graph_tracker.curr_goal_idx == 0
        graph_tracker.curr_goal_idx = num_vertices(graph_tracker.curr_graph)+1 # Since it will be added
        goal_vert = DREAMROpenLoopVertex{US}(cmssp.goal_state, cmssp.mode_actions[1])
        return [goal_vert], []
    else
        # Just return current goal vertex idx
        return [], [graph_tracker.curr_goal_idx]
    end
end

function get_arrival_time_distribution(curr_timestep::Int64, mean_time::Float64, params::Parameters, rng::RNG=Random.GLOBAL_RNG) where {RNG <: AbstractRNG}

    # If very close to time, collapse around estimate
    if mean_time / params.time_params.MDP_TIMESTEP - curr_timestep > 2.0
        time_dist = Normal(mean_time, params.time_params.CAR_TIME_STD)
    else
        time_dist = Normal(mean_time, params.scale_params.EPSILON)
    end

    temp_tp_dict = Dict{Int64,Float64}()

    for j = 1:params.time_params.MC_TIME_NUMSAMPLES
        tval = rand(rng, time_dist) / params.time_params.MDP_TIMESTEP

        low = convert(Int64, floor(tval))
        high = convert(Int64, ceil(tval))
        low_wt = tval - floor(tval)

        cum_low_wt = get(temp_tp_dict, low+1, 0.0)
        temp_tp_dict[low+1] = cum_low_wt + low_wt
        cum_high_wt = get(temp_tp_dict, high+1, 0.0)
        temp_tp_dict[high+1] = cum_high_wt + 1.0 - low_wt
    end

    vals = Vector{Int64}(undef, 0)
    probs = Vector{Float64}(undef, 0)

    for (tval, prob) in temp_tp_dict
        push!(vals, tval)
        push!(probs, prob/params.time_params.MC_TIME_NUMSAMPLES)
    end

    return TPDistribution(vals, probs)
end

function HHPC.generate_bridge_vertex_set!(cmssp::DREAMRCMSSPType, vertex::DREAMROpenLoopVertex,
                                         mode_pair::DREAMRModePair,
                                         graph_tracker::DREAMRGraphTracker,
                                         action::DREAMRModeAction,
                                         rng::RNG=Random.GLOBAL_RNG) where {RNG <: AbstractRNG}

    # NOTE - Will never actually generate anything; only get neighbour idxs
    context_set = cmssp.curr_context_set
    params = cmssp.params

    nbrs_to_add = Vector{Int64}(undef, 0)

    if vertex.state.mode == FLIGHT # If it is a drone or post-hopoff vertex

        @assert mode_pair == (FLIGHT, RIDE)

        # Iterate over cars
        for (car_id, ranges) in graph_tracker.bookkeeping.active_car_to_idx_ranges

            # If this vertex just hopped off a car, don't hop on to another!
            # if car_id == vertex.metadata.car_id
            #     @warn "Hopoff searching for hopon"
            #     @show vertex
            #     readline()
            #     continue
            # end

            hopon_vert_range = ranges[1]

            for next_idx in hopon_vert_range[1]:hopon_vert_range[2]

                # Check if max drone speed not enough
                next_vert = graph_tracker.curr_graph.vertices[next_idx]
                state_dist = point_dist(get_position(vertex.state.continuous), get_position(next_vert.state.continuous))

                if params.scale_params.MAX_DRONE_SPEED*(mean(next_vert.tp) - mean(vertex.tp)) >
                    state_dist*params.scale_params.VALID_FLIGHT_EDGE_DIST_RATIO
                    push!(nbrs_to_add, next_idx)
                end
            end
        end
    else
        # It is a car route vertex
        # Only add hopoff vertices for that car
        @assert mode_pair == (RIDE, FLIGHT)
        @assert vertex.metadata.car_id != ""

        hopoff_vert_range = graph_tracker.bookkeeping.active_car_to_idx_ranges[vertex.metadata.car_id][2]

        # Loop through range and add ones that come after current
        for hopoff_idx = hopoff_vert_range[1]:hopoff_vert_range[2]

            next_vert = graph_tracker.curr_graph.vertices[hopoff_idx]

            # Remember vertex is a hopon vertex, so can only add hopoff ones; compare by tp
            if mean(next_vert.tp) > mean(vertex.tp)
                push!(nbrs_to_add, hopoff_idx)
            end
        end
        # if context_set.car_id != ""
        #     @warn "Popping start ride vertex"
        #     @show vertex
        #     @show nbrs_to_add
        #     readline()
        # end
    end

    return ([], nbrs_to_add)
end

function HHPC.generate_next_valid_modes(cmssp::DREAMRCMSSPType, mode::DREAMR_MODETYPE)

    if mode == FLIGHT
        return [(DREAMRModeAction(HOPON), RIDE)]
    else
        return [(DREAMRModeAction(HOPOFF), FLIGHT)]
    end
end

function HHPC.update_vertices_with_context!(cmssp::DREAMRCMSSPType{US,UA}, graph_tracker::DREAMRGraphTracker,
                                            curr_timestep::Int64) where {US <: UAVState, UA <: UAVAction}

    # Iterate over cars, add if new or update if old
    context_set = cmssp.curr_context_set

    curr_epoch_dict = context_set.epochs_dict[string(context_set.curr_epoch)]
    epoch_cars = curr_epoch_dict["car-info"]
    params = cmssp.params

    keys_to_delete = Vector{String}(undef, 0)

    for (car_id, car_info) in epoch_cars

        # Update car route if it exists
        if haskey(graph_tracker.bookkeeping.active_car_to_idx_ranges, car_id)

            route_info = car_info["route"]

            if route_info != nothing

                # Update times of all future vertices
                sorted_route = sort(collect(route_info), by=x->parse(Float64, x[1]))

                for (id, timept) in sorted_route

                    timeval = timept[2]
                    tp_dist = get_arrival_time_distribution(curr_timestep, timeval, params, graph_tracker.rng)

                    # Update the hopon and hopoff vert tps
                    vert_id = string(car_id, "-", id)
                    graph_tracker.curr_graph.vertices[graph_tracker.bookkeeping.route_vert_id_to_idx_pair[vert_id][1]].tp = tp_dist
                    graph_tracker.curr_graph.vertices[graph_tracker.bookkeeping.route_vert_id_to_idx_pair[vert_id][2]].tp = tp_dist
                end

                # If first waypoint of route updated, change BOTH hopon and hopoff first idxs
                first_route_idx_hopon = graph_tracker.bookkeeping.route_vert_id_to_idx_pair[string(car_id,"-",sorted_route[1][1])][1]
                first_route_idx_hopoff = graph_tracker.bookkeeping.route_vert_id_to_idx_pair[string(car_id,"-",sorted_route[1][1])][2]

                if graph_tracker.bookkeeping.active_car_to_idx_ranges[car_id][1][1] != first_route_idx_hopon
                    # @info "Car ",car_id," has updated its next route point"
                    graph_tracker.bookkeeping.active_car_to_idx_ranges[car_id][1][1] = first_route_idx_hopon
                    graph_tracker.bookkeeping.active_car_to_idx_ranges[car_id][2][1] = first_route_idx_hopoff
                end
            else
                push!(keys_to_delete, car_id)
            end
        else

            # New car
            route_info = car_info["route"]

            if route_info != nothing

                graph_tracker.bookkeeping.active_car_to_idx_ranges[car_id] = [[0.,0.], [0.,0.]]

                # Add two variants - one for flight to hopon and one for hopon to hopoff
                prev_last_route_idx = num_vertices(graph_tracker.curr_graph)
                last_route_idx = prev_last_route_idx

                temp_vertices = Vector{DREAMROpenLoopVertex{US}}(undef, 0)
                temp_vert_ids = Vector{String}(undef, 0)

                # First add hopon variants
                for (id, timept) in sort(collect(route_info), by=x->parse(Float64, x[1]))

                    last_route_idx += 1

                    # Extract time and point respectively
                    timeval = timept[2]
                    waypt = Point(timept[1][1], timept[1][2])

                    # Get the ETA distribution
                    tp_dist = get_arrival_time_distribution(curr_timestep, timeval, params, graph_tracker.rng)

                    # Set pre and post bridge state to waypoint position
                    pre_bridge_state = DREAMRStateType{US}(FLIGHT, get_state_at_rest(US, waypt))
                    post_bridge_state = DREAMRStateType{US}(RIDE, get_state_at_rest(US, waypt))

                    metadata = DREAMRVertexMetadata(car_id, id)

                    new_vertex = DREAMROpenLoopVertex{US}(post_bridge_state, pre_bridge_state, DREAMRModeAction(HOPON), tp_dist, metadata)
                    add_vertex!(graph_tracker.curr_graph, new_vertex)

                    push!(temp_vertices, new_vertex)

                    # Now update bookkeeping vert_id_to_idx
                    vert_id = string(car_id, "-", id)
                    graph_tracker.bookkeeping.route_vert_id_to_idx_pair[vert_id] = [last_route_idx, 0]

                    push!(temp_vert_ids, vert_id)
                end

                # Set hopon range
                graph_tracker.bookkeeping.active_car_to_idx_ranges[car_id][1] = [prev_last_route_idx+1, last_route_idx]

                # NOW DO SAME FOR HOPOFF VERSIONS
                prev_last_route_idx = last_route_idx

                for (i, tmpver) in enumerate(temp_vertices)

                    last_route_idx += 1

                    # can copy data from tmpver, just change mode and action
                    cont_state = tmpver.pre_bridge_state.continuous

                    pre_bridge_state = DREAMRStateType{US}(RIDE, cont_state)
                    post_bridge_state = DREAMRStateType{US}(FLIGHT, cont_state)

                    new_vertex = DREAMROpenLoopVertex{US}(post_bridge_state, pre_bridge_state, DREAMRModeAction(HOPOFF), tmpver.tp, tmpver.metadata)
                    add_vertex!(graph_tracker.curr_graph, new_vertex)

                    # Update hopoff pair member
                    vert_id = temp_vert_ids[i]
                    graph_tracker.bookkeeping.route_vert_id_to_idx_pair[vert_id][2] = last_route_idx
                end

                # Set hopoff range
                graph_tracker.bookkeeping.active_car_to_idx_ranges[car_id][2] = [prev_last_route_idx+1, last_route_idx]

                # @info "Car", car_id, "has been added!"
            end
        end
    end

    # @show graph_tracker.bookkeeping.active_car_to_idx_ranges["car-115"]
    # @show keys_to_delete

    # Delete inactive keys
    for kd in keys_to_delete
        delete!(graph_tracker.bookkeeping.active_car_to_idx_ranges, kd)
    end
end

function HHPC.update_next_target!(cmssp::DREAMRCMSSPType, solver::DREAMRSolverType, curr_timestep::Int64)

    context_set = cmssp.curr_context_set
    next_target = solver.graph_tracker.curr_graph.vertices[solver.graph_tracker.curr_soln_path_idxs[2]]

    car_id = next_target.metadata.car_id
    vertex_id = next_target.metadata.vertex_id

    car_route_info = context_set.epochs_dict[string(context_set.curr_epoch)]["car-info"][car_id]["route"]

    if haskey(car_route_info, vertex_id) == false
        return false # Invalid update as car has disappeared
    end

    timeval = car_route_info[vertex_id][2]
    next_target.tp = get_arrival_time_distribution(curr_timestep, timeval, cmssp.params, solver.graph_tracker.rng)

    return true
end

function HHPC.update_context_set!(cmssp::DREAMRCMSSPType, rng::RNG=Random.GLOBAL_RNG) where {RNG <: AbstractRNG}
    cmssp.curr_context_set.curr_epoch += 1
end

function HHPC.simulate_cmssp!(cmssp::DREAMRCMSSPType{US,UA}, state::DREAMRStateType{US}, a::Union{Nothing,DREAMRActionType{UA}}, t::Int64,
                              rng::RNG=Random.GLOBAL_RNG) where {US <: UAVState, UA <: UAVAction, RNG <: AbstractRNG}

    # Update context
    timeout = false

    context_set = cmssp.curr_context_set
    params = cmssp.params

    if haskey(context_set.epochs_dict, string(context_set.curr_epoch))
        curr_epoch_dict = context_set.epochs_dict[string(context_set.curr_epoch)]
        epoch_cars = curr_epoch_dict["car-info"]
    else
        # No more car actions after 360
        timeout = true
        return (state, 0.0, false, timeout)
    end

    reward = -params.cost_params.TIME_COEFFICIENT*params.time_params.MDP_TIMESTEP
    failed_mode_switch = false

    # FIRST CHECK FOR nothing (abort) - just return same state
    if a == nothing
        return (state, reward, failed_mode_switch, timeout)
    end

    curr_mode = state.mode
    curr_cont_state = state.continuous

    if typeof(a.action) <: UAVAction

        @assert curr_mode == FLIGHT

        new_uav_state = next_state(params, curr_cont_state, a.action, rng)
        reward += -dynamics_cost(params, curr_cont_state, new_uav_state)
        next_dreamr_state = DREAMRStateType{US}(FLIGHT, new_uav_state)

    else
        if a.action.mode_action == HOPON

            @assert curr_mode == FLIGHT

            reward += -cmssp.modeswitch_mdp.R[1, 1]

            hopon_car_id = a.action.car_id

            # Hop on to specific car
            if hopon_car_id != ""

                # Check against car_position
                car_pos = Point(epoch_cars[hopon_car_id]["pos"][1], epoch_cars[hopon_car_id]["pos"][2])

                uav_pos = get_position(curr_cont_state)
                uav_speed = get_speed(curr_cont_state)

                car_uav_dist = point_dist(car_pos, uav_pos)

                if car_uav_dist < 2.0*params.scale_params.HOP_DISTANCE_THRESHOLD*params.time_params.MDP_TIMESTEP &&
                    uav_speed < params.scale_params.XYDOT_HOP_THRESH

                    @warn "Successful hop on to", hopon_car_id, " at epoch ", context_set.curr_epoch
                    new_uav_state = get_state_at_rest(US, car_pos)
                    next_dreamr_state = DREAMRStateType{US}(RIDE, new_uav_state)
                    context_set.car_id = hopon_car_id
                else

                    @warn "Failed to hop on to car!"
                    @debug car_pos
                    failed_mode_switch = true
                    next_dreamr_state = state
                    context_set.car_id = ""
                end
            else
                # Undirected - search over all and update if possible
                next_dreamr_state = state
                hop_success = false
                min_car_uav_dist = Inf

                for (car_id, car_info) in epoch_cars

                    if car_info["route"] != nothing

                        car_pos = Point(car_info["pos"][1], car_info["pos"][2])

                        uav_pos = get_position(curr_cont_state)
                        uav_speed = get_speed(curr_cont_state)

                        car_uav_dist = point_dist(car_pos, uav_pos)

                        if car_uav_dist < min_car_uav_dist
                            min_car_uav_dist = car_uav_dist
                        end

                        if car_uav_dist < HOP_INFLATION_PARAMETER*params.scale_params.HOP_DISTANCE_THRESHOLD*params.time_params.MDP_TIMESTEP &&
                                uav_speed < params.scale_params.XYDOT_HOP_THRESH

                            new_uav_state = get_state_at_rest(US, car_pos)
                            next_dreamr_state = DREAMRStateType{US}(RIDE, new_uav_state)
                            context_set.car_id = car_id

                            hop_success = true
                            @info "Successful hop on to", car_id
                            break
                        end
                    end
                end

                if hop_success == false
                    # @warn "Failed to hop on to car!"
                    reward += -params.cost_params.INVALID_HOP_PENALTY
                end
            end
        elseif a.action.mode_action == HOPOFF

            @assert curr_mode == RIDE

            # Set drone to car's position
            reward += -cmssp.modeswitch_mdp.R[2, 2]

            hopoff_car_id = context_set.car_id
            car_pos = Point(epoch_cars[hopoff_car_id]["pos"][1], epoch_cars[hopoff_car_id]["pos"][2])
            new_uav_state = get_state_at_rest(US, car_pos)

            next_dreamr_state = DREAMRStateType{US}(FLIGHT, new_uav_state)
            context_set.car_id = ""
        else
            reward += -cmssp.modeswitch_mdp.R[2, 3]
            @assert curr_mode == RIDE

            car_pos = Point(epoch_cars[context_set.car_id]["pos"][1], epoch_cars[context_set.car_id]["pos"][2])
            new_uav_state = get_state_at_rest(US, car_pos)

            next_dreamr_state = DREAMRStateType{US}(RIDE, new_uav_state)
        end
    end

    return (next_dreamr_state, reward, failed_mode_switch, timeout)
end

function HHPC.display_context_future(context_set::DREAMRContextSet, future_epoch::Int64)
    println(context_set.epochs_dict[future_epoch])
end

# Return the vertex bridging action and metadata car id
function HHPC.get_bridging_action(vertex::DREAMROpenLoopVertex)
    return DREAMRModeAction(vertex.bridging_action.mode_action, vertex.metadata.car_id)
end

function get_ride_mdp(params::Parameters)
    return DREAMRHopoffMDPType(RIDE, params, Vector{Nothing}(undef, 0))
end

# Create a fake deterministic policy for hopoff - override methods
struct DREAMRDeterministicPolicy <: Policy
    mdp::Union{MDP,POMDP}
end

DREAMRDeterministicPolicy(params::Parameters) = DREAMRDeterministicPolicy(get_ride_mdp(params))

HHPC.get_mdp(ddp::DREAMRDeterministicPolicy) = ddp.mdp

function HHPC.get_best_intramodal_action(policy::DREAMRDeterministicPolicy, curr_timestep::Int64,
                                         tp_dist::TPDistribution, curr_state::US,
                                         target_state::US) where {US <: UAVState}

    # Just blindly return STAY - HOPOFF will be returned based on time
    return ModalAction(DREAMRModeAction(STAY), 0)
end

function HHPC.horizon_weighted_value(policy::DREAMRDeterministicPolicy, curr_timestep::Int64,
                                    tp_dist::TPDistribution, curr_state::US,
                                    target_state::US) where {US <: UAVState}

    # Just the cost corresponding to time difference
    # Return -Inf for minvalue because will never abort
    params = policy.mdp.params
    return -params.cost_params.TIME_COEFFICIENT*(mean(tp_dist) - curr_timestep), -Inf
end

# Based on the current context set, get the position of a car
# k timesteps into the future. Return nothing if car route inactive by then
function get_future_projected_car_position(car_route_info::Dict, curr_time::Float64, future_time::Float64)

    car_pos = Point(car_route_info["pos"][1], car_route_info["pos"][2])
    time_val = curr_time

    route_info = car_route_info["route"]
    sorted_route = sort(collect(route_info), by=x->parse(Float64, x[1]))

    next_car_pos = Point()
    next_time_val = Inf

    for (id, timept) in sorted_route

        timeval = timept[2]
        if timeval > future_time
            next_time_val = timeval
            next_car_pos = Point(timept[1][1], timept[1][2])
        end
    end

    if next_time_val == Inf
        return nothing
    end

    # Now interpolate car position
    frac = (future_time - curr_time)/(next_time_val - curr_time)
    estimated_point = interpolate(car_pos, next_car_pos, frac)

    return estimated_point
end

# Define specific CMSSP type for MCTS (with root depth)
struct DREAMRMCTSState{US <: UAVState}
    cmssp_state::DREAMRStateType{US}
    car_id::String
    depth_root::Int64
end

function DREAMRMCTSState(cmssp_state::DREAMRStateType)
    return DREAMRMCTSState(cmssp_state, "", 0)
end

mutable struct DREAMRMCTSType{US <: UAVState, UA <: UAVAction} <: POMDPs.MDP{DREAMRMCTSState{US}, DREAMRActionType{UA}}
    cmssp::DREAMRCMSSPType{US,UA}
end

POMDPs.actions(dmcts::DREAMRMCTSType) = dmcts.cmssp.actions
POMDPs.actions(dmcts::DREAMRMCTSType, state::DREAMRMCTSState) = actions(dmcts.cmssp, state.cmssp_state)
POMDPs.n_actions(dmcts::DREAMRMCTSType) = length(dmcts.cmssp.actions)
POMDPs.discount(dmcts::DREAMRMCTSType) = 1.0 # SSP - Undiscounted
POMDPs.actionindex(dmcts::DREAMRMCTSType, a::DREAMRActionType) = a.action_idx

function POMDPs.isterminal(dmcts::DREAMRMCTSType, state::DREAMRMCTSState)
    return isterminal(dmcts.cmssp, state.cmssp_state)
end

# Requirement for MCTS-DPW
# NOTE - This uses the
function POMDPs.generate_sr(dmcts::DREAMRMCTSType{US,UA}, state::DREAMRMCTSState{US},
                            a::DREAMRActionType{UA}, rng::RNG=Random.GLOBAL_RNG) where {US <: UAVState, UA <: UAVAction, RNG <: AbstractRNG}

    cmssp = dmcts.cmssp
    params = cmssp.params
    cmssp_state = state.cmssp_state
    curr_cont_state = cmssp_state.continuous

    reward = -params.cost_params.TIME_COEFFICIENT*params.time_params.MDP_TIMESTEP

    if typeof(a.action) <: UAVAction # Just simulate control

        new_uav_state = next_state(params, curr_cont_state, a.action, rng)
        reward += -dynamics_cost(params, curr_cont_state, new_uav_state)
        next_dreamr_state = DREAMRStateType{US}(FLIGHT, new_uav_state)
        next_dmcts_state = DREAMRMCTSState(next_dreamr_state, "", state.depth_root+1)

    else

        # Is a mode-switch action, propagate context and see what to do
        context_set = cmssp.curr_context_set
        curr_time = context_set.curr_epoch*params.time_params.MDP_TIMESTEP
        future_time = curr_time + (state.depth_root+1)*params.time_params.MDP_TIMESTEP

        curr_epoch_dict = context_set.epochs_dict[string(context_set.curr_epoch)]
        epoch_cars = curr_epoch_dict["car-info"]

        # Loop through all current cars and evaluate
        # DOESN'T MATTER if currently it has an ID
        if a.action.mode_action == HOPON

            # If nothing else, copy current state
            next_dmcts_state = DREAMRMCTSState(cmssp_state, "", state.depth_root+1)
            hop_success = false

            reward += -cmssp.modeswitch_mdp.R[1, 1]

            for (car_id, car_info) in epoch_cars

                if car_info["route"] != nothing

                    car_pos = get_future_projected_car_position(car_info, curr_time, future_time)

                    if car_pos != nothing # Not if car not valid at that future timestep

                        uav_pos = get_position(curr_cont_state)
                        uav_speed = get_speed(curr_cont_state)

                        car_uav_dist = point_dist(car_pos, uav_pos)

                        if car_uav_dist < HOP_INFLATION_PARAMETER*params.scale_params.HOP_DISTANCE_THRESHOLD*params.time_params.MDP_TIMESTEP &&
                                uav_speed < params.scale_params.XYDOT_HOP_THRESH

                            new_uav_state = get_state_at_rest(US, car_pos)
                            next_dreamr_state = DREAMRStateType{US}(RIDE, new_uav_state)

                            next_dmcts_state = DREAMRMCTSState(next_dreamr_state, car_id, state.depth_root+1)

                            hop_success = true

                            @info "Successful hop to ", car_id," in sim!"

                            break
                        end
                    end
                end
            end

            # if hop_success == false
            #     reward += -params.cost_params.INVALID_HOP_PENALTY
            # end
        elseif a.action.mode_action == HOPOFF

            @debug "HOPOFF in SIM!"
            reward += -cmssp.modeswitch_mdp.R[2, 2]

            hopoff_car_id = state.car_id
            car_pos = Point(epoch_cars[hopoff_car_id]["pos"][1], epoch_cars[hopoff_car_id]["pos"][2])
            new_uav_state = get_state_at_rest(US, car_pos)
            next_dreamr_state = DREAMRStateType{US}(FLIGHT, new_uav_state)

            next_dmcts_state = DREAMRMCTSState(next_dreamr_state, "", state.depth_root+1)

        else
            reward += -cmssp.modeswitch_mdp.R[2, 3]
            car_info = epoch_cars[state.car_id]
            @debug "STAY in SIM!"

            car_pos = get_future_projected_car_position(car_info, curr_time, future_time)
            if car_pos != nothing
                new_uav_state = get_state_at_rest(US, car_pos)
                next_dreamr_state = DREAMRStateType{US}(RIDE, new_uav_state)
                next_dmcts_state = DREAMRMCTSState(next_dreamr_state, state.car_id, state.depth_root+1)
            else
                # Car has stopped, must get off
                @debug "STAY cannot be simulated as car has ended by this time"
                reward += -params.cost_params.INVALID_HOP_PENALTY
                car_pos = Point(car_info["pos"][1], car_info["pos"][2])

                new_uav_state = get_state_at_rest(US, car_pos)
                next_dreamr_state = DREAMRStateType{US}(FLIGHT, new_uav_state)
                next_dmcts_state = DREAMRMCTSState(next_dreamr_state, "", state.depth_root+1)
            end
        end
    end

    if isterminal(dmcts, next_dmcts_state)
        reward += params.cost_params.FLIGHT_REACH_REWARD
    end

    return (next_dmcts_state, reward)
end

function estimate_value_dreamr(dmcts::DREAMRMCTSType, state::DREAMRMCTSState, depth::Int64)

    goal_pos = get_position(dmcts.cmssp.goal_state.continuous)
    curr_pos = get_position(state.cmssp_state.continuous)

    value = -0.75*dmcts.cmssp.params.cost_params.FLIGHT_COEFFICIENT*point_dist(curr_pos, goal_pos)

    return value
end

function estimate_value_dreamr(policy::P, dmcts::DREAMRMCTSType, state::DREAMRMCTSState{US}, depth::Int64) where {P <: Policy, US <: UAVState}

    if state.cmssp_state.mode == FLIGHT
        relative_state = rel_state(state.cmssp_state.continuous, dmcts.cmssp.goal_state.continuous)
        value = POMDPs.value(policy, relative_state) - dmcts.cmssp.params.cost_params.FLIGHT_REACH_REWARD
        return value
    else
        # How to compute value of ride state?
        # Look at all car waypoints and compute their value to goal
        # and also consider time to get there
        params = dmcts.cmssp.params
        context_set = dmcts.cmssp.curr_context_set

        curr_time = context_set.curr_epoch*params.time_params.MDP_TIMESTEP
        future_time = curr_time + (state.depth_root+1)*params.time_params.MDP_TIMESTEP

        curr_epoch_dict = context_set.epochs_dict[string(context_set.curr_epoch)]
        epoch_cars = curr_epoch_dict["car-info"]
        # @show state
        car_info = epoch_cars[state.car_id]

        best_value = -Inf

        # Now iterate over car waypoints and find their best value
        route_info = car_info["route"]
        sorted_route = sort(collect(route_info), by=x->parse(Float64, x[1]))

        for (id, timept) in sorted_route

            timeval = timept[2]
            if timeval >= future_time # Only choose points in future

                hopoff_state = get_state_at_rest(US, Point(timept[1][1], timept[1][2]))
                relative_state = rel_state(hopoff_state, dmcts.cmssp.goal_state.continuous)

                val = POMDPs.value(policy, relative_state) - params.cost_params.FLIGHT_REACH_REWARD -
                      params.cost_params.TIME_COEFFICIENT*(timeval - future_time)

                best_value = (val > best_value) ? val : best_value
            end
        end

        # @show best_value
        if best_value == -Inf
            relative_state = rel_state(state.cmssp_state.continuous, dmcts.cmssp.goal_state.continuous)
            best_value = POMDPs.value(policy, relative_state) - dmcts.cmssp.params.cost_params.FLIGHT_REACH_REWARD
        end

        return best_value
    end
end

function init_q_dreamr(policy::P, dmcts::DREAMRMCTSType, state::DREAMRMCTSState{US}, a::DREAMRActionType) where {P <: Policy, US <: UAVState}

    params = dmcts.cmssp.params
    curr_cont_state = state.cmssp_state.continuous

    val = -params.time_params.MDP_TIMESTEP*params.cost_params.TIME_COEFFICIENT

    # If control action, use flight policy action value function for goal state
    if typeof(a.action) <: UAVAction

        new_uav_state = next_state(params, curr_cont_state, a.action)

        val += -dynamics_cost(params, curr_cont_state, new_uav_state)

        relative_state = rel_state(new_uav_state, dmcts.cmssp.goal_state.continuous)
        val += POMDPs.value(policy, relative_state) - dmcts.cmssp.params.cost_params.FLIGHT_REACH_REWARD

    end

    return val

    # else
    #     context_set = dmcts.cmssp.curr_context_set
    #     curr_epoch_dict = context_set.epochs_dict[string(context_set.curr_epoch)]
    #     epoch_cars = curr_epoch_dict["car-info"]

    #     hop_success = false

    #     # If hopon and fail, use flight policy + penalty for missed hop
    #     if a.action.mode_action == HOPON

    #         next_dmcts_state = DREAMRMCTSState(state.cmssp_state, "", state.depth_root+1)

    #         val += -dmcts.cmssp.modeswitch_mdp.R[1, 1]

    #         # val += params.cost_params.INVALID_HOP_PENALTY/2.0

    #         for (car_id, car_info) in epoch_cars

    #             if car_info["route"] != nothing

    #                 car_pos = Point(car_info["pos"][1], car_info["pos"][2])

    #                 uav_pos = get_position(curr_cont_state)
    #                 uav_speed = get_speed(curr_cont_state)

    #                 car_uav_dist = point_dist(car_pos, uav_pos)

    #                 if car_uav_dist < HOP_INFLATION_PARAMETER*params.scale_params.HOP_DISTANCE_THRESHOLD*params.time_params.MDP_TIMESTEP &&
    #                         uav_speed < params.scale_params.XYDOT_HOP_THRESH

    #                     new_uav_state = get_state_at_rest(US, car_pos)
    #                     next_dreamr_state = DREAMRStateType{US}(RIDE, new_uav_state)

    #                     next_dmcts_state = DREAMRMCTSState(next_dreamr_state, car_id, state.depth_root+1)

    #                     hop_success = true

    #                     @show next_dmcts_state
    #                     val += estimate_value_dreamr(policy, dmcts, next_dmcts_state, 0)

    #                     @debug "Successful hop to ", car_id," in sim!"
    #                     break
    #                 end
    #             end
    #         end
    #     elseif a.action.mode_action == HOPOFF

    #         val += -dmcts.cmssp.modeswitch_mdp.R[2, 2]

    #         hopoff_car_id = state.car_id

    #         car_pos = Point(epoch_cars[hopoff_car_id]["pos"][1], epoch_cars[hopoff_car_id]["pos"][2])
    #         new_uav_state = get_state_at_rest(US, car_pos)
    #         next_dreamr_state = DREAMRStateType{US}(FLIGHT, new_uav_state)

    #         next_dmcts_state = DREAMRMCTSState(next_dreamr_state, "", state.depth_root+1)
    #         val += estimate_value_dreamr(policy, dmcts, next_dmcts_state, 0)

    #     else
    #         val += -dmcts.cmssp.modeswitch_mdp.R[2, 3]

    #         car_pos = Point(epoch_cars[state.car_id]["pos"][1], epoch_cars[state.car_id]["pos"][2])

    #         new_uav_state = get_state_at_rest(US, car_pos)

    #         next_dreamr_state = DREAMRStateType{US}(RIDE, new_uav_state)
    #         next_dmcts_state = DREAMRMCTSState(next_dreamr_state, state.car_id, state.depth_root+1)
    #         val += estimate_value_dreamr(policy, dmcts, next_dmcts_state, 0)
    #     end

    #     return val
    # end
end


text:
using HashDicts: HashDict

using BenchmarkTools: @benchmark
import Random

function iterated_entry_insert(d::AbstractDict{K, V}, key_iter, value_iter) where {K, V}
    for (k, v) in zip(key_iter, value_iter)
        d[k] = v
    end
    return nothing
end

function safe_get_functions(d::AbstractDict{K, V}, key_iter, default::V) where {K, V}
    default_counter = 0
    for key in key_iter
        v = get(d, key, default)
        if v == default
            default_counter += 1
        end
    end
    return default_counter
end

function unsafe_getindex(d::AbstractDict, key_iter)
    found = 0
    for key in key_iter
        if haskey(d, key)
            v = d[key]
            found += 1
        end
    end
    return found
end

function basic_benchmark(number_items = 1000, rng::Random.AbstractRNG = Random.GLOBAL_RNG)
    @info "Benchmarking with $number_items items"
    key_iter = rand(rng, Int, number_items)
    value_iter = string.(key_iter)
    for d in (Dict{Int,String}(), HashDict{Int,String}())
        @info "benchmarking $(typeof(d))"
        display(@benchmark iterated_entry_insert($d, $key_iter, $value_iter))
        @info "Benchmarking on existing keys"
        display(@benchmark safe_get_functions($d, $key_iter, ""))
        display(@benchmark unsafe_getindex($d, $key_iter))
        new_keys = rand(rng, Int, number_items)
        @info "Benchmarking on new keys"
        display(@benchmark safe_get_functions($d, $new_keys, ""))
        display(@benchmark unsafe_getindex($d, $new_keys))
    end
end

Random.seed!(42)
basic_benchmark(1000)
basic_benchmark(5000)


text:
# included from Potentials.jl
# part of the module JuLIP.Potentials

using JuLIP: zerovecs, JVecsF, JVecF
using JuLIP.ASE.MatSciPy: NeighbourList

export ZeroPairPotential, PairSitePotential,
         LennardJones, lennardjones,
         Morse, morse

# TODO: why is this not in the abstract file?
# a simplified way to calculate gradients of pair potentials
grad(pp::PairPotential, r::Float64, R::JVec) =
            (evaluate_d(p::PairPotential, r) / r) * R

energy(pp::PairPotential, at::AbstractAtoms) =
         sum( pp(r) for (_1,_2,r,_3,_4) in bonds(at, cutoff(pp)) )

function forces(pp::PairPotential, at::AbstractAtoms)
   dE = zerovecs(length(at))
   for (i,j,r,R,_) in bonds(at, cutoff(pp))
      # TODO: this should be equivalent, but for some reason using @GRAD is much slower!
      # dE[j] -= 2 * @GRAD calc.pp(r, R)   # ∇ϕ(|R|) = (ϕ'(r)/r) R
      dE[j] -= 2 * ((@D pp(r))/r) * R
   end
   return dE
end

"""
`LennardJones:` e0 * ( (r0/r)¹² - 2 (r0/r)⁶ )

Constructor: `LennardJonesPotential(r0, e0)`
"""
LennardJones(r0, e0) =
   AnalyticPotential(:($e0 * (($r0/r)^(12) - 2.0 * ($r0/r)^(6))),
                     id = "LennardJones(r0=$r0, e0=$e0)")

LennardJones() = LennardJones(1.0, 1.0)

"""
`lennardjones(; r0=1.0, e0=1.0, rcut = (1.9*r0, 2.7*r0))`

simplified constructor for `LennardJones` (type unstable)
"""
lennardjones(; r0=1.0, e0=1.0, rcut = (1.9*r0, 2.7*r0)) = (
   (rcut == nothing || rcut == Inf)
         ?  LennardJones(r0, e0)
         :  SplineCutoff(rcut[1], rcut[2]) * LennardJones(r0, e0) )

"""
`Morse(A, e0, r0)`: constructs an
`AnalyticPotential` for
   e0 ( exp( -2 A (r/r0 - 1) ) - 2 exp( - A (r/r0 - 1) ) )
"""
Morse(A, e0, r0) =
   AnalyticPotential(:( $e0 * ( exp(-$(2.0*A) * (r/$r0 - 1.0))
                                - 2.0 * exp(-$A * (r/$r0 - 1.0)) ) ),
                     id="MorsePotential(A=$A, e0=$e0, r0=$r0)")
Morse(;A=4.0, e0=1.0, r0=1.0) = Morse(A, e0, r0)

"""
`morse(A=4.0, e0=1.0, r0=1.0, rcut=(1.9*r0, 2.7*r0))`

simplified constructor for `Morse` (type unstable)
"""
morse(A=4.0, e0=1.0, r0=1.0, rcut=(1.9*r0, 2.7*r0)) = (
   (rcut == nothing || rcut == Inf)
         ?  Morse(A, e0, r0)
         :  SplineCutoff(rcut[1], rcut[2]) * Morse(A, e0, r0) )

@pot type ZeroPairPotential end
"""
`ZeroPairPotential()`: creates a potential that just returns zero
"""
ZeroPairPotential
evaluate(p::ZeroPairPotential, r::Float64) = 0.0
evaluate_d(p::ZeroPairPotential, r::Float64) = 0.0
evaluate_dd(p::ZeroPairPotential, r::Float64) = 0.0
cutoff(p::ZeroPairPotential) = 0.0

# ========================================================
# wrapping a pair potential in a site potential

SitePotential(pp::PairPotential) = PairSitePotential(pp)

Base.zero(::Type{JVecF}) = JVecF(0.0,0.0,0.0)

@pot type PairSitePotential{P} <: SitePotential
   pp::P
end
cutoff(psp::PairSitePotential) = cutoff(psp.pp)
evaluate(psp::PairSitePotential, r, R) =
            sum( [psp.pp(s) for s in r] )
evaluate_d(psp::PairSitePotential, r, R) =
            [ ((@D psp.pp(s))/s) * S for (s, S) in zip(r, R) ]

# ======================================================================
#      Special Preconditioner for Pair Potentials
# ======================================================================


text:
using Parameters
using ModelParameters

const SYMHYD_SOIL_ET_CONST = 10.0

"""
"""
Base.@kwdef mutable struct SYMHYDNode{P, A<:Real} <: NetworkNode
    @network_node

    # parameters
    baseflow_coef::P = Param(0.5, bounds=(0.0, 1.0))
    impervious_threshold::P = Param(2.5, bounds=(0.0, 5.0))  # mm
    infiltration_coef::P = Param(200.0, bounds=(0.0, 400.0))
    infiltration_shape::P = Param(5.0, bounds=(0.0, 10.0))
    interflow_coef::P = Param(0.5, bounds=(0.0, 1.0))
    pervious_fraction::P = Param(0.5, bounds=(0.0, 1.0))
    risc::P = Param(2.5, bounds=(0.0, 5.0))  # rainfall interception store capacity (mm)
    recharge_coef::P = Param(0.5, bounds=(0.0, 1.0))
    smsc::P = Param(250.0, bounds=(1.0, 500.0))  # Soil Moisture Store Capacity (mm)

    # stores
    sm_store::Array{A} = [0.0]
    gw_store::Array{A} = [0.0]
    total_store::Array{A} = [0.0]

    # outputs
    outflow::Array{A} = []  # mm
    baseflow::Array{A} = []  # mm
    quickflow::Array{A} = []  # mm
end

"""

Create node from spec.
"""
function SYMHYDNode(name::String, spec::Dict)
    n = create_node(SYMHYDNode, name, spec["area"])
    node_params = spec["parameters"]
    node_params["sm_store"] = [node_params["initial_sm_store"]]
    node_params["gw_store"] = [node_params["initial_gw_store"]]
    node_params["total_store"] = [node_params["initial_total_store"]]

    remove_from_spec = (x) -> delete!(node_params, x)
    map(remove_from_spec, ["initial_sm_store", "initial_gw_store", "initial_total_store"])

    for (k, p) in node_params
        s = Symbol(k)
        if p isa String
            p = eval(Meta.parse(p))
        end

        try
            f = getfield(n, s)
            setfield!(n, s, Param(p, bounds=f.bounds))
        catch err
            msg = sprint(showerror, err, catch_backtrace())
            if occursin("no field bounds", string(msg))
                setfield!(n, s, p)
            else
                throw(err)
            end
        end
    end

    return n
end

"""

Run SYMHYD for all time steps within a climate scenario.
"""
function run_node!(node::SYMHYDNode, climate::Climate; inflow=nothing, extraction=nothing, exchange=nothing)
    timesteps = length(climate)
    for ts in 1:timesteps
        run_timestep!(node, climate, ts; inflow=inflow, extraction=extraction, exchange=exchange)
    end

    return node.outflow
end

"""

Run SYMHYD for a given time step
"""
function run_timestep!(node::SYMHYDNode, climate::Climate, timestep::Int; inflow=nothing, extraction=extraction, exchange=nothing)
    P, E = climate_values(node, climate, timestep)

    run_timestep!(node, P, E, timestep; inflow=inflow, extraction=extraction, exchange=exchange)
end

function run_timestep!(node::SYMHYDNode, rain, et, ts; inflow=nothing, extraction=nothing, exchange=nothing)
    res = run_symhyd(node, rain, et)
    sm_store, gw_store, total_store, total_runoff, baseflow, event_runoff = res

    node_name = node.name
    wo = timestep_value(ts, node_name, "releases", extraction)
    ex = timestep_value(ts, node_name, "exchange", exchange)
    in_flow = timestep_value(ts, node_name, "inflow", inflow)

    if !isnothing(inflow)
        total_runoff = total_runoff + in_flow + ex - wo
    end

    area = node.area
    update_state!(node, sm_store, gw_store, total_store, total_runoff * area, baseflow * area, event_runoff * area)
end

function update_state!(node::SYMHYDNode, sm_store, gw_store, total_store, outflow, baseflow, quickflow)
    push!(node.sm_store, sm_store)
    push!(node.gw_store, gw_store)
    push!(node.total_store, total_store)
    push!(node.outflow, outflow)
    push!(node.baseflow, baseflow)
    push!(node.quickflow, quickflow)
end

function update_params!(node::SYMHYDNode, baseflow_coef::Float64, impervious_threshold::Float64,
                        infiltration_coef::Float64,
                        infiltration_shape::Float64,
                        interflow_coef::Float64,
                        pervious_fraction::Float64,
                        risc::Float64,
                        recharge_coef::Float64,
                        smsc::Float64)::Nothing
    node.baseflow_coef = Param(baseflow_coef, bounds=node.baseflow_coef.bounds)
    node.impervious_threshold = Param(impervious_threshold, bounds=node.impervious_threshold.bounds)
    node.infiltration_coef = Param(infiltration_coef, bounds=node.infiltration_coef.bounds)
    node.infiltration_shape = Param(infiltration_shape, bounds=node.infiltration_shape.bounds)
    node.interflow_coef = Param(interflow_coef, bounds=node.interflow_coef.bounds)
    node.pervious_fraction = Param(pervious_fraction, bounds=node.pervious_fraction.bounds)
    node.risc = Param(risc, bounds=node.risc.bounds)
    node.recharge_coef = Param(recharge_coef, bounds=node.recharge_coef.bounds)
    node.smsc = Param(smsc, bounds=node.smsc.bounds)

    return nothing
end

"""
    reset!(node::SYMHYDNode)::Nothing

Reset node. Clears all states back to their initial values.
"""
function reset!(node::SYMHYDNode)::Nothing
    # stores
    node.sm_store = [node.sm_store[1]]
    node.gw_store = [node.gw_store[1]]
    node.total_store = [node.total_store[1]]

    # outputs
    node.outflow = []
    node.baseflow = []
    node.quickflow = []

    return nothing
end

"""

Run SYMHYD for a single time step with given inputs and state variables.
"""
function run_symhyd(node::SYMHYDNode, P, ET)

    sm_store = node.sm_store[end]
    gw_store = node.gw_store[end]
    total_store = node.total_store[end]

    pervious_incident = P
	impervious_incident = P

    impervious_ET = min(node.impervious_threshold, impervious_incident)
    impervious_runoff = impervious_incident - impervious_ET

    interception_ET = min(pervious_incident, min(ET, node.risc))
    throughfall = pervious_incident - impervious_ET

    smsc = node.smsc
    sm_fraction = sm_store / smsc

    infiltration_capacity = node.infiltration_coef * exp(-node.infiltration_shape*sm_fraction)
    infiltration = min(throughfall, infiltration_capacity)
    infiltration_Xs_runoff = throughfall - infiltration

    interflow_runoff = node.interflow_coef * sm_fraction * infiltration
    infiltration_after_interflow = infiltration - interflow_runoff

    recharge = node.recharge_coef * sm_fraction * infiltration_after_interflow
    soil_input = infiltration_after_interflow - recharge

    # Update states...
    sm_store += soil_input
    sm_fraction = sm_store / smsc
    gw_store += recharge

    if sm_fraction > 1.0
        gw_store += sm_store - smsc
        sm_store = smsc
        sm_fraction = 1.0
    end

    baseflow_runoff = node.baseflow_coef * gw_store
    gw_store -= baseflow_runoff

    soil_ET = min(sm_store, min(ET - interception_ET, sm_fraction*SYMHYD_SOIL_ET_CONST))
    sm_store -= soil_ET

    pervious_frac = node.pervious_fraction
    total_store = (sm_store + gw_store) * pervious_frac

    event_runoff = (1.0 - pervious_frac) * impervious_runoff + pervious_frac * (infiltration_Xs_runoff + interflow_runoff)
    total_runoff = event_runoff + pervious_frac * baseflow_runoff
    baseflow = baseflow_runoff * pervious_frac

	# values for time step...
	return sm_store, gw_store, total_store, total_runoff, baseflow, event_runoff
end


text:
module BKMaxflow

include("bitmask.jl")
export BKStatusBits, BK_EMPTY, BK_FREE, BK_S, BK_T, BK_ACTIVE, BK_ORPHAN
export BK_S_ACTIVE, BK_T_ACTIVE, BK_S_ORPHAN, BK_T_ORPHAN, BK_S_ACTIVE_ORPHAN, BK_T_ACTIVE_ORPHAN

include("boykov_kolmogorov.jl")
export boykov_kolmogorov

abstract type AbstractBKMaxflow{Tv,Ti} end
struct JuliaImpl{Tv<:Real,Ti<:Integer} <: AbstractBKMaxflow{Tv,Ti} end
struct CImpl{Tv<:Real} <: AbstractBKMaxflow{Tv,Integer} end

function create_graph(::Type{JuliaImpl{Tv,Ti}}, vertexNum::Integer) where {Tv<:Real,Ti<:Integer}
    weights = Vector{Tv}()
    neighbors = [Vector{Tuple{Ti,Ti}}() for i = 1:vertexNum]
    return weights, neighbors
end

function add_edge!(weights::AbstractVector{Tv}, neighbors::AbstractVector{Vector{Tuple{Ti,Ti}}},
                   p::Ti, q::Ti, weightp2q::Tv, weightq2p::Tv) where {Tv<:Real,Ti<:Integer}
    idx = length(weights) ÷ 2
    # store forward and reverse weights compactly for cache efficiency,
    # assume a forward weight is the weight from a smaller index to a bigger index
    if p < q
        push!(weights, weightp2q, weightq2p)
        push!(neighbors[p], (q,idx+1))
        push!(neighbors[q], (p,idx+1))
    else
        push!(weights, weightq2p, weightp2q)
        push!(neighbors[q], (p,idx+1))
        push!(neighbors[p], (q,idx+1))
    end
end

export AbstractBKMaxflow, JuliaImpl, CImpl
export create_graph, add_edge!

include("c_wrapper.jl")

end # module

max_stars_repo_name: Gnimuc/BKmaxflow

max_issues_repo_name: Gnimuc/BKmaxflow

max_forks_repo_name: Gnimuc/BKmaxflow


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENCE in the project root.
# ------------------------------------------------------------------

"""
    State

A geological state to be evolved by processes.
"""
abstract type State end

"""
    LandState(land)

A geological state consisting of a landscape map.
"""
struct LandState <: State
  land::Matrix{Float64}
end

@recipe function f(state::LandState)
  aspect_ratio --> :equal
  colorbar --> :false
  seriestype --> :surface
  state.land
end


text:
using Documenter, DoubleFloats

makedocs(
    modules = [DoubleFloats],
    sitename = "DoubleFloats.jl",
    authors = "Jeffrey Sarnoff, Simon Byrne, Sascha Timme, and other contributors",
    pages = Any[
        "Overview" => "index.md",
        "Construction" => "construction.md",
        "Characteristics" => "characteristics.md",
        "Capabilities" => "capabilities.md",
        "Appropriate Uses" => "appropriate.md",
        "Refs" => "references.md"
    ]
)

deploydocs(
    repo = "github.com/JuliaMath/DoubleFloats.jl.git",
    target = "build"
)


text:
### teste 14

function teste(cpf)
    um = 1
    n = 0
    dez = 10
    onze = dez
    zero = 0
    p = zero
    q = p
    m = q

    while onze >= zero
        base = dez^onze
        m = div(cpf, base)
        if onze >= 6
            p = p + m
            elseif onze == 5
                if m == 0
                    n = p
                else
                    n = m
                end
        else
            q = q + m
        end
     cpf = div(cpf, base)
     onze = onze - um

    end

    p = div(p, n)
    q = div(p, n)
    saida = p - q
    return saida

end

teste(57867564534)

x = range(99999991111, 99999999999)
x = teste.(x)

plot(x)


text:
# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
# Since these FMAs can increase the performance of many numerical algorithms,
# we need to opt-in explicitly.
# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
@muladd begin

# du .= zero(eltype(du)) doesn't scale when using multiple threads.
# See https://github.com/trixi-framework/Trixi.jl/pull/924 for a performance comparison.
function reset_du!(du, dg, cache)
  @threaded for element in eachelement(dg, cache)
      du[.., element] .= zero(eltype(du))
  end

  return du
end

#     pure_and_blended_element_ids!(element_ids_dg, element_ids_dgfv, alpha, dg, cache)
#
# Given blending factors `alpha` and the solver `dg`, fill
# `element_ids_dg` with the IDs of elements using a pure DG scheme and
# `element_ids_dgfv` with the IDs of elements using a blended DG-FV scheme.
function pure_and_blended_element_ids!(element_ids_dg, element_ids_dgfv, alpha, dg::DG, cache)
  empty!(element_ids_dg)
  empty!(element_ids_dgfv)

  for element in eachelement(dg, cache)
    # Clip blending factor for values close to zero (-> pure DG)
    dg_only = isapprox(alpha[element], 0, atol=1e-12)
    if dg_only
      push!(element_ids_dg, element)
    else
      push!(element_ids_dgfv, element)
    end
  end

  return nothing
end

function volume_jacobian(element, mesh::TreeMesh, cache)
  return inv(cache.elements.inverse_jacobian[element])^ndims(mesh)
end

# Indicators used for shock-capturing and AMR
include("indicators.jl")
include("indicators_1d.jl")
include("indicators_2d.jl")
include("indicators_3d.jl")

# Container data structures
include("containers.jl")

# 1D DG implementation
include("dg_1d.jl")

# 2D DG implementation
include("dg_2d.jl")
include("dg_2d_parallel.jl")

# 3D DG implementation
include("dg_3d.jl")

# Auxiliary functions that are specialized on this solver
include("dg_2d_compressible_euler.jl")

end # @muladd


text:
using HWconstrained
using Test

@testset "HWconstrained.jl" begin
	@testset "testing components" begin

		@testset "tests gradient of objective function" begin

		end

		@testset "tests gradient of constraint function" begin

		end
	end

	@testset "testing result of both maximization methods" begin

		truth = HWconstrained.DataFrame(a=[0.5;1.0;5.0],
			                            c = [1.008;1.004;1.0008],
			                            omega1=[-1.41237;-0.20618;0.758763],
			                            omega2=[0.801455;0.400728;0.0801455],
			                            omega3=[1.60291;0.801455;0.160291],
			                            fval=[-1.20821;-0.732819;-0.013422])

		@testset "checking result of NLopt maximization" begin

			t1 = table_NLopt()
			for c in names(truth)
				@test all(maximum.(abs.(t1[c].-truth[c])) .< tol2)
			end
		end

		@testset "checking result of NLopt maximization" begin
			t1 = table_JuMP()
			for c in names(truth)
				@test all(maximum.(abs.(t1[c].-truth[c])) .< tol2)
			end
		end
	end

end


text:
using Magrad
using Test

include("test_tensorops.jl")

max_stars_repo_name: leonvol/magrad

max_issues_repo_name: leonvol/magrad

max_forks_repo_name: leonvol/magrad


text:
using DataFrames
using Plots
using StatsBase
activity_symbols=(:sleep,:work,:tv,:edu,:household,:relax,:carehh,:eat,:travel,:personalcare,:sport,:buy,:religion,:telephone,:carenhh,:volunteer)

include("make_dataset.jl")
include("analyze.jl")
include("visualize.jl")

function main(;recreate=false)
  if !isfile("../data/processed/processed.csv") || recreate
    print("Pre-processing data files... ")
    preprocess_datafiles()
    print("done. ")
  end

  print("Reading data file... ")
  data=readtable("../data/processed/processed.csv", nrows=-1)
  println("done")

  print("Analyzing... ")
  work_df, nwork_df =aggregate_by_age(data)
  println("done")

  print("Plotting... ")
  create_plot_workvsfree(work_df,nwork_df)
  create_plot_diff(work_df,nwork_df)
  println("done")
end

#cd to src directory
cd(dirname(@__FILE__))

recreate=false
if(length(ARGS) > 0)
  if(ARGS[1] == "--recreate")
    recreate=true
  end
end

main(recreate=recreate)


text:
module Pathfinder

using Distributions: Distributions
using LinearAlgebra
using Optim: Optim, LineSearches
using PDMats: PDMats
using PSIS: PSIS
using Random
using Statistics: Statistics
using StatsBase: StatsBase
using StatsFuns: log2π
using Nonconvex, NonconvexCore, NonconvexIpopt, ChainRulesCore

export pathfinder, multipathfinder

# Note: we override the default history length to be shorter and the default line search
# to be More-Thuente, which keeps the approximate inverse Hessian positive-definite
const DEFAULT_HISTORY_LENGTH = 6
const DEFAULT_LINE_SEARCH = LineSearches.MoreThuente()
# const DEFAULT_OPTIMIZER = Optim.LBFGS(;
#     m=DEFAULT_HISTORY_LENGTH, linesearch=DEFAULT_LINE_SEARCH
# )
const DEFAULT_OPTIMIZER = NonconvexIpopt.IpoptAlg()

include("woodbury.jl")
include("maximize.jl")
include("inverse_hessian.jl")
include("mvnormal.jl")
include("elbo.jl")
include("resample.jl")
include("singlepath.jl")
include("multipath.jl")

end


text:
function HansenZonalLinear(arg0::jint, arg1::jint)
    return HansenZonalLinear((jint, jint), arg0, arg1)
end

function compute_init_values(obj::HansenZonalLinear, arg0::jdouble)
    return jcall(obj, "computeInitValues", void, (jdouble,), arg0)
end

function get_derivative(obj::HansenZonalLinear, arg0::jint, arg1::jdouble)
    return jcall(obj, "getDerivative", jdouble, (jint, jdouble), arg0, arg1)
end

function get_value(obj::HansenZonalLinear, arg0::jint, arg1::jdouble)
    return jcall(obj, "getValue", jdouble, (jint, jdouble), arg0, arg1)
end


text:
using Test
using Actors
import Actors: hear

# We really want to test with true concurrency. Even if we are only running in
# a VM with one vCPU (i.e. Gitlab CI), we are more likely to see errors with
# atleast two OS threads. Unfortunately Julia limits the number of threads to
# the CPU count though.
if Sys.CPU_THREADS > 1
    @assert Threads.nthreads() > 1 "Set environment variable JULIA_NUM_THREADS to > 1"
else
    @warn "Only testing with one thread"
end

include("TestSet.jl")
include("addressing.jl")

# Sanity check without using our custom TestSet
"Our Play"
struct HelloWorld
    chnl::Channel
end

"Our Actor"
struct Julia end

"Our Message"
struct HelloWorld!
    chnl::Channel
end

function hear(s::Scene{A}, msg::HelloWorld!) where A
    put!(msg.chnl, "Hello, World! I am $(A)!")
    say(s, stage(s), Leave!())
end

function hear(s::Scene{HelloWorld}, ::Genesis!)
    julia = invite!(s, Julia())
    say(s, julia, HelloWorld!(my(s).chnl))
end

@testset "Hello, World!" begin
    p = HelloWorld(Channel(1))

    play!(p)
    @test take!(p.chnl) == "Hello, World! I am Julia!"
    close(p.chnl)
end

include("Minders.jl")
include("Logger.jl")

struct TestSetTest end

function hear(s::Scene{TestSetTest}, ::Genesis!)
    @assert Test.get_testset() isa LuvvyTestSet
    @test true

    say(s, stage(s), Leave!())
end

@testset LuvvyTestSet expect=1 "TestSet Test" begin
    play!(TestSetTest())
end

include("Async.jl")
include("Luvvies.jl")
include("Stack.jl")
include("TypedMessages.jl")


text:
function _outdegree(x::TN) where {TN <: Union{TNDistribution, TransmissionNetwork}}
  infected = [x.external[i] > 0 || any(x.internal[:, i] .> 0) for i = 1:individuals(x)]
  return sum(x.internal[infected, infected], dims=2)
end

@recipe function plot(x::TN) where {TN <: Union{TNDistribution, TransmissionNetwork}}
  seriestype := :histogram
  legend --> :none
  xlab --> "Out degree"
  ylab --> "Frequency"
  _outdegree(x)
end


text:
# Autogenerated wrapper script for UASM_jll for x86_64-linux-musl
export uasm

JLLWrappers.@generate_wrapper_header("UASM")
JLLWrappers.@declare_executable_product(uasm)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_executable_product(
        uasm,
        "bin/uasm",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# Tally and Accumulate

using SimLynx
SimLynx.greet()

tallied = nothing
accumulated = nothing

@process test_process(value_durations) begin
    for (value, duration) in value_durations
        tallied.value = value
        accumulated.value = value
        work(duration)
    end
end

function main(value_durations)
    println("--- Test Tally and Accumulate ---")
    @simulation begin
        # Create tallied and accumulated variables
        global tallied = Variable{Int64}(data=:tally, history=true)
        global accumulated = Variable{Int64}(0, history=true)
        # Schedule the test process and start the simulation
        @schedule at 0.0 test_process(value_durations)
        start_simulation()
        # Print and plot the tallied results
        println("--- Tally ---")
        print_stats(tallied, title="Tallied Statistics")
        plot_history(tallied, title="Tallied History")
        # Print and plot the accumulated results
        println("--- Accumulate ---")
        print_stats(accumulated, title="Accumulated Statistics")
        plot_history(accumulated, title="Accumulated History")
    end
end

main([(1, 2.0), (2, 1.0), (3, 2.0), (4, 3.0)])


text:
"""
Fancy way of updating IntVert. We need to add the vertices that are common and contained.

"""
function PolytopeGeneratingVertices(Simplex1::AbstractArray{T, 2},
                                    Simplex2::AbstractArray{T, 2},
                                    IntVert,#::AbstractArray{T, 2},
                                    ConvexExpIntVert,#::AbstractArray{T, 2},
                                    convexexp1in2,#::AbstractArray{T, 2},
                                    convexexp2in1, #::AbstractArray{T, 2},
                                    ordered_vertices1::AbstractVector{Int},
                                    ordered_vertices2::AbstractVector{Int},
                                    numof1in2::Int,
                                    numof2in1::Int,
                                    Ncomm::Int) where T
    n = size(Simplex1, 1)

    numof1in2NotShared = numof1in2 - Ncomm
    numof2in1NotShared = numof2in1 - Ncomm

    # Deal first with the possibly shared vertices
    if Ncomm > 0
        #println("Deal first with the possibly shared vertices")
        IndexSharedVert1 = round.(Int64, ordered_vertices1[1:Ncomm])
        IndexSharedVert2 = round.(Int64, ordered_vertices2[1:Ncomm])
        IntVert = vcat(
            IntVert,
            transpose(Simplex1[:, IndexSharedVert1])
        )

        common_vertices_expressed_in1 = transpose(convexexp2in1[:, IndexSharedVert2])
        common_vertices_expressed_in2 = transpose(convexexp1in2[:, IndexSharedVert1])

        ConvexExpIntVert = vcat(
            ConvexExpIntVert,
            [common_vertices_expressed_in1 common_vertices_expressed_in2]
        )
    end

    # Dealing with vertices that are contained but not shared

    if numof1in2NotShared * numof2in1NotShared > 0
        # Each simplex contains vertices of the other beyond the shared ones
        #println("# Each simplex contains vertices of the other beyond the shared ones")
        #Id = eye(n+1, n+1)
		Id = Matrix{Float64}(I, n + 1, n + 1)

        indexvert1in2notshared = ordered_vertices1[Ncomm+1:numof1in2]
        indexvert2in1notshared = ordered_vertices2[Ncomm+1:numof2in1]

        IntVert = vcat(
            IntVert,
            copy(transpose(Simplex1[:, indexvert1in2notshared])),
            copy(transpose(Simplex2[:, indexvert2in1notshared]))
        )

        conv_exp_vertof1_expressed_in1 = Id[indexvert1in2notshared, :]
        conv_exp_vertof2_expressed_in2 = Id[indexvert2in1notshared, :]
        conv_exp_vertof1_expressed_in2 = transpose(convexexp1in2[:, indexvert1in2notshared])
        conv_exp_vertof2_expressed_in1 = transpose(convexexp2in1[:, indexvert2in1notshared])

        ConvexExpIntVert = vcat(
            ConvexExpIntVert,
            [conv_exp_vertof1_expressed_in1 conv_exp_vertof1_expressed_in2],
            [conv_exp_vertof2_expressed_in1 conv_exp_vertof2_expressed_in2]
        )

    elseif numof1in2NotShared > 0
        # Only Simplex2 contains vertices of Simplex1 other than the shared ones
        #println("Only Simplex2 contains vertices of Simplex1 other than the shared ones")
        Id = Matrix(1.0I, n+1, n+1)

        indexvert1in2notshared = ordered_vertices1[Ncomm+1:numof1in2]

        IntVert = vcat(
            IntVert,
            transpose(Simplex1[:, indexvert1in2notshared])
        )

        conv_exp_vertof1_expressed_in1 = Id[indexvert1in2notshared, :]
        conv_exp_vertof1_expressed_in2 = transpose(convexexp1in2[:, indexvert1in2notshared])

        ConvexExpIntVert = vcat(
            ConvexExpIntVert,
            [conv_exp_vertof1_expressed_in1 conv_exp_vertof1_expressed_in2]
        )

    elseif numof2in1NotShared > 0
        # Only Simplex1 contains vertices of Simplex2 other than the shared ones
        #println("Only Simplex1 contains vertices of Simplex2 other than the shared ones")

        Id = Matrix(1.0I, n + 1, n + 1)

        indexvert2in1notshared = ordered_vertices2[Ncomm+1:numof2in1]
        IntVert = vcat(
            IntVert,
            transpose(Simplex2[:, indexvert2in1notshared])
        )

        conv_exp_vertof2_expressed_in2 = Id[indexvert2in1notshared, :]
        conv_exp_vertof2_expressed_in1 = transpose(convexexp2in1[:, indexvert2in1notshared])
        ConvexExpIntVert = vcat(ConvexExpIntVert,
                                    [conv_exp_vertof2_expressed_in1 conv_exp_vertof2_expressed_in2])
    end

    return IntVert, ConvexExpIntVert
end


text:
using BinDeps
@BinDeps.setup

opencv_core = library_dependency("opencv_core", aliases = ["libopencv_core","libopencv_core.3.0.0"],
    validate = function(p,h)
        # Check that we don't have opencv version 2
        # Not actually sure if this is a valid way to check this,
        # but it happened to be true of the opecv I have on my system
        Libdl.dlsym_e(h,"cvLine") == C_NULL
    end)
opencv_highgui = library_dependency("opencv_highgui", aliases = ["libopencv_highgui","libopencv_highgui.3.0.0"])
opencv_improc = library_dependency("opencv_imgproc", aliases = ["libopencv_imgproc","libopencv_imgproc.3.0.0"])
opencv_objdetect = library_dependency("opencv_objdetect", aliases = ["libopencv_objdetect","libopencv_objdetect.3.0.0"])
opencv_text = library_dependency("opencv_text", aliases = ["libopencv_text","opencv_text.3.0.0"])

@osx_only begin
  using Homebrew
  provides( Homebrew.HB, "opencv3", opencv_core, os = :Darwin )
end

@BinDeps.install Dict(:opencv_core => :opencv_core, :opencv_highgui => :opencv_highgui,
                      :opencv_imgproc => :opencv_imgproc, :opencv_objdetect => :opencv_objdetect,
                      :opencv_text => :opencv_text)


text:
import Base.BLAS: gemm!, gemv!, syr!

function takelast!{F <: AbstractFloat, T <: Integer}(
    dest::Array{F}, src::Array{F},
    idxes::Union{Tuple{Vararg{T}}, Array{T}})
    pre = ntuple(i->Colon(), ndims(src)-1)
    for (i, k) in enumerate(idxes)
        dest[pre..., i] = src[pre..., k]
    end
    dest
end

takelast!(dest, src, idx::Integer) = takelast!(dest, src, (idx,))

function viewlast{F <: AbstractFloat}(src::Array{F}, idx::Integer)
    pre = ntuple(i->Colon(), ndims(src)-1)
    @view src[pre..., idx]
end

function splitarray(a::Vector, shape1::Tuple, shape2::Tuple)
    l1 = prod(shape1)
    l2 = prod(shape2)
    @assert length(a) == l1 + l2
    p1 = @view a[1:l1]
    p2 = @view a[l1+1:end]
    return (reshape(p1, shape1), reshape(p2, shape2))
end

if BACKEND == "GPU"
    include("gpu_array.jl")
end


text:
using Statistics
using NCDatasets

using ArgParse
using JSON
using Formatting

include("LinearRegression.jl")
include("nanop.jl")
include("CESMReader.jl")

using .CESMReader

function parse_commandline()

    s = ArgParseSettings()
    @add_arg_table s begin

        "--data-file-prefix"
            help = "Data filename prefix including folder and path until the timestamp. File extension `nc` is assumed."
            arg_type = String
            required = true

        "--data-file-timestamp-form"
            help = "Data filename timestamp form. Either `YEAR` or `YEAR_MONTH`."
            arg_type = String
            required = true

        "--domain-file"
            help = "Domain file."
            arg_type = String
            required = true

        "--output-file"
            help = "Output file."
            arg_type = String
            required = true

        "--beg-year"
            help = "Year of begin."
            arg_type = Int64
            required = true

        "--end-year"
            help = "Year of end."
            arg_type = Int64
            required = true

        "--varname"
            help = "Variable name. If it is 2D then it should be (lon, lat, time), if it is 3D then it should be (lon, lat, lev/depth, time)"
            arg_type = String
            required = true

        "--dims"
            help = "The form of dimensions. Now can be `XYZT`, `XYT`, `YZT`"
            arg_type = String
            required = true

    end

    return parse_args(ARGS, s)
end

parsed = parse_commandline()
print(json(parsed, 4))

output_file = parsed["output-file"]

if parsed["data-file-timestamp-form"] == "YEAR"
    filename_format = format("{:s}{{:04d}}.nc", joinpath(parsed["data-file-prefix"]))
    form = :YEAR
elseif parsed["data-file-timestamp-form"] == "YEAR_MONTH"
    filename_format = format("{:s}{{:04d}}-{{:02d}}.nc", joinpath(parsed["data-file-prefix"]))
    form = :YEAR_MONTH
end

Dataset(parsed["domain-file"], "r") do ds
    global mask = replace(ds["mask"], missing=>NaN)
end

if ! ( parsed["dims"] in ["XYZT", "XYT", "YZT", "YZ"] )
    ErrorException("Unknown dimension: " * parsed["dims"]) |> throw
end

# Creating preknowledge of dimensions
Dataset(format(filename_format, parsed["beg-year"], 1), "r") do ds

    var = ds[parsed["varname"]]
    dims = size(var)

    global beg_t = (parsed["beg-year"] - 1) * 12 + 1
    global end_t = (parsed["end-year"] - 1) * 12 + 12

    if parsed["dims"] == "XYZT"
        global (Nx, Ny, Nz, _ ) = size(var)
    elseif parsed["dims"] == "XYT"
        global (Nx, Ny, _ ) = size(var)
        Nz = 1
    elseif parsed["dims"] == "YZT"
        global (Ny, Nz, _ ) = size(var)
        Nx = 1
    elseif parsed["dims"] == "YZ"
        global (Ny, Nz ) = size(var)
        Nx = 1
    end

    global Nt = end_t - beg_t + 1

    if mod(Nt, 12) != 0
        ErrorException("Time record is not multiple of 12") |> throw
    end

    global Na = Int64(Nt / 12)
    global Ns = Int64(Nt / 3)
end

months  = collect(Float64, 1:Nt)
years   = collect(Float64, 1:Na)

fh = FileHandler(filename_format=filename_format, form=form)

if parsed["dims"] == "XYZT"

    spatial_rng = (:, :, :)
    global data = reshape( getData(fh, parsed["varname"], (parsed["beg-year"], parsed["end-year"]), spatial_rng), Nx, Ny, Nz, Nt)

elseif  parsed["dims"] == "XYT"

    spatial_rng = (:, :)
    global data = reshape( getData(fh, parsed["varname"], (parsed["beg-year"], parsed["end-year"]), spatial_rng), Nx, Ny, Nz, Nt)

end

data_MM = nanmean( data, dims=(1,) )

data_AM = zeros(Float64, 1, Ny, Nz, Na)

for a=1:Na
    data_AM[:, :, :, a] = mean( view(data_MM, :, :, :, ((a-1)*12+1):(a*12) ), dims=4 )[:, :, :, 1]
end

Dataset(output_file, "c") do ds

    defDim(ds, "years",  Na)
    defDim(ds, "months", Nt)
    defDim(ds, "Nx", 1)
    defDim(ds, "Ny", Ny)
    defDim(ds, "Nz", Nz)

    datas =  convert(Array{Any}, [

        (format("{:s}_MM",    parsed["varname"]),       data_MM,    ("Nx", "Ny", "Nz", "months"), Dict()),
        (format("{:s}_AM",    parsed["varname"]),       data_AM,    ("Nx", "Ny", "Nz", "years"), Dict()),

    ])

    if  parsed["dims"] in [ "XYT", "XYZT" ]
        push!(datas, ("mask",  mask,    ("Nx", "Ny"),   Dict()))
    end

    for (varname, vardata, vardim, attrib) in datas
        if ! haskey(ds, varname)
            var = defVar(ds, varname, Float64, vardim)
            var.attrib["_FillValue"] = 1e20
        end

        println("Writing variable:  ", varname)
        var = ds[varname]

        for (k, v) in attrib
            var.attrib[k] = v
        end

        rng = []
        for i in 1:length(vardim)-1
            push!(rng, Colon())
        end
        push!(rng, 1:size(vardata)[end])
        var[rng...] = vardata

    end

end

max_stars_repo_name: meteorologytoday/CESM-diagnostic

max_issues_repo_name: meteorologytoday/CESM-diagnostic

max_forks_repo_name: meteorologytoday/CESM-diagnostic


text:
using Cxx
using Libdl

#const path_to_lib = "/home/jarv1s/pixhawk_sensor_data"
#const path_to_header = joinpath(@__DIR__, "..", "deps", "usr", "include")

#addHeaderDir(path_to_lib, kind=C_System)
#Libdl.dlopen(path_to_lib * "/libPixhawkData.so", Libdl.RTLD_GLOBAL)

cxx"#include <iostream>"
cxx"#include <iostream>"
cxx"#include <stdio.h>"
cxx"#include <cstdlib>"
cxx"#include <unistd.h>"
cxx"#include <cmath>"
cxx"#include <string.h>"
cxx"#include <inttypes.h>"
cxx"#include <fstream>"
cxx"#include <signal.h>"
cxx"#include <time.h>"
cxx"#include <sys/time.h>"

cxx"#include <mavlink/include/mavlink/v2.0/common/mavlink.h>"

cxxinclude("autopilot_interface.h")
cxxinclude("serial_port.h")
cxxinclude("mavlink_control.h")

cxx"""
    int main(int argc, char **argv);
    int top(int argc, char **argv);

    void commands(Autopilot_Interface &autopilot_interface);
    void parse_commandline(int argc, char **argv, char *&uart_name, int &baudrate);

    void quit_handler( int sig );
"""

cxx"""
int
top (int argc, char **argv)
{
#ifdef __APPLE__
	char *uart_name = (char*)"/dev/tty.usbmodem1";
#else
	char *uart_name = (char*)"/dev/ttyUSB0";
#endif
	int baudrate = 57600;

	parse_commandline(argc, argv, uart_name, baudrate);

	Serial_Port serial_port(uart_name, baudrate);

	Autopilot_Interface autopilot_interface(&serial_port);

	serial_port_quit         = &serial_port;
	autopilot_interface_quit = &autopilot_interface;
	signal(SIGINT,quit_handler);

	serial_port.start();
	autopilot_interface.start();

	commands(autopilot_interface);

	autopilot_interface.stop();
	serial_port.stop();

	return 0;
}
"""

top() = @cxx top()

cxx"""
void
commands(Autopilot_Interface &api)
{
	Mavlink_Messages messages = api.current_messages;

	mavlink_highres_imu_t imu = messages.highres_imu;

	printf("acc_x =  % f, acc_y = %f, acc_z = %f \n", imu.xacc, imu.yacc, imu.zacc); // (m/s^2)
	printf("gyro_x = % f, gyro_y = %f, gyro_z = %f \n", imu.xgyro, imu.ygyro, imu.zgyro); // (rad/s)
	printf("mag_x = % f, mag_y = %f, mag_z = %f \n", imu.xmag, imu.ymag, imu.zmag); // (Ga)

	return;

}
"""

commands() = @cxx commands()

cxx"""
void
parse_commandline(int argc, char **argv, char *&uart_name, int &baudrate)
{

	// string for command line usage
	const char *commandline_usage = "usage: mavlink_serial -d <devicename> -b <baudrate>";

	// Read input arguments
	for (int i = 1; i < argc; i++) { // argv[0] is "mavlink"

		// Help
		if (strcmp(argv[i], "-h") == 0 || strcmp(argv[i], "--help") == 0) {
			printf("%s\n",commandline_usage);
			throw EXIT_FAILURE;
		}

		// UART device ID
		if (strcmp(argv[i], "-d") == 0 || strcmp(argv[i], "--device") == 0) {
			if (argc > i + 1) {
				uart_name = argv[i + 1];

			} else {
				printf("%s\n",commandline_usage);
				throw EXIT_FAILURE;
			}
		}

		// Baud rate
		if (strcmp(argv[i], "-b") == 0 || strcmp(argv[i], "--baud") == 0) {
			if (argc > i + 1) {
				baudrate = atoi(argv[i + 1]);

			} else {
				printf("%s\n",commandline_usage);
				throw EXIT_FAILURE;
			}
		}

	}

	return;
}
"""

parse_commandline() = @cxx parse_commandline()

cxx"""
void
quit_handler( int sig )
{
	printf("\n");
	printf("TERMINATING AT USER REQUEST\n");
	printf("\n");

	// autopilot interface
	try {
		autopilot_interface_quit->handle_quit(sig);
	}
	catch (int error){}

	// serial port
	try {
		serial_port_quit->handle_quit(sig);
	}
	catch (int error){}

	// end program here
	exit(0);

}
"""

quit_handler() = @cxx quit_handler()

cxx"""
int
main(int argc, char **argv)
{
	// This program uses throw, wrap one big try/catch here
	try
	{
		int result = top(argc,argv);
		return result;
	}

	catch ( int error )
	{
		fprintf(stderr,"mavlink_control threw exception %i \n" , error);
		return error;
	}

}
"""

main() = @cxx main()


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

abstract type BATExecutor end

default_executor() = MultiThreadedExec()

struct SequentialExec <: BATExecutor end

function exec_map!(f::Base.Callable, executor::SequentialExec, Y::AbstractVector, X::AbstractVector)
    @argcheck length(eachindex(X)) == length(eachindex(X))
    for i in 0:(length(eachindex(Y)) - 1)
        Y[firstindex(Y) + i] = f(X[firstindex(X) + i])
    end
    return Y
end

struct MultiThreadedExec <: BATExecutor end

function exec_map!(f::Base.Callable, executor::MultiThreadedExec, Y::AbstractVector, X::AbstractVector)
    @argcheck length(eachindex(X)) == length(eachindex(X))
    @threads for i in 0:(length(eachindex(Y)) - 1)
        Y[firstindex(Y) + i] = f(X[firstindex(X) + i])
    end
    return Y
end

@with_kw struct DistributedExec{WP<:AbstractWorkerPool} <: BATExecutor
    workers::WP = WorkerPool(workers())
    batchsize::Int = 1
end

function exec_map!(f::Base.Callable, executor::DistributedExec, Y::AbstractVector, X::AbstractVector)
    Y .= pmap(f, executor.workers, X, distributed = true, batch_size = executor.batchsize)
    return Y
end