text:
module _FuzzEn2D
export FuzzEn2D
using Statistics: mean, std
    """
        Fuzz2D = FuzzEn2D(Mat)

    Returns the bidimensional fuzzy entropy estimate (`Fuzz2D`) estimated for
    the data matrix (`Mat`) using the default parameters: time delay = 1,
    fuzzy function (Fx) = 'default', fuzzy function parameters (r) = [0.2, 2],
    logarithm = natural, template matrix size = [floor(H/10) floor(W/10)],
    (where H and W represent the height and width of the data matrix 'Mat') \n
    ** The minimum dimension size of Mat must be > 10.**

        Fuzz2D = FuzzEn2D(Mat::AbstractArray{T,2} where T<:Real; m::Union{Int,Tuple{Int,Int}}=floor.(Int, size(Mat)./10),
                            tau::Int=1, r::Union{Real,Tuple{Real,Real}}=(.2*std(Mat, corrected=false),2),
                                Fx::String="default", Logx::Real=exp(1), Lock::Bool=true)

    Returns the bidimensional fuzzy entropy (`Fuzz2D`) estimates for the data
    matrix (`Mat`) using the specified 'keyword' arguments:

    # Arguments:
    `m`     - Template submatrix dimensions, an integer scalar (i.e. the same
              height and width) or a two-element vector of integers
              [height, width] with a minimum value > 1.
              (default: [floor(H/10) floor(W/10)])  \n
    `tau`   - Time Delay, a positive integer   (default: 1)  \n
    `Fx`    - Fuzzy funtion name, one of the following:
              {`"sigmoid", "modsampen", "default", "gudermannian", "linear"`}  \n
    `r`     - Fuzzy function parameters, a 1 element scalar or a 2 element
              vector of positive values. The 'r' parameters for each fuzzy
              function are defined as follows:\n
              sigmoid:      r(1) = divisor of the exponential argument
                            r(2) = value subtracted from argument (pre-division)
              modsampen:    r(1) = divisor of the exponential argument
                            r(2) = value subtracted from argument (pre-division)
              default:      r(1) = divisor of the exponential argument
                            r(2) = argument exponent (pre-division)
              gudermannian:   r  = a scalar whose value is the numerator of
                                   argument to gudermannian function:
                                   GD(x) = atan(tanh(r/x))
              linear:         r  = an integer value. When r = 0, the
                                   argument of the exponential function is
                                   normalised between [0 1]. When r = 1,
                                   the minimuum value of the exponential
                                   argument is set to 0. \n
    `Logx`  - Logarithm base, a positive scalar    (default: natural)\n
    `Lock`  - By default, FuzzEn2D only permits matrices with a maximum
              size of 128 x 128 to prevent memory errors when storing data on RAM.
              e.g. For Mat = [200 x 200], m = 3, and tau = 1, FuzzEn2D
              creates a vector of 753049836 elements. To enable matrices
              greater than [128 x 128] elements, set `Lock` to false.
              (default: true)\n
              ` WARNING: unlocking the permitted matrix size may cause
              your Julia IDE to crash.` \n

    # See also `SampEn2D`, `FuzzEn`, `XFuzzEn`

    # References:
        [1] Luiz Fernando Segato Dos Santos, et al.,
            "Multidimensional and fuzzy sample entropy (SampEnMF) for
            quantifying H&E histological images of colorectal cancer."
            Computers in biology and medicine
            103 (2018): 148-160.

        [2] Mirvana Hilal and Anne Humeau-Heurtier,
            "Bidimensional fuzzy entropy: Principle analysis and biomedical
            applications."
            41st Annual International Conference of the IEEE (EMBC) Society
            2019.

    """
    function FuzzEn2D(Mat::AbstractArray{T,2} where T<:Real;
        m::Union{Int,Tuple{Int,Int}}=floor.(Int, size(Mat)./10),
        tau::Int=1, r::Union{Real,Tuple{Real,Real}}=(.2*std(Mat,corrected=false),2),
        Fx::String="default", Logx::Real=exp(1), Lock::Bool=true)

    NL, NW = size(Mat)
    ((NL > 128 || NW > 128) && Lock) ?
        error("To prevent memory errors, matrix width & length must have <= 128 elements.
            To estimate FuzzEn2D  for the current matrix ($NL,$NW) change Lock to 'false'.
            Caution: unlocking the safe matrix size may cause the Julia IDE to crash.") :
            nothing

    length(m)==1 ?  (mL = m; mW = m) : (mL = m[1];  mW = m[2])
    (NL > 10 && NW > 10) ? nothing :
    error("Number of rows and columns in Mat must be > 10")
    (minimum(m)>1) ? nothing : error("m:  must be an integer > 1, or a 2 element integer tuple > 1")
    (tau > 0) ? nothing : error("tau:   must be an integer > 0")
    (minimum(r)>=0) ? nothing : error("r:     must be a positive value")
        #error("r:     must be 2 element tuple of positive values")
    (lowercase(Fx) in ["default","sigmoid","modsampen","gudermannian","linear"]) ?
        nothing : error("Fx:    must be one of the following strings -
        'default', 'sigmoid', 'modsampen', 'gudermannian', 'linear'")
    (Logx>0) ? nothing :  error("Logx:     must be a positive number > 0")

    if length(r) == 2 && lowercase(Fx)=="linear"
        r = 0;
        print("Multiple values for r entered. Default value (0) used.\n")
    elseif length(r) == 2 && lowercase(Fx)=="gudermannian"
        r = r[1]
        print("Multiple values for r entered. First value used.\n")
    end

    Fun = getfield(_FuzzEn2D,Symbol(lowercase(Fx)))
    NL = NL - mL*tau
    NW = NW - mW*tau
    X1 = zeros(mL,mW,NL*NW)
    X2 = zeros(mL+1,mW+1,NL*NW)
    p = 0;
    for k = 1:NL
        for n = 1:NW
            p += 1
            Temp2 = Mat[k:tau:(mL)*tau+k,n:tau:(mW)*tau+n]
            Temp1 = Temp2[1:end-1,1:end-1]
            X1[:,:,p] = Temp1 .- mean(Temp1)
            X2[:,:,p] = Temp2 .- mean(Temp2)
        end
    end
    p = size(X1,3)
    p != NL*NW ? @warn("Potential error with submatrix division.") : nothing
    Ny = p*(p-1)/2;
    Ny > 300000000 ? @warn("Number of pairwise distance calculations is $Ny") : nothing
    Y1 = zeros(p-1)
    Y2 = zeros(p-1)
    for k = 1:p-1
        Temp1 = maximum(abs.(X1[:,:,k+1:end] .- X1[:,:,k]),dims=(1,2))
        Y1[k] = sum(Fun(Temp1[:], r))
        Temp2 = maximum(abs.(X2[:,:,k+1:end] .- X2[:,:,k]),dims=(1,2))
        Y2[k] = sum(Fun(Temp2[:], r))
    end

    Fuzz2D = -log(Logx, sum(Y2)/sum(Y1))
    return Fuzz2D
    end

    function sigmoid(x,r)
        if length(r) == 1
            error("When Fx = 'Sigmoid', r must be a two-element vector.")
        end
        y = inv.(1 .+ exp.((x.-r[2])/r[1]))
        return y
    end

    function modsampen(x,r)
        if length(r) == 1
            error("When Fx = 'Modsampen', r must be a two-element vector.")
        end
        y = inv.(1 .+ exp.((x.-r[2])/r[1]))
        return y
    end

    function default(x,r)
        if length(r) == 1
            error("When Fx = 'Default', r must be a two-element vector.")
        end
        y = exp.(-(x.^r[2])/r[1])
        return y
    end

    function gudermannian(x,r)
        if r <= 0
            error("When Fx = 'Gudermannian', r must be a scalar > 0.")
        end
        y = atan.(tanh.(r[1]./x))
        y ./= maximum(y)
        return y
    end

    function linear(x,r)
        if r == 0 && length(x)>1
            y = exp.(-(x .- minimum(x))/(maximum(x)-minimum(x)))
        elseif r == 1
            y = exp.(-(x .- minimum(x)))
        elseif r == 0 && length(x)==1
            y = [0]
        else
            error("When Fx = 'Linear', r must be 0 or 1.")
        end
        return y
    end

end

"""
Copyright 2021 Matthew W. Flood, EntropyHub

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

    http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.

For Terms of Use see https://github.com/MattWillFlood/EntropyHub
"""


text:
using AMQPClient
using JSON
using Unmarshal
using DataFrames

const DLQ_EXCHANGE = "DLQ"

struct QueueDeclarationError <: Exception
    queue_name::String
end

function declare_queue!(chan, queue)
    success, queue_name, message_count, consumer_count = queue_declare(chan, queue)
    success || throw(QueueDeclarationError(queue_name))
    queue_name
end

serialize(serializable::Vector{UInt8}) = serializable
serialize(serializable::String) = Vector{UInt8}(serializable)
serialize(serializable::Int) = serialize(string(serializable))

compose(message::T where T) = Message(serialize(message), content_type="text/plain", delivery_mode=PERSISTENT)

send!(chan,exchange::String, message::T where T, routing_key::String) =
    basic_publish(chan, compose(message); exchange=exchange, routing_key=routing_key)

send!(chan, exchange::String, message::T where T) =
    send!(chan,exchange, message, "*")

AckUnit(::Type{Message}, src) =
    msg::Message -> AckUnit{Message}(msg, () -> src.dlFn(msg.data,msg.delivery_tag), () -> src.ackFn(msg.delivery_tag))

struct QueueSource <: Source{Message}
    name::String
    dlq::String
    pollFn::Function
    dlFn::Function
    ackFn::Function
    QueueSource(chan, name) = begin
        consumer_amqp_config!(chan, name)
        dlq_amqp_config!(chan, name)
        dlq_sink = dlq_sink!(name)
        pollFn = () -> basic_get(chan, name, false)
        ackFn = ack_handle -> basic_ack(chan, ack_handle)
        dlFn = (msg_data, ack_handle) -> (dlq_sink.sendFn(msg_data); ackFn(ack_handle))
        new(name, dlq, pollFn, dlFn, ackFn)
    end
end

source!(chan, name::String) = ack_source!(QueueSource(chan, name))
source!(chan, name::String, type::Type{String}) = source!(chan, name) |> map_ack(String, msg -> String(msg.data))
source!(chan, name::String, type::Type{Dict{String,Any}}) = source!(chan, name, String) |> safe_map_ack!(Dict{String,Any}, msg -> JSON.parse(msg))
source!(chan, name::String, type::Type{DataFrame}) = source!(chan, name, Dict{String,Any}) |> safe_map_ack!(DataFrame, msg -> DataFrame(msg))
source!(chan, name::String, unmar::Type{T}) where T = source!(chan, name, Dict{String,Any}) |> safe_map_ack!(type, msg -> Unmarshal.unmarshal(type, msg))

_queues = []

queue!(queue::Function) = push!(_queues, queue)

function execute_queues!()
    tasks = Task[]
    @async begin
        amqp_connection!(_conn_def) do conn
            @sync for queue in _queues
                 @async push!(tasks,
                    Threads.@spawn amqp_channel!(conn) do chan
                        queue(chan)
                    end)
            end

        end
    end
    return tasks
end

ext: jl


text:
module ParticlePlot

export kde, heatBoxPlot, densityplot, pplot, pploti
using RecipesBase, StatsBase
"""
Weighted kernel density estimate of the data `x` ∈ ℜN with weights `w` ∈ ℜN
`xi, densityw, density = kde(x,w)`
The number of grid points is chosen automatically and will approximately be equal to N/3
The bandwidth of the Gaussian kernel is chosen based on Silverman's rule of thumb
returns both weighted and non-weighted densities `densityw, density`
"""
function kde(x,w)
    @assert all(w .>= 0) "All weights must be non-negative"
    x   = x[:]
    e   = StatsBase.histrange(x,ceil(Int,length(x)/3))
    nb  = length(e)-1
    np  = length(x)
    I   = sortperm(x)
    s   = (rand()/nb+0):1/nb:e[end]
    j   = zeros(Float64,nb)
    bo = 1
    ilast = 0
    for i = 1:np
        ii = I[i]
        for b = bo:nb
            if x[ii] <= e[b+1]
                j[b] += w[ii]
                bo = b
                ilast = i
                break
            end
        end
    end
    j[end] += sum(w[I[ilast+1:np]])
    @assert all(j .>= 0) "j"
    xi      = e[1:end-1]+0.5(e[2]-e[1])
    σ       = std(x)
    if σ == 0 # All mass on one point
        return x,ones(x),ones(x)
    end
    h        = 1.06σ*np^(-1/5)
    K(x)     = exp(-(x/h).^2/2) / √(2π)
    densityw = [1/h*sum(j.*K(xi.-xi[i])) for i=1:nb] # This density is effectively normalized by nb due to sum(j) = 1
    density  = [1/(np*h)*sum(K(x.-xi[i])) for i=1:nb]

    # @assert all(densityw .>= 0) "densityw"
    # @assert all(density  .>= 0) "density"

    return xi, densityw, density

end

@userplot DensityPlot
@recipe function f(dp::DensityPlot)
    x = dp.args[1]
    w = length(dp.args) > 1 ? dp.args[2] : ones(x)/length(x)
    xi, densityw, density = kde(x,w)
    if maximum(x)-minimum(x) > 0
        title --> "Kernel density estimate"
        seriestype := :path
        @series begin
            label --> "Weighted density"
            linecolor --> :blue
            xi,densityw
        end
        @series begin
            label --> "Non-weighted density"
            linecolor --> :red
            xi,density
        end
    end
end

function heatBoxPlot(plt, x, t, minmax, nbinsy=30)
    if maximum(x)-minimum(x) > 0
        heatmap!(plt, x=(t-1):1/length(x):t,y=x',nbins=(1,nbinsy), ylims=tuple(minmax...))
    end
    return
end

""" `pplot(x, w, y, yhat, N, a, t, xreal, xhat, xOld, pdata)`
To be called inside a particle filter, plots either particle density (`density=true`) or individual particles (`density=false`) \n
Will plot all the real states in `xIndices` as well as the expected vs real measurements of `yIndices`.
Arguments: \n
* `x`: `Array(M,N)`. The states for each patyicle where `M` number of states, `N` number of Particles
* `w`: `Array(N)`. weight of each particle
* `y`: `Array(R,T)`. All true outputs. `R` is number of outputs, `T` is total number of time steps (will only use index `t`)
* `yhat`: `Array(R,N)` The expected output per particle. `R` is number of outputs, `N` number of Particles
* `N`, Number of particles
* `a`, `Array(N)`, reorderng of particles (e.g. `1:N`)
* `t`, Current time step
* `xreal`: `Array(M,T)`. All true states. `R` is number of states, `T` is total number of time steps (will only use index `t`)
* `xhat`: Not used
* `xOld`: Same as `x`, but for previous time step, only used when `!density` to show states origins
* `pdata`: Persistant data for plotting. Set to void in first call and pdataOut on remaining \n
* `density = true` To only plot the particle trajectories, set (`leftOnly=false`)\n
* `leftOnly = false`\n
* `xIndices = 1:size(x,1)`\n
* `yIndices = 1:size(y,1)`\n
Returns: `pdataOut`
"""
function pplot(x, w, y, yhat, N, a, t, xreal, xhat, xOld, pdata; density = true, leftOnly = false, xIndices = 1:size(x,1), yIndices = 1:size(y,1), slidef=0.9)
    immerse()
    cols = leftOnly?1:2
    grd = (r,c) -> (r-1)*cols+c
    println("Surviving: "*string((N-length(setdiff(Set(1:N),Set(a))))/N))
    vals = [x;yhat]
    realVals = [xreal;y]
    if !density
        valsOld = [xOld;yhat]
    end

    pltIdx = [xIndices; size(x,1)+yIndices]
    if pdata == Void
        pdata = (subplot(layout=cols*ones(Int,length(pltIdx))), zeros(length(pltIdx),2))
        gui(pdata[1])
    end
    p, minmax = pdata
    dataMin = minimum(vals[pltIdx,:],2)
    dataMax = maximum(vals[pltIdx,:],2)
    minmax = [min(minmax[:,1], dataMin)*slidef+(1-slidef)*dataMin max(minmax[:,2], dataMax)*slidef+(1-slidef)*dataMax]
    #c = (w[:]-minimum(w))*3
    c = w[:]*5*N

    for (i, val) in enumerate(pltIdx)
        if !leftOnly
            oldFig = p.plts[grd(i,2)].o[1]
            newPlot = plot()
        end
        #Plot the heatmap on the left plot
        density && heatBoxPlot(p.plts[grd(i,1)], vals[val,:], t, minmax[i,:])

        for j = 1:N
            if !density
                #Plot the line on the left plot
                plot!(p.plts[grd(i,1)], [t-1.5,t-1], [valsOld[val,j], valsOld[val,j]], legend=false)
                plot!(p.plts[grd(i,1)], [t-1,t-0.5], [valsOld[val,a[j]], vals[val,j]], legend=false)
            end
            #Plot each of the dots on the right side
            !leftOnly && !density && plot!(newPlot, [j,j], [vals[val,j]',vals[val,j]'], marker=(:circle, :red, c[j])  , legend=false)
        end
        if !leftOnly
            density && densityplot!(newPlot, vals[val,:], w , ylims = tuple(minmax[i,:]...), c=:blue)
            #Fix the bug with updating plot
            p.plts[grd(i,2)] = copy(newPlot)
            p.plts[grd(i,2)].o = (oldFig, p.plts[grd(i,2)].o[2])
        end
        #Plot Real State Here
        plot!(p.plts[grd(i,1)], (t-1):t, [realVals[val,t], realVals[val,t]], legend=false, color=:black, linewidth=5)
    end
    gui(p)
    (p, minmax)
end

# """ `pploti(x, w, y, yhat, N, a, t, xreal, xhat, xOld, pdata)`
# Same function as pplot but with options to skip, wait and quit.
# """
function pploti(x, w, y, yhat, N, a, t, xreal, xhat, xOld, pdata; density = true, leftOnly = false, xIndices = 1:size(x,1), yIndices = 1:size(y,1), slidef=0.9)
    cols = leftOnly?1:2
    pltIdx = [xIndices; size(x,1)+yIndices]
    if pdata == Void
        pdata = (subplot(layout=cols*ones(Int,length(pltIdx))), Array{Float64,2}(length(pltIdx),2),0,0)
        gui(pdata[1])
    end
    p, minmax, skip, skipped = pdata
    if skip > 0 && skipped < skip
        skipped += 1
        return (p, minmax, skip, skipped)
    elseif skip > 0
        skipped = 0
        skip = 0
    end
    p, minmax = pplot(x, w, y, yhat, N, a, t, xreal, xhat, xOld, (p, minmax), density = density, leftOnly = leftOnly, xIndices = xIndices, yIndices = yIndices)
    print("Waiting for command. q to Quit, ^D to run all, s NN to skip NN steps:\n")
    line = readline(STDIN)
    if line == "q\n"
        error("Quitting")
    elseif contains(line, "s")
        ss = split(strip(line,'\n'))
        skip = parse(Int,ss[2])
    end
    (p, minmax, skip, skipped)
end
end # module


text:
using CSP
using Base.Test

c1 = Channel(Int)
c2 = Channel(Int)

t1 = @async begin
    sleep(0.5)
    for i=6:10
        @test recv(c1) == i
    end
    for i=1:5
        send(c1, i)
    end

    # time million writes
    println("sending a million integers")
    tic()
    for i = 1:1000000
        send(c2, i)
    end
    toc()
end

t2 = @async begin
    for i=6:10
        send(c1, i)
    end
    for i=1:5
        @test recv(c1) == i
    end

    for i = 1:1000000
        recv(c2)
    end
end

wait(t1)
wait(t2)


text:
module ArgDecomposition

export flatten_to_tuple, reassemble_tuple

struct StaticType{T} end
gettype(::StaticType{T}) where {T} = T

function _append_fields!(t::Expr, body::Expr, sym::Symbol, ::Type{T}) where {T}
  for f ∈ 1:fieldcount(T)
    TF = fieldtype(T, f)
    Base.issingletontype(TF) && continue
    gfcall = Expr(:call, getfield, sym, f)
    if fieldcount(TF) ≡ 0
      push!(t.args, gfcall)
    elseif TF <: DataType
      push!(t.args, :(StaticType{$gfcall}()))
    else
      newsym = gensym(sym)
      push!(body.args, Expr(:(=), newsym, gfcall))
      _append_fields!(t, body, newsym, TF)
    end
  end
  return nothing
end
@generated function flatten_to_tuple(r::T) where {T}
  body = Expr(:block, Expr(:meta,:inline))
  t = Expr(:tuple)
  if Base.issingletontype(T)
    nothing
  elseif fieldcount(T) ≡ 0
    push!(t.args, :r)
  elseif T <: DataType
    push!(t.args, :(StaticType{r}()))
  else
    _append_fields!(t, body, :r, T)
  end
  push!(body.args, t)
  body
end
function rebuild_fields(offset::Int, ::Type{T}) where {T}
  call = (T <: Tuple) ? Expr(:tuple) : Expr(:new, T)
  for f ∈ 1:fieldcount(T)
    TF = fieldtype(T, f)
    if Base.issingletontype(TF)
      push!(call.args, TF.instance)
    elseif fieldcount(TF) ≡ 0
      push!(call.args, :($getfield(t, $(offset += 1), false)))
    elseif TF <: DataType
      push!(call.args, :($gettype($getfield(t, $(offset += 1), false))))
    else
      arg, offset = rebuild_fields(offset, TF)
      push!(call.args, arg)
    end
  end
  return call, offset
end
@generated function reassemble_tuple(::Type{T}, t::Tuple) where {T}
  if Base.issingletontype(T)
    return T.instance
  elseif fieldcount(T) ≡ 0
    call = :($getfield(t, 1, false))
  elseif T <: DataType
    call = :($gettype($getfield(t, 1, false)))
  else
    call, _ = rebuild_fields(0, T)
  end
  Expr(:block, Expr(:meta,:inline), call)
end

end


text:
using Test
using PlayingCards
using PlayingCards: rank_string

@testset "Ranks" begin
    for v in ranks()
        v==1 && continue
        @test high_value(v*♣) == low_value(v*♣) == v
    end
    @test low_value(1*♣) == 1
end

@testset "Color" begin
    @test color(J♣) == :black
    @test color(A♠) == :black
    @test color(♣) == :black
    @test color(♠) == :black
    @test color(J♢) == :red
    @test color(A♡) == :red
    @test color(♢) == :red
    @test color(♡) == :red
end

@testset "Card" begin
    @test rank(J♣) == 11
    @test suit(J♣) == ♣
    @test suit(J♡) == ♡
    @test suit(J♢) == ♢
    @test suit(J♠) == ♠
    @test suit(A♣) == ♣
    @test suit(A♡) == ♡
    @test suit(A♢) == ♢
    @test suit(A♠) == ♠
    @test_throws ArgumentError 14*♣
    @test_throws ArgumentError 0*♣
end

@testset "strings" begin
    @test sprint(show, 2♣) == "2♣"
    @test sprint(show, J♣) == "J♣"
    @test string(2♣) == "2♣"
    @test string(T♣) == "T♣"
    @test string(J♣) == "J♣"
    @test string(Q♣) == "Q♣"
    @test string(K♣) == "K♣"
    @test string(A♣) == "A♣"
    @test string(♣) == "♣"
    @test string(♠) == "♠"
    @test string(♡) == "♡"
    @test string(♢) == "♢"
    @test_throws AssertionError rank_string(Int8(-1))
end

@testset "Deck" begin
    deck = ordered_deck()
    @test length(deck) == 52
    @test iterate(deck) == iterate(deck.cards)
    shuffle!(deck)
    cards = pop!(deck, 2)
    @test length(cards)==2
    @test length(deck)==50
    @test length(full_deck())==52

    s="A♣ 2♣ 3♣ 4♣ 5♣ 6♣ 7♣ 8♣ 9♣ T♣ J♣ Q♣ K♣
A♠ 2♠ 3♠ 4♠ 5♠ 6♠ 7♠ 8♠ 9♠ T♠ J♠ Q♠ K♠
A♡ 2♡ 3♡ 4♡ 5♡ 6♡ 7♡ 8♡ 9♡ T♡ J♡ Q♡ K♡
A♢ 2♢ 3♢ 4♢ 5♢ 6♢ 7♢ 8♢ 9♢ T♢ J♢ Q♢ K♢
"
    @test sprint(show, ordered_deck()) == s

    deck = ordered_deck()
    five_spades = pop!(deck, 5♠)
    @test five_spades === 5♠
    @test length(deck) == 51
    @test findfirst(x->x==5♠, deck.cards) == nothing
    @test_throws ErrorException pop!(deck, 5♠)
end

@testset "Allocations" begin
    alloc = @allocated ordered_deck()
    @test alloc == 304
end


text:
@testset "VGG ($model)" for model in [vgg11, vgg11bn, vgg13, vgg13bn, vgg16, vgg16bn, vgg19, vgg19bn]
  imsize = (224, 224)
  m = model(imsize)

  @test size(m(rand(Float32, imsize..., 3, 50))) == (1000, 50)
end


text:
function _get_expected_hashes(julia_version_to_expected_hashes::AbstractDict, julia_version_number::VersionNumber)
    up_to_patch = VersionNumber(julia_version_number.major, julia_version_number.minor, julia_version_number.patch)
    up_to_minor = VersionNumber(julia_version_number.major, julia_version_number.minor)
    up_to_major = VersionNumber(julia_version_number.major)
    if haskey(julia_version_to_expected_hashes, julia_version_number)
        return julia_version_to_expected_hashes[julia_version_number]
    elseif haskey(julia_version_to_expected_hashes, up_to_patch)
        return julia_version_to_expected_hashes[up_to_patch]
    elseif haskey(julia_version_to_expected_hashes, up_to_minor)
        return julia_version_to_expected_hashes[up_to_minor]
    else
        return julia_version_to_expected_hashes[up_to_major]
    end
end


text:
using Weave
weave("weave-source/index.jmd", out_path="docs", doctype="github")

max_stars_repo_name: mpastell/weave-ghpages-sample

max_issues_repo_name: mpastell/weave-ghpages-sample

max_forks_repo_name: mpastell/weave-ghpages-sample


text:
function CombinedEigenfunction(
    domain::AbstractPlanarDomain,
    us::Vector{<:AbstractPlanarEigenfunction},
    orders::Vector{Int} = ones(Int, length(us));
    us_to_boundary = [active_boundaries(domain, u) for u in us],
    even_boundaries = BitSet(),
)
    boundary_to_us = OrderedDict(
        i => BitSet(findall(B -> i ∈ B, us_to_boundary)) for i in boundaries(domain)
    )

    return CombinedEigenfunction(
        domain,
        us,
        boundary_to_us,
        BitSet(even_boundaries),
        orders,
    )
end

function CombinedEigenfunction(
    domain::IntersectedDomain,
    orders::Vector{Int} = ones(Int, length(boundaries(domain))),
)
    boundary_to_us = OrderedDict(i => BitSet() for i in boundaries(domain))
    us = AbstractPlanarEigenfunction[]

    for i in exterior_boundaries(domain)
        # TODO: Handle other types
        @assert typeof(domain.exterior) <: Triangle ||
                typeof(domain.exterior) <: TransformedDomain{<:Triangle}

        v = StandaloneVertexEigenfunction(domain.exterior, i)
        push!(us, v)
        idx = length(us)
        push!(boundary_to_us[i], idx)
        for j in interior_boundaries(domain)
            push!(boundary_to_us[j], idx)
        end
    end

    for i in interior_boundaries(domain)
        _, ip = get_domain_and_boundary(domain, i)
        # TODO: Handle other types
        @assert typeof(domain.exterior) <: Triangle ||
                typeof(domain.exterior) <: TransformedDomain{<:Triangle}

        v = StandaloneVertexEigenfunction(domain.interior, ip, outside = true)
        push!(us, v)

        idx = length(us)
        push!(boundary_to_us[i], length(us))
        for j in exterior_boundaries(domain)
            push!(boundary_to_us[j], idx)
        end
    end

    # TODO: Add interior eigenfunction
    v = StandaloneInteriorEigenfunction(domain)
    push!(us, v)
    idx = length(us)
    for j in boundaries(domain)
        push!(boundary_to_us[j], idx)
    end

    return CombinedEigenfunction(domain, us, boundary_to_us, BitSet(), orders)
end

function Base.show(io::IO, ::MIME"text/plain", u::CombinedEigenfunction)
    print(
        io,
        "Combined eigenfunction on $(typeof(u.domain)) with $(length(coefficients(u))) coefficients",
    )
    if !get(io, :compact, false)
        println(io, "")
        recur_io = IOContext(io, :compact => true)
        for v in u.us
            print(io, "- ")
            show(recur_io, v)
        end
    end
end

Base.show(io::IO, u::CombinedEigenfunction) = print(
    io,
    "Combined eigenfunction on $(typeof(u.domain)) with $(length(coefficients(u))) coefficients",
)

function set_domain!(u::CombinedEigenfunction, domain::AbstractDomain)
    u.domain = domain
    for v in u.us
        set_domain!(v, domain)
    end
    return u
end

function recompute!(u::CombinedEigenfunction)
    for v in u.us
        recompute!(v)
    end
    return u
end

function active_boundaries(domain::AbstractPlanarDomain, u::CombinedEigenfunction)
    if domain === u.domain
        return findall(!isempty, u.boundary_to_us)
    else
        return boundaries(domain)
    end
end

"""
    basis_function(u::CombinedEigenfunction, k::Integer)

Return the index for the eigenfunctions that `u` is a combination of
that should be used and the index to use for that eigenfunction.

This depends on the relative order between the eigenfunctions. It
returns `(i, j)` which corresponds to calling the `j`-th
basis-function of `u.us[i]`.

TODO: This could be optimized a lot more.
"""
function basis_function(u::CombinedEigenfunction, k::Integer)
    k > 0 || throw(ArgumentError("k must be positive not $k"))
    i = 1
    j = 0
    while k > u.orders[i]
        k -= u.orders[i]

        i += 1
        if i > length(u.orders)
            i = 1
            j += 1
        end
    end

    return i, u.orders[i] * j + k
end

"""
    basis_function(u::CombinedEigenfunction, ks::UnitRange{Int})

Return a vector with elements of type `UnitRange{Int}` where element
`i` is the range of indices that should be used for the respective
eigenfunctions, i.e. the eigenfunction `u.us[i]` should take the range
given by the `i`th element in this vector.
"""
function basis_function(u::CombinedEigenfunction, ks::UnitRange{Int})
    ks.start == 1 || throw(ArgumentError("ks must start with 1, got ks = $ks"))
    ks.stop == 0 && return fill(1:0, length(u.us))

    A = [0; cumsum(u.orders)]
    fullcycles = div(ks.stop, A[end])
    R = ks.stop % A[end]
    remaining = [max(min(u.orders[i], R - A[i]), 0) for i in eachindex(u.us)]
    return [1:(fullcycles*u.orders[i]+remaining[i]) for i in eachindex(u.us)]
end

function coefficients(u::CombinedEigenfunction)
    coeffs = [coefficients(v) for v in u.us]
    N = sum(length, coeffs)
    return [coeffs[i][j] for (i, j) in [basis_function(u, k) for k = 1:N]]
end

function set_eigenfunction!(u::CombinedEigenfunction, coefficients::Vector)
    is = [basis_function(u, k)[1] for k = 1:length(coefficients)]
    for i = 1:length(u.us)
        set_eigenfunction!(u.us[i], coefficients[is.==i])
    end
    return u
end

function (u::CombinedEigenfunction)(
    xy::AbstractVector{T},
    λ::arb,
    k::Integer;
    boundary = nothing,
    notransform::Bool = false,
) where {T<:Union{arb,arb_series}}
    i, j = basis_function(u, k)
    if !isnothing(boundary) && !(i ∈ u.boundary_to_us[boundary])
        if T == arb
            return u.domain.parent(0)
        else
            return 0 * xy[1]
        end
    end

    return u.us[i](xy, λ, j, boundary = boundary, notransform = notransform)
end

function (u::CombinedEigenfunction)(
    xy::AbstractVector{T},
    λ::arb,
    ks::UnitRange{Int};
    boundary = nothing,
    notransform::Bool = false,
) where {T<:Union{arb,arb_series}}
    if isnothing(boundary)
        indices = eachindex(u.us)
    else
        indices = u.boundary_to_us[boundary]
    end

    ks_per_index = basis_function(u, ks)

    res_per_index = similar(u.us, Vector{T})
    for i in eachindex(u.us)
        if i in indices
            res_per_index[i] =
                one(λ) .* u.us[i](
                    xy,
                    λ,
                    ks_per_index[i],
                    boundary = boundary,
                    notransform = notransform,
                )
        else
            z = T == arb ? zero(λ) : 0 * xy[1]
            res_per_index[i] = fill(z, length(ks_per_index[i]))
        end
    end

    # TODO: This can be done more efficiently
    res = similar(ks, T)
    for i in eachindex(ks)
        j, l = basis_function(u, ks[i])
        res[i] = res_per_index[j][l]
    end

    return res
end

function (u::CombinedEigenfunction)(
    xy::AbstractVector{T},
    λ::arb;
    boundary = nothing,
    notransform::Bool = false,
) where {T<:Union{arb,arb_series}}
    if isnothing(boundary)
        indices = eachindex(u.us)
    else
        indices = u.boundary_to_us[boundary]
    end

    res = u.domain.parent(0)
    for i in indices
        res += u.us[i](xy, λ, boundary = boundary, notransform = notransform)
        if (T == arb && !isfinite(res)) || (T == arb_series && !isfinite(res[end]))
            return res
        end
    end

    return res
end


text:
using Oceananigans
using OceanTurbulenceParameterEstimation
using Oceananigans.Units
using DataDeps
using GLMakie
using Printf

suite = "four_day_suite"
case_path(case) = @datadep_str("$(suite)_1m/$(case)_instantaneous_statistics.jld2")

cases = [
         "strong_wind_no_rotation",
         "free_convection",
         "weak_wind_strong_cooling",
         "strong_wind_weak_cooling",
         "strong_wind",
        ]

observations = [SyntheticObservations(case_path(case), field_names=(:u, :v, :T, :e)) for case in cases]

u_catke_time_series = FieldTimeSeries("catke_simulation_$suite.jld2", "u", boundary_conditions=nothing)
v_catke_time_series = FieldTimeSeries("catke_simulation_$suite.jld2", "v", boundary_conditions=nothing)
T_catke_time_series = FieldTimeSeries("catke_simulation_$suite.jld2", "T", boundary_conditions=nothing)
e_catke_time_series = FieldTimeSeries("catke_simulation_$suite.jld2", "e", boundary_conditions=nothing)

#####
##### Animate!
#####

fig = Figure(resolution=(1800, 800))

T_axs = []
u_axs = []
e_axs = []

for (i, case) in enumerate(cases)
    ax_l = Label(fig[1, i+1], replace(case, "_" => "\n"), tellwidth=false)
    ax_T = Axis(fig[2, i+1])
    ax_u = Axis(fig[3, i+1])
    ax_e = Axis(fig[4, i+1])

    xlims!(ax_e, 0, 0.005)
    xlims!(ax_u, -0.2, 0.3)

    hidedecorations!(ax_T)
    hidedecorations!(ax_u)
    hidedecorations!(ax_e)

    hidespines!(ax_T)
    hidespines!(ax_u)
    hidespines!(ax_e)

    push!(T_axs, ax_T)
    push!(u_axs, ax_u)
    push!(e_axs, ax_e)
end

Nt = length(first(observations).times)
slider = Slider(fig[5, :], range=1:Nt, startvalue=1)
n = slider.value

###
### Vertical velocity
###

z = znodes(T_catke_time_series)

for (j, case) in enumerate(cases)
    obs = observations[j]

    u_les = @lift interior(obs.field_time_serieses.u[$n])[1, 1, :]
    v_les = @lift interior(obs.field_time_serieses.v[$n])[1, 1, :]
    T_les = @lift interior(obs.field_time_serieses.T[$n])[1, 1, :]
    e_les = @lift interior(obs.field_time_serieses.e[$n])[1, 1, :]

    ax_T = T_axs[j]
    ax_u = u_axs[j]
    ax_e = e_axs[j]

    linewidth = 5
    color = :indigo
    lines!(ax_T, T_les, z; linewidth, color=(color, 0.8))
    lines!(ax_u, u_les, z; linewidth, color=(color, 0.8))
    lines!(ax_u, v_les, z; linewidth, color=(color, 0.5))
    lines!(ax_e, e_les, z; linewidth, color=(color, 0.8))

    linewidth = 2
    colors = [:darkorange3, :deeppink, :lightseagreen]
    for i = 1:3
        u_catke = @lift interior(u_catke_time_series[$n])[i, j, :]
        v_catke = @lift interior(v_catke_time_series[$n])[i, j, :]
        T_catke = @lift interior(T_catke_time_series[$n])[i, j, :]
        e_catke = @lift interior(e_catke_time_series[$n])[i, j, :]

        color = colors[i]
        lines!(ax_T, T_catke, z; linewidth, color)
        lines!(ax_u, u_catke, z; linewidth, color)
        lines!(ax_u, v_catke, z; linewidth, color=(color, 0.5))
        lines!(ax_e, e_catke, z; linewidth, color)
    end
end

Label(fig[2, 1], "Temperature (ᵒC)", tellheight=false)
Label(fig[3, 1], "Velocity (m s⁻¹)", tellheight=false)
Label(fig[4, 1], "Turbulent \n kinetic \n energy", tellheight=false)

display(fig)


text:
module Buoyancy

export
    BuoyancyTracer, SeawaterBuoyancy, buoyancy_perturbation,
    LinearEquationOfState, RoquetIdealizedNonlinearEquationOfState, TEOS10,
    ∂x_b, ∂y_b, ∂z_b, buoyancy_perturbation, buoyancy_frequency_squared,
    BuoyancyField

using Printf
using Oceananigans.Grids
using Oceananigans.Operators

import SeawaterPolynomials: ρ′, thermal_expansion, haline_contraction

# Physical constants
# https://en.wikipedia.org/wiki/Gravitational_acceleration#Gravity_model_for_Earth (30 Oct 2019)
const g_Earth = 9.80665

"""
    AbstractBuoyancy{EOS}

Abstract supertype for buoyancy models.
"""
abstract type AbstractBuoyancy{EOS} end

"""
    AbstractEquationOfState

Abstract supertype for equations of state.
"""
abstract type AbstractEquationOfState end

function validate_buoyancy(buoyancy, tracers)
    req_tracers = required_tracers(buoyancy)

    all(tracer ∈ tracers for tracer in req_tracers) ||
        error("$(req_tracers) must be among the list of tracers to use $(typeof(buoyancy).name.wrapper)")

    return nothing
end

validate_buoyancy(::Nothing, tracers) = nothing

required_tracers(::Nothing) = ()

@inline buoyancy_perturbation(i, j, k, grid::AbstractGrid{FT}, ::Nothing, C) where FT = zero(FT)

@inline ∂x_b(i, j, k, grid::AbstractGrid{FT}, ::Nothing, C) where FT = zero(FT)
@inline ∂y_b(i, j, k, grid::AbstractGrid{FT}, ::Nothing, C) where FT = zero(FT)
@inline ∂z_b(i, j, k, grid::AbstractGrid{FT}, ::Nothing, C) where FT = zero(FT)

"""
    BuoyancyTracer <: AbstractBuoyancy{Nothing}

Type indicating that the tracer `b` represents buoyancy.
"""
struct BuoyancyTracer <: AbstractBuoyancy{Nothing} end

required_tracers(::BuoyancyTracer) = (:b,)

@inline buoyancy_perturbation(i, j, k, grid, ::BuoyancyTracer, C) = @inbounds C.b[i, j, k]

@inline ∂x_b(i, j, k, grid, ::BuoyancyTracer, C) = ∂xᶠᵃᵃ(i, j, k, grid, C.b)
@inline ∂y_b(i, j, k, grid, ::BuoyancyTracer, C) = ∂yᵃᶠᵃ(i, j, k, grid, C.b)
@inline ∂z_b(i, j, k, grid, ::BuoyancyTracer, C) = ∂zᵃᵃᶠ(i, j, k, grid, C.b)

include("seawater_buoyancy.jl")
include("linear_equation_of_state.jl")
include("nonlinear_equation_of_state.jl")

Base.show(io::IO, b::SeawaterBuoyancy{FT}) where FT =
    println(io, "SeawaterBuoyancy{$FT}: g = $(b.gravitational_acceleration)", '\n',
                "└── equation of state: $(b.equation_of_state)")

Base.show(io::IO, eos::LinearEquationOfState{FT}) where FT =
    println(io, "LinearEquationOfState{$FT}: ", @sprintf("α = %.2e, β = %.2e", eos.α, eos.β))

include("buoyancy_field.jl")

end


text:
module TestDifferentiate

using ModiaBase
using Test

removeBlock(ex) = ex
function removeBlock(ex::Expr)
    if ex.head in [:block]
        ex.args[1]
    else
        Expr(ex.head, [removeBlock(arg) for arg in ex.args]...)
    end
end

function showDifferentiate(ex, timeInvariants=[])
    Base.remove_linenums!(ex)
    print(removeBlock(ex), "   :   ")
    der = ModiaBase.derivative(ex, timeInvariants)
    Base.remove_linenums!(der)
    der = string(removeBlock(der))
    println(der)
    string(der)
end

function testDifferentiate()
    println("\nTest differentiate")

    @test showDifferentiate(10) == "0"
    @test showDifferentiate(10.0) == "0"
    @test showDifferentiate(:time) == "1"
    @test showDifferentiate(:x) == "der(x)"
    @test showDifferentiate(:(der(x))) == "der(der(x))"
    println()

    @test showDifferentiate(:(x.y.z)) == "der(x.y.z)"
    @test showDifferentiate(:(x = y)) == "der(x) = der(y)"
    println()

    @test showDifferentiate(:(+ y)) == "+(der(y))"
    @test showDifferentiate(:(x + y)) == "der(x) + der(y)"
    @test showDifferentiate(:(x + y + z)) == "der(x) + der(y) + der(z)"
    @test showDifferentiate(:(- y)) == "-(der(y))"
    @test showDifferentiate(:(x - y)) == "der(x) - der(y)"
    @test showDifferentiate(:(x - y - z)) == "(der(x) - der(y)) - der(z)"

    @test showDifferentiate(:(x * y)) == "der(x) * y + x * der(y)"
    @test showDifferentiate(:(x * y * z)) == "der(x) * y * z + x * der(y) * z + x * y * der(z)"
    @test showDifferentiate(:(x / y)) == "(one(x) / y) * der(x) + -((x / y) / y) * der(y)"
#    @test showDifferentiate(:(x / y / z)) == "(one(x / y) / z) * ((one(x) / y) * der(x) + -((x / y) / y) * der(y)) + -(((x / y) / z) / z) * der(z)"
    @test showDifferentiate(:(x ^ 2)) == "(2 * x ^ (2 - 1)) * der(x)"
    @test showDifferentiate(:(x ^ y)) == "(y * x ^ (y - 1)) * der(x) + (x ^ y * log(x)) * der(y)"
#    @test showDifferentiate(:(x ^ y ^ z)) == "(y ^ z * x ^ (y ^ z - 1)) * der(x) + (x ^ (y ^ z) * log(x)) * ((z * y ^ (z - 1)) * der(y) + (y ^ z * log(y)) * der(z))"

    @test showDifferentiate(:(sin(x))) == "cos(x) * der(x)"

    @test showDifferentiate(:([x, y])) == "[der(x), der(y)]"

    @test showDifferentiate(:(if x; y else z end)) == "if x\n    der(y)\nelse\n    der(z)\nend"
    @test showDifferentiate(:(if x; y elseif z; v else w end)) == "if x\n    der(y)\nelseif z\n    der(v)\nelse\n    der(w)\nend"

    @test showDifferentiate(:(1+2+3)) == "0"
    @test showDifferentiate(:(1-2-3)) == "0"
    @test showDifferentiate(:(x-2-3)) == "der(x)"
    @test showDifferentiate(:(1-x-3)) == "-(der(x))"

    @test showDifferentiate(:(2*x + y + z)) == "2 * der(x) + der(y) + der(z)"

    # Previous test suit
    der = showDifferentiate(:(x + 5 + z = w))
    @test der == "der(x) + der(z) = der(w)"

#    der = showDifferentiate(ModiaBase.derivative(:(x + 5 + z = w)))
#    @test der == "der(der(x)) + der(der(z)) = der(der(w))"

    der = showDifferentiate(Expr(:(=), Expr(:call, :+, :x), :w))
    @test der == "+(der(x)) = der(w)"

    der = showDifferentiate(:(2 + 3 = w))
    @test der == "0 = der(w)"

    der = showDifferentiate(Expr(:(=), Expr(:call, :-, :x), :w))
    @test der == "-(der(x)) = der(w)"

    der = showDifferentiate(:(x - 5 - z = w))
    @test der == "der(x) - der(z) = der(w)"

    der = showDifferentiate(:(5 - x - z = w))
    @test der == "-(der(x)) - der(z) = der(w)"

    der = showDifferentiate(:(5x = w))
    @test der == "5 * der(x) = der(w)"

    der = showDifferentiate(:(x * 5 * z = w))
    @test der == "der(x) * 5 * z + x * 5 * der(z) = der(w)"

    der = showDifferentiate(:(4 * 5 * 6 = w))
    @test der == "0 = der(w)"

    der = showDifferentiate(:(y = x/5))
    @test der == "der(y) = (one(x) / 5) * der(x)"

    der = showDifferentiate(:(y = 5/y))
    @test der == "der(y) = -((5 / y) / y) * der(y)"

    der = showDifferentiate(:(y = [1, x]))
    @test der == "der(y) = [0, der(x)]"

    der = showDifferentiate(:(y = [2x 3x; 4x 5x]))
    @test der == "der(y) = [2 * der(x) 3 * der(x); 4 * der(x) 5 * der(x)]"

    der = showDifferentiate(:(y = [2*x 3x; 4x 5x]*[1, x]))
    @test der == "der(y) = [2 * der(x) 3 * der(x); 4 * der(x) 5 * der(x)] * [1, x] + [2x 3x; 4x 5x] * [0, der(x)]"
    #=
    der = showDifferentiate(:(y = transpose(B) + B´))
    @test der == "der(y) = transpose(der(B)) + der(B´)"
    =#
    der = showDifferentiate(:(y = x[5, 6]))
    @test der == "der(y) = (der(x))[5, 6]"

    der = showDifferentiate(:(y = x[5:7]))
    @test der == "der(y) = (der(x))[5:7]"

  #=    der = showDifferentiate(:(y = [x for x in z]))
    @test der == "der(y) = [x for x = der(z)]"

    der = showDifferentiate(:(y = [x[i] for i in 1:5]))
    @test der == "der(y) = [x[i] for i = nothing]"
  =#
    der = showDifferentiate(:(y = sin(x)))
    @test der == "der(y) = cos(x) * der(x)"

    der = showDifferentiate(:(y = cos(x)))
    @test der == "der(y) = -(sin(x)) * der(x)"

    #=
    der = showDifferentiate(:(y = tan(x)))
    @test der == "der(y) = (1 / cos(x) ^ 2) * der(x)"
    =#
    der = showDifferentiate(:(y = exp(x)))
    @test der == "der(y) = exp(x) * der(x)"

    der = showDifferentiate(:(z = x^y))
    @test der == "der(z) = (y * x ^ (y - 1)) * der(x) + (x ^ y * log(x)) * der(y)"

    der = showDifferentiate(:(y = log(x)))
    @test der == "der(y) = inv(x) * der(x)"

    der = showDifferentiate(:(y = asin(x)))
    @test der == "der(y) = inv(sqrt(1 - x ^ 2)) * der(x)"

    der = showDifferentiate(:(y = acos(x)))
    @test der == "der(y) = -(inv(sqrt(1 - x ^ 2))) * der(x)"

    der = showDifferentiate(:(y = atan(x)))
    @test der == "der(y) = inv(1 + x ^ 2) * der(x)"

    #=
    der = showDifferentiate(:(y = f(x, 5, z)))
    @test der == "der(y) = f_der_1(x, 5, z) * der(x) + f_der_3(x, 5, z) * der(z)"

    der = showDifferentiate(:(y = f(x, 5, g(z))))
    @test der == "der(y) = f_der_1(x, 5, g(z)) * der(x) + f_der_3(x, 5, g(z)) * (g_der_1(z) * der(z))"
    =#
    der = showDifferentiate(:(y = true ? x : y))
      @test der == "der(y) = if true\n        der(x)\n    else\n        der(y)\n    end"

    #=
    der = showDifferentiate(:(y = if b; x elseif false; y else z end))
    @test der == "der(y) = if b\n        der(x)\n    else\n        if false\n            der(y)\n        else \n            der(z)\n        end\n    end"
    =#
    der = showDifferentiate(:(y = time))
    @test der == "der(y) = 1"

    der = showDifferentiate(:(y = a*x), [:a])
    @test der == "der(y) = a * der(x)"
  end

  @testset "Symbolic" begin

    @testset "Differentiate" begin
      testDifferentiate()
    end

  end

end

max_stars_repo_name: ModiaSim/ModiaBase

max_issues_repo_name: ModiaSim/ModiaBase

max_forks_repo_name: ModiaSim/ModiaBase


text:
# Box constrained problems

using Test, LinearAlgebra, SparseArrays, Random
using COSMO

@testset "Box Problems - feasible" begin
    A = sparse([1. 0; 0. 1.])
    b = [0.;0.]

    P = Matrix(Diagonal(ones(2)))
    q = [1.;-1.]
    l = [0.;0]
    u = [1.;1.]
    constraint = COSMO.Constraint(A,b, COSMO.Box(l,u))
    constraints = [constraint]

    model = COSMO.Model()
    assemble!(model, P, q, constraints)

    res = COSMO.optimize!(model)
    @test res.status == :Solved
    @test isapprox(res.obj_val, -0.5, atol=1e-5)
end
nothing

@testset "Box Problems - primal infeasible" begin
    A = sparse([1. 0; 1. 0.])
    b = [2.;0.]

    P = Matrix(Diagonal(ones(2)))
    q = [1.;-1.]
    l = [0.;0]
    u = [1.;1.]
    constraint = COSMO.Constraint(A,b, COSMO.Box(l,u))
    constraints = [constraint]

    model = COSMO.Model()
    assemble!(model, P, q, constraints)

    res = COSMO.optimize!(model)
    @test res.status == :Primal_infeasible
end

@testset "Box Problems - primal infeasible" begin
    A = sparse([1. 0; 1. 0.])
    b = [0.;0.]

    P = Matrix(Diagonal(ones(2)))
    q = [1.;-1.]
    l = [0.;2.]
    u = [1.;3.]
    constraint = COSMO.Constraint(A,b, COSMO.Box(l,u))
    constraints = [constraint]

    model = COSMO.Model()
    assemble!(model, P, q, constraints)

    res = COSMO.optimize!(model)
    @test res.status == :Primal_infeasible
end

@testset "Box Problems - dual infeasible (unscaled)" begin
    A = sparse([1. 0; 0. 1.])
    b = [1.;1.]

    P = Matrix(Diagonal(zeros(2)))
    q = [1.;1.]
    l = [0.;-Inf]
    u = [1.;3.]
    constraint = COSMO.Constraint(A,b, COSMO.Box(l,u))
    constraints = [constraint]

    model = COSMO.Model()
    settings = COSMO.Settings(check_infeasibility=20,scaling=0)
    assemble!(model, P, q, constraints,settings  = settings)

    res = COSMO.optimize!(model)
    @test res.status == :Dual_infeasible
end

@testset "Box Problems - dual infeasible (scaled)" begin
    A = sparse([1. 0; 0. 1.])
    b = [1.;1.]

    P = Matrix(Diagonal(zeros(2)))
    q = [1.;1.]
    l = [0.;-Inf]
    u = [1.;3.]
    constraint = COSMO.Constraint(A,b, COSMO.Box(l,u))
    constraints = [constraint]

    model = COSMO.Model()
    settings = COSMO.Settings(check_infeasibility=40,scaling=10)
    assemble!(model, P, q, constraints,settings = settings)

    res = COSMO.optimize!(model)
    @test res.status == :Dual_infeasible
end

nothing


text:
using PDMats: PDMat
using GaussianProcesses: make_posdef!

function chistat(gpT::GPE, gpC::GPE, Xb::AbstractMatrix)
    pred_T = predict_full(gpT, Xb)
    pred_C = predict_full(gpC, Xb)
    μ = pred_T[1].-pred_C[1]
    Σ = pred_T[2]+pred_C[2]
    return dot(μ, Σ \ μ)
end
function chistat(gpT::GPE, gpC::GPE, Xb::AbstractMatrix,
                 Σcliff::PDMat, cK_T::AbstractMatrix, cK_C::AbstractMatrix)
    μT = predict_mu(gpT, Xb, cK_T)
    μC = predict_mu(gpC, Xb, cK_C)
    μ = μT - μC
    return dot(μ, Σcliff \ μ)
end

function sim_chi_null!(
        gpT::GPE, gpC::GPE, gpNull::GPE,
        treat::BitVector, Xb::AbstractMatrix,
        Σcliff::PDMat, cK_T::AbstractMatrix, cK_C::AbstractMatrix)
    Ysim = prior_rand(gpNull)

    gpT.y = Ysim[treat]
    gpC.y = Ysim[.!treat]

    update_alpha!(gpT)
    update_alpha!(gpC)

    return chistat(gpT, gpC, Xb, Σcliff, cK_T, cK_C)
end
function nsim_chi(gpT::GPE, gpC::GPE, gpNull::GPE, Xb::AbstractMatrix, nsim::Int)
    gpT_mod = modifiable(gpT)
    gpC_mod = modifiable(gpC)

    _, Σcliff = cliff_face(gpT, gpC, Xb)
    cK_T = cov(gpT.kernel, Xb, gpT.x)
    cK_C = cov(gpC.kernel, Xb, gpC.x)

    treat = BitVector(undef, gpNull.nobs)
    treat[:] .= false
    treat[1:gpT.nobs] .= true
    kernel = gpT_mod.kernel
    mT = gpT_mod.mean
    mC = gpC_mod.mean
    Σraw, chol = make_posdef!(copy(Σcliff))
    PDcliff = PDMat(Σraw, chol)
    t_sims = [sim_chi_null!(gpT_mod, gpC_mod, gpNull, treat, Xb, PDcliff, cK_T, cK_C)
        for _ in 1:nsim];
    return t_sims
end

function boot_chi2test(gpT::GPE, gpC::GPE, Xb::AbstractMatrix, nsim::Int)
    gpNull = make_null(gpT, gpC)
    chi_obs = chistat(gpT, gpC, Xb)
    chi_sims = GeoRDD.nsim_chi(gpT, gpC, gpNull, Xb, nsim)
    return mean(chi_sims .> chi_obs)
end

function placebo_chi(angle::Float64, X::AbstractMatrix, Y::Vector,
                 kern::Kernel, mean::Mean, logNoise::Float64,
                 nsim::Int)
    shift = shift_for_even_split(angle, X)
    left = left_points(angle, shift, X)
    gp_left  = GPE(X[:,left],   Y[left], mean, kern, logNoise)
    gp_right = GPE(X[:,.!left], Y[.!left], mean, kern, logNoise)
    Xb = placebo_sentinels(angle, shift, X, 100)
    pval = boot_chi2test(gp_left, gp_right, Xb, nsim)
    return pval
end
function placebo_chi(angle::Float64, gp::GPE, nsim::Int; kwargs...)
    logNoise = convert(Float64, gp.logNoise)
    return placebo_chi(angle, gp.x, gp.y, gp.kernel, gp.mean, logNoise, nsim; kwargs...)
end


text:
neurons_all=Vector{Int64}(undef,0)
ee_all=Vector{Float64}(undef,0)
maetes=Vector{Float64}(undef,0)
for i=7000
    global neurons_all,ee_all,maetes
    for j=0.01:0.01:0.1
        tol1=copy(j)
        num_to_test=0
        yy_train=copy(yy_train_all)
        i1=abs.(yy_train.-num_to_test).<0.01;i2=abs.(yy_train.-num_to_test).>0.01;yy_train[i1].=1;;yy_train[i2].=0;
        mi=minimum(yy_train);yy_train.-=mi;ma=maximum(yy_train);yy_train./=ma;yy_train.*=1.0-2tol1;yy_train.+=tol1

        yy_test=copy(yy_test_all)
        i1=abs.(yy_test.-num_to_test).<0.01;i2=abs.(yy_test.-num_to_test).>0.01;yy_test[i1].=1;yy_test[i2].=0;
        yy_test.-=mi;yy_test./=ma;yy_test.*=1.0-2tol1;yy_test.+=tol1

        inds_all=1:i_train
        neurons=copy(i);items_per_neuron=(Int64(floor(i_train/(neurons))))*ones(Int64,neurons);n_per_part=[0;cumsum(items_per_neuron)];n_per_part[end]=i_train
        @time a_all,a_layer1,layer1=ANNBN.train_layer_1_sigmoid_fast(neurons,vars,i_train,n_per_part,inds_all,xx_train,yy_train)
        predl1=[layer1 ones(i_train)]*a_layer1
        maetr=mean(abs.(yy_train-predl1))
        maximum(abs.(yy_train-predl1))
        predl1_test,layer1_test=ANNBN.predict_new(a_all,a_layer1,xx_test,i_test,neurons);
        maete=mean(abs.(yy_test-predl1_test))
        maximum(abs.(yy_test-predl1_test))
        i1=predl1.<=0.5;predl1[i1].=tol1;i1=predl1.>0.5;predl1[i1].=1.0-tol1;i2=abs.(predl1.-yy_train).<1e-10;maetr=100sum(i2)/length(i2)
        i1=predl1_test.<=0.5;predl1_test[i1].=tol1;i1=predl1_test.>0.5;predl1_test[i1].=1.0-tol1;i2=abs.(predl1_test.-yy_test).<1e-10;
        maete=100sum(i2)/length(i2)

        println(i," ",tol1," ",maete)
        push!(neurons_all,i)
        push!(ee_all,tol1)
        push!(maetes,maete)

    end
end


text:
# Command-line interface for Pepper.jl

module CLI

using DocOpt, REPL.TerminalMenus
using Pepper

import Base: show

const logfile = Ref{Union{IO, Nothing}}()

log(args...) = logfile[] ≡ nothing || println(logfile[], args...)

function printlnlog(io, args...)
    println(io, args...)
    log(args...)
end

const doc = """
Usage:
    pepper [options] <criteria>

Options:
    -h, --help          Show this screen
    -d, --degree=<num>  Maximum number of criteria per question [default: 2]
    -t, --total=<num>   The maximum total score attainable
    -l, --logfile=<log> Log history of decisions
"""

const PROMPTS = [
    "They're equally good"
    "Alternative 1 is better"
    "Alternative 2 is better"
]

const MENU = RadioMenu(PROMPTS)

clearscr() = print("\x1B[2J\x1B[0f")

print_alternatives(criteria, levels, base, lhs, rhs; indent="") =
    print_alternatives(stdout, criteria, levels, base, lhs, rhs, indent=indent)

function print_alternatives(io::IO, criteria, levels, base, lhs, rhs; indent="")

    for (i, c) in enumerate([lhs, rhs])
        printlnlog(io, indent, "Alternative ", i, ":\n")
        for (criterion, level) in zip(base, c)
            printlnlog(io, indent, " "^3, criteria[criterion], ": ",
                    levels[criterion][level])
        end
        printlnlog(io)
    end

end

function ask_user(s, base, lhs, rhs)

    clearscr()

    print_alternatives(Pepper.criteria(s), Pepper.levels(s),
                       base, lhs, rhs)

    choice = request("Which is better?\n", MENU)
    if choice ∉ 1:3
        println()
        exit()
    end

    log("Choice: ", PROMPTS[choice])
    log()
    return choice

end

function parse(lines::Vector, delegate)
    criterion = nothing
    for line in lines
        line = strip(line)
        startswith(line, "#") && continue
        err() = error("malformed line: $line")
        if endswith(line, ":")
            criterion = rstrip(chop(line))
            Pepper.addcriterion(delegate, criterion)
        elseif startswith(line, "-")
            level = lstrip(chop(line, head=1, tail=0))
            pts = nothing
            parts = split(level, ":")
            1 ≤ length(parts) ≤ 2 || err()
            if length(parts) > 1
                level = strip(parts[1])
                pts = tryparse(Int, parts[2])
                pts ≡ nothing && err()
            end
            Pepper.addlevel(delegate, level)
            if pts ≢ nothing
                Pepper.setpts(delegate, criterion => level, pts)
            end
        elseif !isempty(line)
            err()
        end
    end
end

function show(io::IO, s::Pepper.System)
    first = true
    for (criterion, levels) in zip(Pepper.criteria(s), Pepper.levels(s))
        if first
            first = false
        else
            println()
        end
        print(io, criterion, ":")
        for level in levels
            println()
            print(io, " "^4, "- ", level, ": ",
                  something(Pepper.getpts(s, criterion=>level), "~"))
        end
    end
end

function cli()

    args = docopt(doc)

    s = Pepper.System()

    parse(readlines(args["<criteria>"]), s)

    if args["--logfile"] ≡ nothing
        logfile[] = nothing
    else
        logfile[] = open(args["--logfile"], "w")
    end

    if args["--total"] ≡ nothing
        total = nothing
    else
        total = Base.parse(Int, args["--total"])
    end

    Pepper.paprika(ask_user, s,
        maxdegree = Base.parse(Int, args["--degree"]),
        total = total)

    clearscr()

    println(s)

    println()

    close(logfile[])

end

__init__() = TerminalMenus.config(ctrl_c_interrupt=false)

end


text:
module Day07

using Statistics: median, mean

export day07

function readInput(inputfile="inputs/input07.txt")
    return parse.(Int64, split(readline(inputfile), ","))
end

function day07(inputdata = readInput("inputs/input07.txt"))
    return (day07Part1(inputdata), day07Part2(inputdata))
end

function day07Part1(values)
    mid = convert(Int64, median(values))
    return sum(abs.(values .- mid))
end

function day07Part2(values)
    mid = mean(values)
    return min( _total_fuel(values, round(Int64, mid, RoundDown)),
                _total_fuel(values, round(Int64, mid, RoundUp)))
end

function _total_fuel(values, point)
    Δs = abs.(values .- point)
    return sum(div.(Δs .* ( Δs .+ 1),  2))
end

end #module


text:
module BrokemiaHelperDashSpring

using ..Ahorn, Maple

@mapdef Entity "BrokemiaHelper/dashSpring" DashSpring(x::Integer, y::Integer, playerCanUse::Bool=true)
@mapdef Entity "BrokemiaHelper/dashSpringDown" DashSpringDown(x::Integer, y::Integer, playerCanUse::Bool=true)
@mapdef Entity "BrokemiaHelper/wallDashSpringRight" DashSpringRight(x::Integer, y::Integer)
@mapdef Entity "BrokemiaHelper/wallDashSpringLeft" DashSpringLeft(x::Integer, y::Integer)

const placements = Ahorn.PlacementDict(
    "Dash Spring (Up, Brokemia Helper)" => Ahorn.EntityPlacement(
        DashSpring
    ),
    "Dash Spring (Down, Brokemia Helper)" => Ahorn.EntityPlacement(
        DashSpringDown
    ),
    "Dash Spring (Left, Brokemia Helper)" => Ahorn.EntityPlacement(
        DashSpringRight
    ),
    "Dash Spring (Right, Brokemia Helper)" => Ahorn.EntityPlacement(
        DashSpringLeft
    ),
)

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

    return Ahorn.Rectangle(x - 6, y - 3, 12, 5)
end

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

    return Ahorn.Rectangle(x - 6, y, 12, 5)
end

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

    return Ahorn.Rectangle(x - 1, y - 6, 5, 12)
end

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

    return Ahorn.Rectangle(x - 4, y - 6, 5, 12)
end

sprite = "objects/BrokemiaHelper/dashSpring/00.png"

Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::DashSpring, room::Maple.Room) = Ahorn.drawSprite(ctx, sprite, 0, -8)
Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::DashSpringDown, room::Maple.Room) = Ahorn.drawSprite(ctx, sprite, 16, 24, rot=pi)
Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::DashSpringLeft, room::Maple.Room) = Ahorn.drawSprite(ctx, sprite, 24, 0, rot=pi / 2)
Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::DashSpringRight, room::Maple.Room) = Ahorn.drawSprite(ctx, sprite, -8, 16, rot=-pi / 2)

end

max_stars_repo_name: Brokemia/BrokemiaHelper

max_issues_repo_name: Brokemia/BrokemiaHelper

max_forks_repo_name: Brokemia/BrokemiaHelper


text:
import GitLab
using BenchmarkTools

myauth = GitLab.authenticate(ENV["GITLAB_AUTH"]) # don't hardcode your access tokens!
println("Authentication successful")
options = Dict("private_token" => myauth.token)

## Get the BenchmarkResults repo!
benchmark_results_repo = GitLab.repo_by_name("BenchmarkResults"; headers=options)

# CommentListener settings
trigger = r"`runbenchmarks\(.*?\)`"

# We can use Julia's `do` notation to set up the listener's handler function.
# Note that, in our example case
listener = GitLab.CommentListener(trigger; auth = myauth) do event, phrase
    ## pkg_names = matchall(r"\".*?\"", phrase.match)
    ## pkg_name = pkg_names[1]
    ## Take the current package !!
    pkg_name = event.payload["project"]["name"]
    ## @show pkg_name

    ## @show event.payload["object_attributes"]
    # Parse the original comment event for all the necessary reply info
    comment = GitLab.Comment(event.payload["object_attributes"])
    ## @show comment

    results = include(Pkg.dir(pkg_name, "benchmarks", "runbenchmarks.jl"))

    ## Send the results to the BenchmarkResults repository and get the link for the same.
    issue = GitLab.create_issue(benchmark_results_repo; headers=options, params = Dict("title" => "Benchmark Results"))
    benchmark_reply_to = get(issue.id)

    str_buf = IOBuffer(true, true)
    print(str_buf, get(benchmark_results_repo.web_url))
    issue_url = UTF8String(str_buf.data) * "/issues/$(get(issue.id))"

    ## @show event.payload["object_attributes"]
    if event.payload["object_attributes"]["noteable_type"] == "Issue"
        comment_params = Dict("body" => "
Your benchmark results are available here - $(issue_url)
")
        comment_kind = :issue
        reply_to = event.payload["object_attributes"]["noteable_id"]
    elseif event.payload["object_attributes"]["noteable_type"] == "Commit"
        comment_params = Dict("note" => "
Your benchmark results are available here - $(issue_url)
")
        comment_kind = :commit
        reply_to = get(comment.commit_id)
    elseif event.payload["object_attributes"]["noteable_type"] == "MergeRequest"
        comment_params = Dict("body" => "
Your benchmark results are available here - $(issue_url)
")
        comment_kind = :review
        reply_to = event.payload["object_attributes"]["noteable_id"]
        # load required query params for review comment creation
        comment_params["commit_id"] = "$(comment.id)"
        comment_params["path"] = "$(comment.url)"
        comment_params["position"] = "$(comment.id)"
    end
    ## @show reply_to, comment_kind

    ## @show benchmark_comment_params
    ## Update the issue with results
    benchmark_comment_params = Dict("body" => "
Your benchmark results are available !
 $results")
    GitLab.create_comment(benchmark_results_repo, benchmark_reply_to, :issue; headers = options, params = benchmark_comment_params)

    # send the comment creation request to GitLab
    GitLab.create_comment(event.repository, reply_to, comment_kind; headers = options, params = comment_params)

    return HttpCommon.Response(200)
end

# Start the listener on localhost at port 8000
GitLab.run(listener, host=IPv4(0,0,0,0), port=8000)


text:
import ODINN

processes = 16
python_path = "/home/jovyan/.conda/envs/oggm_env/bin/python3.9"
ODINN.initialize(processes, python_path)

function run()
    # Configure OGGM settings in all workers
    @everywhere oggm_config()

    ###############################################################
    ###########################  MAIN #############################
    ###############################################################

    # Defining glaciers to be modelled with RGI IDs
    rgi_ids = ["RGI60-11.03638", "RGI60-11.01450", "RGI60-08.00213", "RGI60-04.04351", "RGI60-01.02170",
                "RGI60-02.05098", "RGI60-01.01104", "RGI60-01.09162", "RGI60-01.00570", "RGI60-04.07051",
                "RGI60-07.00274", "RGI60-07.01323", "RGI60-03.04207", "RGI60-03.03533", "RGI60-01.17316"]

    ### Initialize glacier directory to obtain DEM and ice thickness inversion  ###
    gdirs = init_gdirs(rgi_ids)

    #########################################
    ###########  CLIMATE DATA  ##############
    #########################################

    # Process climate data for glaciers
    gdirs_climate = get_gdirs_with_climate(gdirs, true)

    # Run forward model for selected glaciers
    if create_ref_dataset
        println("Generating reference dataset for training...")

        # Compute reference dataset in parallel
        @everywhere solver = Ralston()
        H_refs, V̄x_refs, V̄y_refs = @time generate_ref_dataset(gdirs_climate)

        println("Saving reference data")
        jldsave(joinpath(root_dir, "data/PDE_refs.jld2"); H_refs, V̄x_refs, V̄y_refs)
    end

    # Load stored PDE reference datasets
    PDE_refs = load(joinpath(root_dir, "data/PDE_refs.jld2"))
    # Plot training dataset of glaciers
    plot_glacier_dataset(gdirs_climate, PDE_refs)

    #######################################################################################################
    #############################             Train UDE            ########################################
    #######################################################################################################

    # Training setup
    global current_epoch = 1

    # Train iceflow UDE in parallel
    # First train with ADAM to move the parameters into a favourable space
    @everywhere solver = ROCK4()

    if retrain
        println("Retraining from previous NN weights...")
        trained_weights = load(joinpath(root_dir, "data/trained_weights.jld2"))
        current_epoch = trained_weights["current_epoch"]
        θ_trained = trained_weights["θ_trained"]
        # train_settings = (BFGS(initial_stepnorm=0.05), 20) # optimizer, epochs
        train_settings = (ADAM(0.01), 20) # optimizer, epochs
        iceflow_trained, UA = @time train_iceflow_UDE(gdirs_climate, train_settings, PDE_refs, θ_trained) # retrain
        θ_trained = iceflow_trained.minimizer

        # Save trained NN weights
        println("Saving NN weights...")
        jldsave(joinpath(root_dir, "data/trained_weights.jld2"); θ_trained, current_epoch)
    else
        println("Training from scratch...")
        train_settings = (ADAM(0.05), 10) # optimizer, epochs
        iceflow_trained, UA = @time train_iceflow_UDE(gdirs_climate, train_settings, PDE_refs)
        θ_trained = iceflow_trained.minimizer
        println("Saving NN weights...")
        jldsave(joinpath(root_dir, "data/trained_weights.jld2"); θ_trained, current_epoch)

        # Continue training with BFGS
        # train_settings = (BFGS(initial_stepnorm=0.02f0), 20) # optimizer, epochs
        train_settings = (ADAM(0.02), 20) # optimizer, epochs
        iceflow_trained, UA = @time train_iceflow_UDE(gdirs_climate, train_settings, PDE_refs, θ_trained) # retrain
        θ_trained = iceflow_trained.minimizer
        # Save trained NN weights
        println("Saving NN weights...")
        jldsave(joinpath(root_dir, "data/trained_weights.jld2"); θ_trained, current_epoch)
    end

    ##########################################
    ####  Plot the final trained model  ######
    ##########################################
    data_range = -20.0:0.0
    pred_A = predict_A̅(UA, θ_trained, collect(data_range)')
    pred_A = [pred_A...] # flatten
    true_A = A_fake(data_range)

    Plots.scatter(true_A, label="True A")
    train_final = Plots.plot!(pred_A, label="Predicted A")
    Plots.savefig(train_final,joinpath(root_plots,"training","final_model.png"))
end

# Run main
run()


text:
function writemesh(mesh,filename,ifile)

nbe = size(mesh.eblks,2);
nbf = size(mesh.fblks,2);
neb = maximum(eblks[2,:] .- eblks[1,:])+1;
nfb = maximum(fblks[2,:] .- fblks[1,:])+1;
# nfe = size(mesh.t2f,1);
# nv = maximum(mesh.t[:]));

ndims = zeros(20,1);
ndims[1] = mesh.nd;
ndims[2] = mesh.ne;
ndims[3] = mesh.nf;
ndims[4] = mesh.nv;
ndims[5] = mesh.nfe;
ndims[6] = nbe;
ndims[7] = neb;
ndims[8] = nbf;
ndims[9] = nfb;

nsize = zeros(30,1);
nsize[1] = length(ndims);
nsize[2] = length(mesh.facecon);
nsize[3] = length(mesh.eblks);
nsize[4] = length(mesh.fblks);
nsize[5] = length(mesh.nbsd);
nsize[6] = length(mesh.elemsend);
nsize[7] = length(mesh.elemrecv);
nsize[8] = length(mesh.elemsendpts);
nsize[9] = length(mesh.elemrecvpts);
nsize[10] = length(mesh.elempart);
nsize[11] = length(mesh.elempartpts);
nsize[12] = length(mesh.cgelcon);
nsize[13] = length(mesh.rowent2elem);
nsize[14] = length(mesh.cgent2dgent);
nsize[15] = length(mesh.colent2elem);
nsize[16] = length(mesh.rowe2f1);
nsize[17] = length(mesh.cole2f1);
nsize[18] = length(mesh.ent2ind1);
nsize[19] = length(mesh.rowe2f2);
nsize[20] = length(mesh.cole2f2);
nsize[21] = length(mesh.ent2ind2);

if ifile>0
    display(["Writing mesh into file " * string(ifile)]);
    fileID = open([filename * string(ifile) * ".bin"],"w");
else
    display(["Writing mesh into file "]);
    fileID = open([filename * ".bin"],"w");
end

write(fileID,Float64(length(nsize[:])));
write(fileID,Float64.(nsize[:]));

if nsize[1]>0
    write(fileID,Float64.(ndims[:]));
end
#reshape(permute(mesh.facecon,[2 1 3]),[2*master.npf(1)*ndims[3) 1])
if nsize[2]>0
    write(fileID,Float64.(mesh.facecon[:]));
end
if nsize[3]>0
    write(fileID,Float64.(mesh.eblks[:]));
end
if nsize[4]>0
    write(fileID,Float64.(mesh.fblks[:]));
end
if nsize[5]>0
    write(fileID,Float64.(mesh.nbsd[:]));
end
if nsize[6]>0
    write(fileID,Float64.(mesh.elemsend[:]));
end
if nsize[7]>0
    write(fileID,Float64.(mesh.elemrecv[:]));
end
if nsize[8]>0
    write(fileID,Float64.(mesh.elemsendpts[:]));
end
if nsize[9]>0
    write(fileID,Float64.(mesh.elemrecvpts[:]));
end
if nsize[10]>0
    write(fileID,Float64.(mesh.elempart[:]));
end
if nsize[11]>0
    write(fileID,Float64.(mesh.elempartpts[:]));
end
if nsize[12]>0
    write(fileID,Float64.(mesh.cgelcon[:]-1));
end
if nsize[13]>0
    write(fileID,Float64.(mesh.rowent2elem[:]));
end
if nsize[14]>0
    write(fileID,Float64.(mesh.cgent2dgent[:]-1));
end
if nsize[15]>0
    write(fileID,Float64.(mesh.colent2elem[:]-1));
end
if nsize[16]>0
    write(fileID,Float64.(mesh.rowe2f1[:]));
end
if nsize[17]>0
    write(fileID,Float64.(mesh.cole2f1[:]-1));
end
if nsize[18]>0
    write(fileID,Float64.(mesh.ent2ind1[:]-1));
end
if nsize[19]>0
    write(fileID,Float64.(mesh.rowe2f2[:]));
end
if nsize[20]>0
    write(fileID,Float64.(mesh.cole2f2[:]-1));
end
if nsize[21]>0
    write(fileID,Float64.(mesh.ent2ind2[:]-1));
end

close(fileID);

end

max_stars_repo_name: glwagner/Exasim

max_issues_repo_name: glwagner/Exasim

max_forks_repo_name: glwagner/Exasim


text:
## Generative Merging MDP model

const HARD_BRAKE = 1
const RELEASE = 7
const WRAP_AROUND_TOL = 2.0 # when vehicles reach the end of the main lane - WRAP_AROUND_TOL they are spawned at the beginning

"""
    AugScene
Driving scene augmented with information about the ego vehicle
"""
struct AugScene
    scene::Scene
    ego_info::NamedTuple{(:acc,), Tuple{Float64}}
end

# Dummy type since the action is given by the policy
mutable struct EgoDriver{A} <: DriverModel{A}
    a::A
end

Base.rand(model::EgoDriver) = model.a
AutomotiveDrivingModels.observe!(m::EgoDriver, s::EntityFrame{S,D,I}, roadway::R, egoid::Int64) where {S,D,I,R} = m

"""
    GenerativeMergingMDP

A simulation environment for a highway merging scenario. Implemented using POMDPs.jl

# Parameters
    - `env::MergingEnvironment = MergingEnvironment(main_lane_angle = 0.0, merge_lane_angle = pi/7)`
    - `n_cars_main::Int64 = 1`
    - `n_cars_merge::Int64 = 1`
    - `n_agents::Int64 = n_cars_main + n_cars_merge`
    - `max_cars::Int64 = 16`
    - `min_cars::Int64 = 0`
    - `car_def::VehicleDef = VehicleDef()`
    - `dt::Float64 = 0.5 # time step`
    - `jerk_levels::SVector{5, Float64} = SVector(-1, -0.5, 0, 0.5, 1.0)`
    - `accel_levels::SVector{6, Float64} = SVector(-4.0, -2.0, -1.0, 0.0, 1.0, 2.0)`
    - `max_deceleration::Float64 = -4.0`
    - `max_acceleration::Float64 = 3.5`
    - `comfortable_acceleration::Float64 = 2.0`
    - `discount_factor::Float64 = 0.95`
    - `ego_idm::IntelligentDriverModel = IntelligentDriverModel(σ=0.0, v_des=env.main_lane_vmax)`
    - `default_driver_model::DriverModel{LaneFollowingAccel} = IntelligentDriverModel(v_des=env.main_lane_vmax)`
    - `observe_cooperation::Bool = false`
    - `observe_speed::Bool = true`
    - `traffic_speed::Symbol = :mixed`
    - `random_n_cars::Bool = false`
    - `driver_type::Symbol = :random` can be :binary, or :aggressive, or :cooperative
    - `max_burn_in::Int64 = 20`
    - `min_burn_in::Int64 = 10`
    - `initial_ego_velocity::Float64 = 10.0`
    - `initial_velocity::Float64 = 5.0`
    - `initial_velocity_std::Float64 = 1.0`
    - `main_lane_slots::LinRange{Float64} = LinRange(0.0, env.main_lane_length + env.after_merge_length, max_cars)`
    - `collision_cost::Float64 = -1.0`
    - `goal_reward::Float64 = 1.0`
    - `hard_brake_cost::Float64 = 0.0`
"""
@with_kw mutable struct GenerativeMergingMDP <: MDP{AugScene, Int64}
    env::MergingEnvironment = MergingEnvironment(main_lane_angle = 0.0, merge_lane_angle = pi/7)
    n_cars_main::Int64 = 1
    n_cars_merge::Int64 = 1
    n_agents::Int64 = n_cars_main + n_cars_merge
    max_cars::Int64 = 16
    min_cars::Int64 = 0
    car_def::VehicleDef = VehicleDef()
    dt::Float64 = 0.5 # time step
    jerk_levels::SVector{5, Float64} = SVector(-1, -0.5, 0, 0.5, 1.0)
    accel_levels::SVector{6, Float64} = SVector(-4.0, -2.0, -1.0, 0.0, 1.0, 2.0)
    max_deceleration::Float64 = -4.0
    max_acceleration::Float64 = 3.5
    comfortable_acceleration::Float64 = 2.0
    discount_factor::Float64 = 0.95
    ego_idm::IntelligentDriverModel = IntelligentDriverModel(σ=0.0, v_des=env.main_lane_vmax)
    default_driver_model::DriverModel{LaneFollowingAccel} = IntelligentDriverModel(v_des=env.main_lane_vmax)
    observe_cooperation::Bool = false
    observe_speed::Bool = true
    # initial state params
    traffic_speed::Symbol = :mixed
    random_n_cars::Bool = false
    driver_type::Symbol = :random
    max_burn_in::Int64 = 20
    min_burn_in::Int64 = 10
    initial_ego_velocity::Float64 = 10.0
    initial_velocity::Float64 = 5.0
    initial_velocity_std::Float64 = 1.0
    main_lane_slots::LinRange{Float64} = LinRange(0.0,
                                            env.main_lane_length + env.after_merge_length,
                                            max_cars)
    # reward params
    collision_cost::Float64 = -1.0
    goal_reward::Float64 = 1.0
    hard_brake_cost::Float64 = 0.0

    # internal states
    driver_models::Dict{Int64, DriverModel} = Dict{Int64, DriverModel}(EGO_ID=>EgoDriver(LaneFollowingAccel(0.0)))
end

POMDPs.discount(mdp::GenerativeMergingMDP) = mdp.discount_factor
POMDPs.actions(mdp::GenerativeMergingMDP) = 1:7
POMDPs.n_actions(mdp::GenerativeMergingMDP) = 7
POMDPs.actionindex(mdp::GenerativeMergingMDP, a::Int64) = a

function POMDPs.initialstate(mdp::GenerativeMergingMDP, rng::AbstractRNG)
    s0 = Scene()
    if mdp.random_n_cars
        mdp.n_cars_main = rand(rng, mdp.min_cars:mdp.max_cars)
    end
    mdp.driver_models = Dict{Int64, DriverModel}(EGO_ID=>EgoDriver(LaneFollowingAccel(0.0)))
    start_positions = sample(rng, mdp.main_lane_slots, mdp.n_cars_main, replace=false)

    start_velocities = mdp.initial_velocity .+ mdp.initial_velocity_std*randn(rng, mdp.n_cars_main)
    ego = initial_merge_car_state(mdp, rng, EGO_ID)
    ego_acc_0 = 0.0
    scene = s0
    for i=EGO_ID+1:EGO_ID+mdp.n_cars_main
        veh_state = vehicle_state(start_positions[i - EGO_ID], main_lane(mdp.env),
        start_velocities[i - EGO_ID], mdp.env.roadway)
        veh = Vehicle(veh_state, mdp.car_def, i)
        push!(s0, veh)
        mdp.driver_models[i] = CooperativeIDM()
        if mdp.traffic_speed == :mixed
            v_des = sample(rng, [4.0, 5., 6.0], Weights([0.2, 0.3, 0.5]))
        elseif mdp.traffic_speed == :fast
            v_des = 15.0
        end
        set_desired_speed!(mdp.driver_models[i], v_des)
        if mdp.driver_type == :random
            mdp.driver_models[i].c = rand(rng, 0:0.01:1) # change cooperativity
        elseif mdp.driver_type == :binary
            mdp.driver_models[i].c = sample(rng, [0,1], Weights([0.9, 0.1])) # change cooperativity
        elseif mdp.driver_type == :aggressive
            mdp.driver_models[i].c = 0.0
        elseif mdp.driver_type == :cooperative
            mdp.driver_models[i].c = 1.0
        end
    end
    # burn in
    acts = Vector{LaneFollowingAccel}(undef, length(scene))
    burn_in =rand(rng, mdp.min_burn_in:mdp.max_burn_in)
    scene = s0
    for t=1:burn_in
        get_actions!(acts, scene, mdp.env.roadway, mdp.driver_models)
        tick!(scene, mdp.env.roadway, acts, mdp.dt, true)
        for (i, veh) in enumerate(scene)
            scene[i] = wrap_around(mdp.env, scene[i])
        end
    end
    s = Scene()
    push!(s, ego)
    for veh in scene
        push!(s, veh)
    end
    return AugScene(s, (acc=ego_acc_0,))
end

function POMDPs.reward(mdp::GenerativeMergingMDP, s::AugScene, a::Int64, sp::AugScene)
    egop = get_by_id(sp.scene, EGO_ID)
    r = 0.0
    if reachgoal(mdp, egop)
       r += mdp.goal_reward
    elseif collision_checker(sp.scene, EGO_ID)
        r += mdp.collision_cost
    end
    if caused_hard_brake(mdp, sp.scene)
        r += mdp.hard_brake_cost
    end
    return r
end

function POMDPs.reward(pomdp::FullyObservablePOMDP{AugScene,Int64}, s::AugScene, a::Int64, sp::AugScene)
    return reward(pomdp.mdp, s, a, sp)
end

function POMDPs.isterminal(mdp::GenerativeMergingMDP, s::AugScene)
    return collision_checker(s.scene, EGO_ID) || reachgoal(mdp, get_by_id(s.scene, EGO_ID))
end

function POMDPs.gen(mdp::GenerativeMergingMDP, s::AugScene, a::Int64, rng::AbstractRNG)
    scene = deepcopy(s.scene)
    mdp.driver_models[EGO_ID].a = action_map(mdp, s.ego_info.acc, a)
    ego_acc = mdp.driver_models[EGO_ID].a.a
    acts = Vector{LaneFollowingAccel}(undef, length(scene))

    # call driver models
    for i=EGO_ID+1:EGO_ID+mdp.n_cars_main
        mdp.driver_models[i].other_acc = s.ego_info.acc
    end
    get_actions!(acts, scene, mdp.env.roadway, mdp.driver_models)

    # update scene
    tick!(scene, mdp.env.roadway, acts, mdp.dt, true)

    # clamp speed
    for (i, veh) in enumerate(scene)
        # scene[i] = clamp_speed(mdp.env, veh)
        scene[i] = wrap_around(mdp.env, scene[i])
    end

    return (sp=AugScene(scene, (acc=ego_acc,)), )
end

"""
    extract_features(mdp::GenerativeMergingMDP, s::AugScene)
extract a feature vector from AugScene
"""
function extract_features(mdp::GenerativeMergingMDP, s::AugScene)
    # @show keys(mdp.driver_models)
    scene = s.scene
    env = mdp.env
    features = -3*ones(15)
    ego_ind = findfirst(EGO_ID, scene)
    ego = scene[ego_ind]

    # distance to merge point
    s_ego = dist_to_merge(env, ego)
    v_ego = ego.state.v
    features[1] = s_ego
    features[2] = v_ego
    features[3] = s.ego_info.acc
    # get neighbors
    fore, merge, fore_main, rear_main = get_neighbors(env, scene, EGO_ID)

    # front neighbor
    features[6] = 0.0
    if fore.ind != nothing
        fore_id = scene[fore.ind].id
        v_oth = scene[fore.ind].state.v
        headway = fore.Δs
        features[4] = headway
        if mdp.observe_speed
            features[5] = v_oth
        else
            features[5] = 0.0
        end
        if mdp.observe_cooperation
            features[6] = mdp.driver_models[fore_id].c
        else
            features[6] = 0.
        end
    end

    # two closest cars in main lane
    features[9] = 0.0
    if fore_main.ind != nothing
        fore_main_id = scene[fore_main.ind].id
        v_oth_main_fore = scene[fore_main.ind].state.v
        headway_main_fore = fore_main.Δs
        features[7] = headway_main_fore
        if mdp.observe_speed
            features[8] = v_oth_main_fore
        else
            features[8] = 0.0
        end
        if mdp.observe_cooperation
            features[9] = mdp.driver_models[fore_main_id].c
        else
            features[9] = 0.
        end
    end
    features[12] = 0.0
    if rear_main.ind != nothing
        rear_main_id = scene[rear_main.ind].id
        v_oth_main_rear = scene[rear_main.ind].state.v
        headway_main_rear = rear_main.Δs
        features[10] = headway_main_rear
        if mdp.observe_speed
            features[11] = v_oth_main_rear
        else
            features[11] = 0.0
        end
        if mdp.observe_cooperation
            features[12] = mdp.driver_models[rear_main_id].c
        else
            features[12] = 0.
        end
    end

    # rear car to the merge point
    features[15] = 0.0
    if merge.ind != nothing
        merge_id = scene[merge.ind].id
        v_oth = scene[merge.ind].state.v
        headway = merge.Δs
        features[13] = headway
        if mdp.observe_speed
            features[14] = v_oth
        else
            features[14] = 0.0
        end
        if mdp.observe_cooperation
            features[15] = mdp.driver_models[merge_id].c
        else
            features[15] = 0.
        end
    end
    return features
end

"""
    normalize_features!(mdp::GenerativeMergingMDP, features::Vector{Float64})
normalize a feature vector extracted from a scene
"""
function normalize_features!(mdp::GenerativeMergingMDP, features::Vector{Float64})
    features[1] /=  mdp.env.main_lane_length
    features[2] /=  mdp.env.main_lane_vmax
    features[3] /=   mdp.max_deceleration
    features[4] /= mdp.env.main_lane_length
    features[5] /=  mdp.env.main_lane_vmax
    features[6] /= 1.
    features[7] /= mdp.env.main_lane_length
    features[8] /=   mdp.env.main_lane_vmax
    features[9] /= 1.
    features[10] /= mdp.env.main_lane_length
    features[11] /=  mdp.env.main_lane_vmax
    features[12] /= 1.
    features[13] /= mdp.env.main_lane_length
    features[14] /=  mdp.env.main_lane_vmax
    features[15] /= 1.

    return features
end

"""
     unnormalize_features!(mdp::GenerativeMergingMDP, features::Vector{Float64})
rescale feature vector
"""
function unnormalize_features!(mdp::GenerativeMergingMDP, features::Vector{Float64})
    features[1]  *=  mdp.env.main_lane_length
    features[2]  *=  mdp.env.main_lane_vmax
    features[3]  *=   mdp.max_deceleration
    features[4]  *= mdp.env.main_lane_length
    features[5]  *=  mdp.env.main_lane_vmax
    features[6]  *= 1.
    features[7]  *= mdp.env.main_lane_length
    features[8]  *=   mdp.env.main_lane_vmax
    features[9]  *= 1.
    features[10] *= mdp.env.main_lane_length
    features[11] *=  mdp.env.main_lane_vmax
    features[12] *= 1.
    features[13] *= mdp.env.main_lane_length
    features[14] *=  mdp.env.main_lane_vmax
    features[15] *= 1.

    return features
end

function POMDPs.convert_s(::Type{V}, s::AugScene, mdp::GenerativeMergingMDP) where V<:AbstractArray
    feature_vec = extract_features(mdp, s)
    return normalize_features!(mdp, feature_vec)
end

function POMDPs.convert_s(::Type{AugScene}, o::V, mdp::GenerativeMergingMDP) where V<:AbstractArray
    feature_vec = deepcopy(o)
    unnormalize_features!(mdp, feature_vec)
    s_ego, v_ego, s_front, v_front, s_fm, v_fm, s_rm, v_rm, s_m, v_m, a_ego = feature_vec[1:11]

    # reconstruct the scene
    scene = Scene()
    if s_ego < 0.0
        lane_ego = merge_lane(mdp.env)
        s_ego = get_end(lane_ego) + s_ego
    else
        lane_ego = main_lane(mdp.env)
        s_ego = mdp.env.roadway[mdp.env.merge_index].s + s_ego
    end
    ego = Vehicle(vehicle_state(s_ego, lane_ego, v_ego, mdp.env.roadway), VehicleDef(), EGO_ID)
    push!(scene, ego)

    # front neighbor
    if lane_ego == main_lane(mdp.env)
        s_front = ego.state.posF.s + s_front
    else
        s_front = mdp.env.main_lane_length + s_front
    end
    front = Vehicle(vehicle_state(s_front, main_lane(mdp.env), v_front, mdp.env.roadway), VehicleDef(), EGO_ID+1)
    push!(scene, front)

    # front neighbor to projection
    proj_lane = main_lane(mdp.env)
    main_lane_proj = proj(ego.state.posG, proj_lane, mdp.env.roadway)
    s_main = proj_lane[main_lane_proj.curveproj.ind, mdp.env.roadway].s
    s_fm = s_main + s_fm
    frontmain = Vehicle(vehicle_state(s_fm, main_lane(mdp.env), v_fm, mdp.env.roadway), VehicleDef(), EGO_ID+2)
    push!(scene, frontmain)

    s_rm = s_main - s_rm
    rearmain = Vehicle(vehicle_state(s_rm, main_lane(mdp.env), v_rm, mdp.env.roadway), VehicleDef(), EGO_ID+3)
    push!(scene, rearmain)

    if lane_ego == main_lane(mdp.env)
        s_m = ego.state.posF.s - s_m
    else
        s_m = mdp.env.main_lane_length - s_m
    end
    merge = Vehicle(vehicle_state(s_m, main_lane(mdp.env), v_m, mdp.env.roadway), VehicleDef(), EGO_ID+4)
    push!(scene, merge)
    return AugScene(scene, (acc=a_ego,))
end

## helpers

"""
    initial_merge_car_state(mdp::GenerativeMergingMDP, rng::AbstractRNG, id::Int64)
returns a Vehicle, at the initial state of the merging car.
"""
function initial_merge_car_state(mdp::GenerativeMergingMDP, rng::AbstractRNG, id::Int64)
    v0 = mdp.initial_velocity + mdp.initial_velocity_std*randn(rng)
    v0 = mdp.initial_ego_velocity
    veh_state = vehicle_state(0.0, merge_lane(mdp.env), v0, mdp.env.roadway)
    return Vehicle(veh_state, mdp.car_def, id)
end

"""
    reset_main_car_state(mdp::GenerativeMergingMDP, veh::Vehicle)
initialize a car at the beginning of the main lane
"""
function reset_main_car_state(mdp::GenerativeMergingMDP, veh::Vehicle)
    v0 = mdp.initial_velocity + mdp.initial_velocity_std*randn(rng)
    veh_state = vehicle_state(0.0, main_lane(mdp.env), v0, mdp.env.roadway)
    return Vehicle(veh_state, mdp.car_def, veh.id)
end

"""
    reachgoal(mdp::GenerativeMergingMDP, ego::Vehicle)
return true if `ego` reached the goal position
"""
function reachgoal(mdp::GenerativeMergingMDP, ego::Vehicle)
    lane = get_lane(mdp.env.roadway, ego)
    s = ego.state.posF.s
    return lane.tag == main_lane(mdp.env).tag && s >= get_end(lane)
end

"""
    caused_hard_brake(mdp::GenerativeMergingMDP, scene::Scene)
returns true if the ego vehicle caused its rear neighbor to hard brake
"""
function caused_hard_brake(mdp::GenerativeMergingMDP, scene::Scene)
    ego_ind = findfirst(EGO_ID, scene)
    fore_res = get_neighbor_rear_along_lane(scene, ego_ind, mdp.env.roadway)
    if fore_res.ind == nothing
        return false
    else
        return mdp.driver_models[fore_res.ind].a <= mdp.driver_models[fore_res.ind].idm.d_max
    end
end

"""
    action_map(mdp::GenerativeMergingMDP, acc::Float64, a::Int64)
maps integer to a LaneFollowingAccel
"""
function action_map(mdp::GenerativeMergingMDP, acc::Float64, a::Int64)
    if a == HARD_BRAKE
        return LaneFollowingAccel(mdp.max_deceleration)
    elseif a == RELEASE
        return LaneFollowingAccel(0.0)
    else
        return LaneFollowingAccel(clamp(acc + mdp.jerk_levels[a-1], mdp.max_deceleration, mdp.max_acceleration))
    end
end

"""
    vehicle_state(s::Float64, lane::Lane, v::Float64, roadway::Roadway)
convenient constructor for VehicleState
"""
function vehicle_state(s::Float64, lane::Lane, v::Float64, roadway::Roadway)
    posF = Frenet(lane, s)
    return VehicleState(posF, roadway, v)
end

"""
    wrap_around(env::MergingEnvironment, veh::Vehicle)
respawn vehicle at the beginning of the main lane
"""
function wrap_around(env::MergingEnvironment, veh::Vehicle)
    lane = get_lane(env.roadway, veh)
    s_end = get_end(lane)
    s = veh.state.posF.s
    if s >= s_end - WRAP_AROUND_TOL && lane == main_lane(env) && veh.id != EGO_ID
        veh_state = vehicle_state(0.0, main_lane(env), veh.state.v, env.roadway)
        return Vehicle(veh_state,veh.def, veh.id)
    end
    return veh
end

"""
    clamp_speed(env::MergingEnvironment, veh::Vehicle)
clamp the speed of `veh` between 0. and main lane vmax
"""
function clamp_speed(env::MergingEnvironment, veh::Vehicle)
    v = clamp(veh.state.v, 0.0, env.main_lane_vmax)
    vehstate = VehicleState(veh.state.posG, veh.state.posF, v)
    return Vehicle(vehstate, veh.def, veh.id)
end

"""
    spread_out_initialization(mdp::GenerativeMergingMDP, rng::AbstractRNG)
spread out vehicles on the main lane
"""
function spread_out_initialization(mdp::GenerativeMergingMDP, rng::AbstractRNG)
    start_positions = zeros(mdp.n_cars_main)
    start_positions[1] = rand(rng, mdp.main_lane_slots)
    gap_length = div(mdp.env.main_lane_length + mdp.env.after_merge_length, mdp.n_cars_main)
    main_roadway = StraightRoadway(mdp.env.main_lane_length + mdp.env.after_merge_length)
    for i=2:mdp.n_cars_main
        start_positions[i] = mod_position_to_roadway(start_positions[i-1] + gap_length, main_roadway)
    end
    return start_positions
end

"""
    global_features(mdp::GenerativeMergingMDP, s::AugScene)
extract a vector with the states of all the vehicles in AugScene
"""
function global_features(mdp::GenerativeMergingMDP, s::AugScene)
    n_features = 2*mdp.n_cars_main + 3 + mdp.n_cars_main
    features = zeros(n_features)
    ego = get_by_id(s.scene, EGO_ID)
    features[1] = dist_to_merge(mdp.env, ego)
    features[2] = ego.state.v
    features[3] = s.ego_info.acc
    @assert s.scene[1].id == EGO_ID
    for i=2:length(s.scene)
        veh = s.scene[i]
        features[3*i-2] = veh.state.posF.s
        features[3*i-1] = veh.state.v
        obs_c = 0.5
        if mdp.observe_cooperation
            obs_c = mdp.driver_models[i].c
        end
        features[3*i] = obs_c
    end
    return features
end

function normalize_global_features!(mdp::GenerativeMergingMDP, features::Vector{Float64})
    features[1] /= mdp.env.main_lane_length
    features[2] /= mdp.env.main_lane_vmax
    features[3] /= mdp.max_deceleration
    for i=2:mdp.n_cars_main+1
        features[3*i - 2] /= mdp.env.main_lane_length
        features[3*i - 1] /= mdp.env.main_lane_vmax
    end
    return features
end

function unnormalize_global_features!(mdp::GenerativeMergingMDP, features::Vector{Float64})
    features[1] *= mdp.env.main_lane_length
    features[2] *= mdp.env.main_lane_vmax
    features[3] *= mdp.max_deceleration
    for i=2:mdp.n_cars_main+1
        features[3*i-2] *= mdp.env.main_lane_length
        features[3*i-1] *= mdp.env.main_lane_vmax
    end
    return features
end


text:
#A whole bunch of splatting...
u(x; p...) = x^(1-p[:γ])/(1-p[:γ])
g(; a, γ) = a + γ #Splatting the exact named parameters!
function f(x; p...)
    x2 = g(; p...) + x * p[:a]
    return u(x2; p...)
end

#Trying
u2(x, γ) = x^(1-γ)/(1-γ)
function f2(x; p...)
    x2 = g(;p...) + x * p[:a]
    return u2(x2, p[:γ])
end
kw = (a=1, γ=0.5)
f(2.0;kw...)

@btime f(2.0;$kw...)
@btime f2(2.0;$kw...) #The same.

#With a structure and no types
struct Params
    a
    γ
end
param = Params(1, 0.5)

## There is something oddly inconsistent about splatting named tuples in v0.7
myp = (a=2,b=4)
myp.a #Note typical named tuple access
f(;p...) = p[:a] #Note access as a dictionary!
f(;myp...)
#Base.getproperty(p::Base.Iterators.Pairs{<:Any,<:Any,<:Any,<:NamedTuple}, prop) = getproperty(p.data, prop)
#The above used to work, but doesn't now...


text:
"""
    hausdorff_metric(set1, set2)
    hausdorff_metric(set1, set2, p)

Given two sets of points `set1` & `set2`, compute the Hausdorff distance between the two sets.
By default the returned Hausdorff is the mean Hausdorff. If given a parameter `p`, Returns a
tuple for the percentile of each set.

# Arguments
- set1: set of points in 2D or 3D space
- set2: set of points in 2D or 3D space
- p: percentile to be computed based on `StatsBase.percentile`
"""
function hausdorff_metric(set1, set2)
    min_euc_list_u = []
    min_euc_list_v = []

    # Loop through every point in `set1` and find its corresponding closest point to `set2`
    for points1 in set1
        euc_list_1 = []
        for points2 in set2
            euclidean_dist = euc(points1, points2)
            append!(euc_list_1, euclidean_dist)
        end
        append!(min_euc_list_u, minimum(euc_list_1))
    end

    # Loop through every point in `set2` and find its corresponding closest point to `set1`
    for points1 in set2
        euc_list_1 = []
        for points2 in set1
            euclidean_dist = euc(points1, points2)
            append!(euc_list_1, euclidean_dist)
        end
        append!(min_euc_list_v, minimum(euc_list_1))
    end

    # Take the mean of each of these points to return the mean Hausdorff distance
    return Statistics.mean([
        Statistics.mean(min_euc_list_u), Statistics.mean(min_euc_list_v)
    ])
end

function hausdorff_metric(set1, set2, p)
    min_euc_list_u = []
    min_euc_list_v = []

    # Loop through every point in `set1` and find its corresponding closest point to `set2`
    for points1 in set1
        euc_list_1 = []
        for points2 in set2
            euclidean_dist = euc(points1, points2)
            append!(euc_list_1, euclidean_dist)
        end
        append!(min_euc_list_u, minimum(euc_list_1))
    end

    # Loop through every point in `set2` and find its corresponding closest point to `set1`
    for points1 in set2
        euc_list_1 = []
        for points2 in set1
            euclidean_dist = euc(points1, points2)
            append!(euc_list_1, euclidean_dist)
        end
        append!(min_euc_list_v, minimum(euc_list_1))
    end

    return (
        StatsBase.percentile(min_euc_list_u, p), StatsBase.percentile(min_euc_list_v, p)
    )
end

"""

    mean_hausdorff_2D(u, v, d, f)
Extract the edges of two images and then find the average Hausdorff distance along those edges.
Both arrays u and v are required to be the same size and they must be binary or boolean images.
The argument d corresponds to the distance metric used for computing. Typically, Euclidean() is
a common metric for use in Hausdorff distance computation. More options can be found in
Distances.jl. The argument `f` is the algorithm used for edge detection. More options can be
found in ImageEdgeDetection.jl.
"""

function mean_hausdorff_2D(u, v, d, f)
    edges_1 = detect_edges(u, f)
    edges_2 = detect_edges(v, f)

    sz_1 = size(edges_1, 1)
    sz_2 = size(edges_2, 1)

    # Loop through every point on `edge_1` and find its corresponding closest point to `edge_2`
    min_u =
        minimum.([
            map(x -> evaluate(d, edges_1[i, :], edges_2[x, :]), 1:sz_2) for i in 1:sz_1
        ])

    # Loop through every point on `edge_2` and find its corresponding closest point to `edge_1`
    min_v =
        minimum.([
            map(x -> evaluate(d, edges_2[i, :], edges_1[x, :]), 1:sz_1) for i in 1:sz_2
        ])

    # Take the average of each of these points to return the average Hausdorff distance
    return mean([mean(min_u), mean(min_v)])
end


text:
TEST_TIME = false
if !isdefined(Main,:TEST_TIME)
    TEST_TIME = true
end

# Double Integrator
solver = ALTROSolver(Problems.DoubleIntegrator()...)
b = benchmark_solve!(solver)
TEST_TIME && @test minimum(b).time / 1e6 < 0.5
@test max_violation(solver) < 1e-6
@test iterations(solver) <= 8 # 8

# Pendulum
solver = ALTROSolver(Problems.Pendulum()...)
b = benchmark_solve!(solver)
TEST_TIME && @test minimum(b).time / 1e6 < 1
@test max_violation(solver) < 1e-6
@test iterations(solver) <= 19 # 19

# Cartpole
solver = ALTROSolver(Problems.Cartpole()...)
b = benchmark_solve!(solver)
TEST_TIME && @test minimum(b).time / 1e6 < 5
@test max_violation(solver) < 1e-6
@test iterations(solver) <= 40 # 40

# Acrobot
solver = ALTROSolver(Problems.Acrobot()...)
b = benchmark_solve!(solver)
TEST_TIME && @test minimum(b).time / 1e6 < 10
@test max_violation(solver) < 1e-6
@test iterations(solver) <= 50 # 50

# Parallel Park
solver = ALTROSolver(Problems.DubinsCar(:parallel_park)...)
b =  benchmark_solve!(solver)
TEST_TIME && @test minimum(b).time /1e6 < 10
@test max_violation(solver) < 1e-6
@test iterations(solver) <= 13 # 13

# Three Obstacles
solver = ALTROSolver(Problems.DubinsCar(:three_obstacles)...)
b = benchmark_solve!(solver)
TEST_TIME && @test minimum(b).time /1e6 < 10
@test max_violation(solver) < 1e-6
@test iterations(solver) <= 20 # 20

# Escape
solver = ALTROSolver(Problems.DubinsCar(:escape)..., infeasible=true, R_inf=0.1)
b = benchmark_solve!(solver)
TEST_TIME && @test minimum(b).time / 1e6 < 20
@test max_violation(solver) < 1e-6
@test iterations(solver) <= 13 # 13

# Zig-zag
solver = ALTROSolver(Problems.Quadrotor(:zigzag)...)
b = benchmark_solve!(solver)
TEST_TIME && @test minimum(b).time / 1e6 < 50
@test max_violation(solver) < 1e-6
@test iterations(solver) <= 15 # 14

# Barrell Roll
solver = ALTROSolver(Problems.YakProblems()...)
b = benchmark_solve!(solver)
TEST_TIME && @test minimum(b).time / 1e6 < 100
@test max_violation(solver) < 1e-6
@test iterations(solver) <= 18 # 18


text:
using Core.Compiler: AbstractInterpreter, NativeInterpreter, InferenceState,
    OptimizationState, CodeInfo, CodeInstance, InferenceResult, WorldRange,
    IRCode, SSAValue, inlining_policy

struct InferredSource
    src::CodeInfo
    stmt_info::Vector{Any}
    rt
    InferredSource(src::CodeInfo, stmt_info::Vector{Any}, @nospecialize(rt)) =
        new(src, stmt_info, rt)
end

struct OptimizedSource
    ir::IRCode
    src::CodeInfo
    isinlineable::Bool
end

const Remarks = Vector{Pair{Int, String}}

struct CthulhuInterpreter <: AbstractInterpreter
    native::AbstractInterpreter

    unopt::Dict{Union{MethodInstance, InferenceResult}, InferredSource}
    opt::Dict{MethodInstance, CodeInstance}

    remarks::Dict{MethodInstance, Remarks}
end

CthulhuInterpreter(interp::AbstractInterpreter=NativeInterpreter()) =
    CthulhuInterpreter(
        interp,
        Dict{MethodInstance, InferredSource}(),
        Dict{MethodInstance, CodeInstance}(),
        Dict{MethodInstance, Remarks}()
    )

import Core.Compiler: InferenceParams, OptimizationParams, get_world_counter,
    get_inference_cache, code_cache, method_table,
    WorldView, lock_mi_inference, unlock_mi_inference, InferenceState
using Base: @invoke

Compiler.InferenceParams(interp::CthulhuInterpreter) = InferenceParams(interp.native)
Compiler.OptimizationParams(interp::CthulhuInterpreter) = OptimizationParams(interp.native)
Compiler.get_world_counter(interp::CthulhuInterpreter) = get_world_counter(interp.native)
Compiler.get_inference_cache(interp::CthulhuInterpreter) = get_inference_cache(interp.native)
Compiler.method_table(interp::CthulhuInterpreter, sv::InferenceState) = method_table(interp.native, sv)

# No need to do any locking since we're not putting our results into the runtime cache
Compiler.lock_mi_inference(interp::CthulhuInterpreter, mi::MethodInstance) = nothing
Compiler.unlock_mi_inference(interp::CthulhuInterpreter, mi::MethodInstance) = nothing

struct CthulhuCache
    cache::Dict{MethodInstance, CodeInstance}
end

Compiler.code_cache(interp::CthulhuInterpreter) = WorldView(CthulhuCache(interp.opt), WorldRange(get_world_counter(interp)))
Compiler.get(wvc::WorldView{CthulhuCache}, mi::MethodInstance, default) = get(wvc.cache.cache, mi, default)
Compiler.haskey(wvc::WorldView{CthulhuCache}, mi::MethodInstance) = haskey(wvc.cache.cache, mi)
Compiler.setindex!(wvc::WorldView{CthulhuCache}, ci::CodeInstance, mi::MethodInstance) = setindex!(wvc.cache.cache, ci, mi)

Compiler.may_optimize(interp::CthulhuInterpreter) = true
Compiler.may_compress(interp::CthulhuInterpreter) = false
Compiler.may_discard_trees(interp::CthulhuInterpreter) = false
Compiler.verbose_stmt_info(interp::CthulhuInterpreter) = true

function Compiler.add_remark!(interp::CthulhuInterpreter, sv::InferenceState, msg)
    if !haskey(interp.remarks, sv.linfo)
        interp.remarks[sv.linfo] = Remarks()
    end
    push!(interp.remarks[sv.linfo], sv.currpc => msg)
end

function Compiler.finish(state::InferenceState, interp::CthulhuInterpreter)
    r = @invoke Compiler.finish(state::InferenceState, interp::AbstractInterpreter)
    interp.unopt[Core.Compiler.any(state.result.overridden_by_const) ? state.result : state.linfo] = InferredSource(
        copy(state.src),
        copy(state.stmt_info),
        state.result.result)
    return r
end

function Compiler.transform_result_for_cache(interp::CthulhuInterpreter, linfo::MethodInstance,
        valid_worlds::WorldRange, @nospecialize(inferred_result))
    if isa(inferred_result, OptimizationState)
        opt = inferred_result
        ir = opt.ir
        if ir !== nothing
            return OptimizedSource(ir, opt.src, opt.src.inlineable)
        end
    end
    return inferred_result
end

# branch on https://github.com/JuliaLang/julia/pull/41328
@static if isdefined(Compiler, :is_stmt_inline)
function Compiler.inlining_policy(
    interp::CthulhuInterpreter, @nospecialize(src), stmt_flag::UInt8,
    mi::MethodInstance, argtypes::Vector{Any})
    @assert isa(src, OptimizedSource) || isnothing(src)
    if isa(src, OptimizedSource)
        if Compiler.is_stmt_inline(stmt_flag) || src.isinlineable
            return src.ir
        end
    else
        # the default inlining policy may try additional effor to find the source in a local cache
        return @invoke Compiler.inlining_policy(
            interp::AbstractInterpreter, nothing, stmt_flag::UInt8,
            mi::MethodInstance, argtypes::Vector{Any})
    end
    return nothing
end
else # @static if isdefined(Compiler, :is_stmt_inline)
function Compiler.inlining_policy(interp::CthulhuInterpreter)
    function (src)
        @assert isa(src, OptimizedSource)
        return src.isinlineable ? src.ir : nothing
    end
end
end # @static if isdefined(Compiler, :is_stmt_inline)

function Compiler.finish!(interp::CthulhuInterpreter, caller::InferenceResult)
    src = caller.src
    if isa(src, OptimizationState)
        opt = src
        ir = opt.ir
        if ir !== nothing
            caller.src = OptimizedSource(ir, opt.src, opt.src.inlineable)
        end
    end
end


text:
module TestWeights

using Test
import GraphIdx

@testset "Weights                        " begin
    @test GraphIdx.Ones <: GraphIdx.Weights
    @test GraphIdx.Const <: GraphIdx.Weights
    @test GraphIdx.Vec <: GraphIdx.Weights
end

@testset "Weights: is_const              " begin
    @test GraphIdx.is_const(GraphIdx.Ones)
    @test GraphIdx.is_const(GraphIdx.Ones{Int8})
    let c = GraphIdx.Weights(Float64)
        @test GraphIdx.is_const(c)
    end

    @test GraphIdx.is_const(GraphIdx.Const)
    @test GraphIdx.is_const(GraphIdx.Const{Float64})
    let c = GraphIdx.Weights(2.0)
        @test GraphIdx.is_const(c)
    end

    @test !GraphIdx.is_const(GraphIdx.Vec)
    @test !GraphIdx.is_const(GraphIdx.Vec{Float64})
    let c = GraphIdx.Weights([2.0, 1.0])
        @test !GraphIdx.is_const(c)
    end
end

end

max_stars_repo_name: EQt/graphidx

max_issues_repo_name: EQt/graphidx

max_forks_repo_name: EQt/graphidx