text:
@testset "Lexer tests" begin
    @testset "Old PList Lexing" begin
        l = lex_plist("{eggs = spam; foo = (bar, \"foo bar\");}")
        @test next_token(l) == Token(LBRACE, "{")
        @test next_token(l) == Token(IDENT, "eggs")
        @test next_token(l) == Token(EQUAL, "=")
        @test next_token(l) == Token(IDENT, "spam")
        @test next_token(l) == Token(SEMICOLON, ";")
        @test next_token(l) == Token(IDENT, "foo")
        @test next_token(l) == Token(EQUAL, "=")
        @test next_token(l) == Token(LPAREN, "(")
        @test next_token(l) == Token(IDENT, "bar")
        @test next_token(l) == Token(COMMA, ",")
        @test next_token(l) == Token(STRING, "foo bar")
        @test next_token(l) == Token(RPAREN, ")")
        @test next_token(l) == Token(SEMICOLON, ";")
        @test next_token(l) == Token(RBRACE, "}")
    end

    @testset "Numbers and Identifiers Lexing" begin
        l = lex_plist("one 2 three 3 five")
        @test next_token(l) == Token(IDENT, "one")
        @test next_token(l) == Token(NUMBER, "2")
        @test next_token(l) == Token(IDENT, "three")
        @test next_token(l) == Token(NUMBER, "3")
        @test next_token(l) == Token(IDENT, "five")
    end
end

@testset "Parser tests" begin
    @testset "Numbers and Identifiers Parser" begin
        l = lex_plist("{eggs = spam; foo = (bar, \"foo bar\");}")
        p = Parser(l)
        @test peek_token_type(p) == LBRACE
        next_token(p)
        @test peek_token_type(p) == IDENT
        next_token(p)
        @test peek_token_type(p) == EQUAL
        next_token(p)
        @test peek_token_type(p) == IDENT
        next_token(p)
        @test peek_token_type(p) == SEMICOLON
        next_token(p)
        @test peek_token_type(p) == IDENT
        next_token(p)
        @test peek_token_type(p) == EQUAL
        next_token(p)
        @test peek_token_type(p) == LPAREN
        next_token(p)
        @test peek_token_type(p) == IDENT
        next_token(p)
        @test peek_token_type(p) == COMMA
        next_token(p)
        @test peek_token_type(p) == STRING
        next_token(p)
        @test peek_token_type(p) == RPAREN
        next_token(p)
        @test peek_token_type(p) == SEMICOLON
        next_token(p)
        @test peek_token_type(p) == RBRACE
        next_token(p)
        @test peek_token_type(p) == EOF

        # strange_tokens = token_producer("\"Dot(u Vector2D) float64\" \"()\$0\"")
        # @test pop!(strange_tokens) == Token(STRING, "Dot(u Vector2D) float64")
        # @test pop!(strange_tokens) == Token(STRING, "()\$0")
    end

    @testset "PList decoding tests" begin
        @test readplist_string("1234") == 1234
        @test readplist_string("foobar") == "foobar"
        @test readplist_string("(1, 2, 3, 4)") == [1, 2, 3, 4]
        @test readplist_string("(one, two, three)") == ["one", "two", "three"]
        @test readplist_string("{foo = bar;}") == Dict("foo" => "bar")
        @test readplist_string("{one = 1;}") == Dict("one" => 1)
        @test readplist_string("<0fbd77 1c2735ae>") == UInt8[0x0f, 0xbd, 0x77, 0x1c, 0x27, 0x35, 0xae]
        @test readplist_string("<ffff>") == UInt8[0xff, 0xff]
        @test readplist_string("<00 01 02 03>") == UInt8[0x00, 0x01, 0x02, 0x03]
        @test readplist_string("{bin = <abcdef>;}") == Dict("bin" => UInt8[0xab, 0xcd, 0xef])

        dict = readplist("example.plist")

        @test !isempty(dict)
        @test haskey(dict, "Dogs")
        @test length(dict["Dogs"][1]) == 3
        @test dict["Dogs"][1]["Name"] == "Scooby Doo"
        @test dict["Dogs"][1]["Age"] == 43
        @test dict["Dogs"][1]["Colors"] == ["Brown", "Black"]
        @test dict["BinaryData"] == UInt8[0x0f, 0xbd, 0x77, 0x1c, 0x27, 0x35, 0xae]
    end
end

@testset "PList encoding tests" begin
    @test writeplist_string(Dict()) == "{}"
    @test writeplist_string([]) == "()"
    @test writeplist_string("foobar") == "foobar"
    @test writeplist_string("foo bar") == "\"foo bar\""
    @test writeplist_string(1234) == "1234"
    @test writeplist_string([1, 2, 3, 4]) == "(1, 2, 3, 4)"
    @test writeplist_string(["one", "two", "three"]) == "(one, two, three)"
    @test writeplist_string(["first number", "second number"]) == "(\"first number\", \"second number\")"
    @test writeplist_string(Dict("foo" => "bar")) == "{foo = bar;}"
    @test writeplist_string(Dict("one" => 1)) == "{one = 1;}"

    @test writeplist_string([1, [2, 3], 4]) == "(1, (2, 3), 4)"
    @test writeplist_string(Dict("colors" => ["red", "green", "blue"])) == "{colors = (red, green, blue);}"
end

@testset "XML parser tests" begin
    @testset "Test Element Node" begin
        doc = parsexml("<key>Author</key>", ignore_declaration=true)
        @test hasroot(doc)
        r = root(doc)
        @test nodename(r) == "key"
        n = first(nodes(r))
        @test istext(n)
        @test nodecontent(n) == "Author"
        @test countnodes(r) == 1
    end

    @testset "Test Element Nodes Hierarchy" begin
        doc = parsexml("<numbers><one>en</one><two>to</two><three>tre</three></numbers>", ignore_declaration=true)
        @test hasroot(doc)
        r = root(doc)
        @test nodename(r) == "numbers"
        @test nodecontent(r) == "entotre"
        @test countnodes(r) == 3

        @test nodecontent(nodes(r)[1]) == "en"
        @test nodecontent(nodes(r)[2]) == "to"
        @test nodecontent(nodes(r)[3]) == "tre"

        @test nodename(nodes(r)[1]) == "one"
        @test nodename(nodes(r)[2]) == "two"
        @test nodename(nodes(r)[3]) == "three"
    end

    @testset "Test Small Example" begin
        doc = parsexml("""
        <primates>
            <genus name="Homo">
                <species name="sapiens">Human</species>
            </genus>
            <genus name="Pan">
                <species name="paniscus">Bonobo</species>
                <species name="troglodytes">Chimpanzee</species>
            </genus>
        </primates>
        """, ignore_declaration=true)
        primates = root(doc)
        @test nodename(primates) == "primates"
        genus = nodes(primates)
        @test countnodes(primates) == 2
        @test nodecontent(genus[1]) == "Human"
        homo_genus_attr = attributes(genus[1])
        homo_attr = first(homo_genus_attr)
        @test homo_attr.name == "name"
        @test homo_attr.value == "Homo"
    end
end

@testset "XML PList tests" begin
    @testset "PList decoding tests" begin
        dict = read_xml_plist("example.plist.xml")

        @test !isempty(dict)
        @test haskey(dict, "Author")
        @test dict["Author"] == "William Shakespeare"
        @test haskey(dict, "Lines")
        @test length(dict["Lines"]) == 2
        @test dict["Lines"][1]  == "It is a tale told by an idiot,"
        @test dict["Lines"][2]  == "Full of sound and fury, signifying nothing."
        @test dict["Birthdate"] == 1564
    end

    @testset "Simple PList tests" begin
            dict = read_xml_plist_string(
            """
            <plist version="1.0">
            <dict>
                <key>egg</key>
                <string>spam</string>
                <key>numbers</key>
                <array>
                    <string>one</string>
                    <string>two</string>
                </array>
            </dict>
            </plist>""")
            @test dict["egg"] == "spam"
            @test dict["numbers"][1] == "one"
            @test dict["numbers"][2] == "two"
    end

    @testset "Test PList parsing independent from XML parsing" begin
        root = ElementNode("plist")
        root["version"] = "1.0"
        strings  = Node[ElementNode("string", "one"), ElementNode("string", "two")]
        children = Node[ElementNode("key", "egg"),
                        ElementNode("string", "spam"),
                        ElementNode("key", "numbers"),
                        ElementNode("array", strings)]
        addchild!(root, ElementNode("dict", children))

        dict = PLists.parse_obj(root)
        @test dict["egg"] == "spam"
        @test dict["numbers"][1] == "one"
        @test dict["numbers"][2] == "two"
    end
end


text:
#####################################################################
# IMPLEMENTACIÓN DEL ALGORITMOS DE INICIALIZACIÓN DE PESOS
# Basado en capítulo 7, algoritmo 6
#####################################################################

# Tamaño de la red neuronal y conjunto de datos

#Caso h:R -> R
include("single-input-single-output.jl")
include("utils.jl")
#Caso h:R^d -> R
include("multiple-input-single-ouput.jl")
#Caso h:R^d -> R^s con r,s> 1
include("multiple-input-multiple-output.jl")

ext: jl


text:
@testset "a simple arXiv API call" begin
	one_entry_bib = request("LaBollita",
				field=author,
				max_results=1,
				filename="test"
				)
	url  = "http://arxiv.org/abs/2008.02237v1"
	year = "2020"
	key  = "Krishna2020Effe"
	@test length(one_entry_bib) == 1
	@test one_entry_bib[1].url == url
	@test one_entry_bib[1].year == year
	@test one_entry_bib[1].key == key
	@test isfile("test.bib") == true

	if isfile("test.bib")
		rm("test.bib")
	end
end


text:
module MatrixChainMultiply

export matrixchainmultiply

### mcm = matrix chain multiply

"""
  `matrixchainmultiply(A,B,C,D,...)`

Algorithm for **Matrix-chain multiplication** (where matrices are also
generalizations of scalars and vectors). Using *Cormen* Ed. 3, p. 371.

Arguments
---------
`A,B,C,...` are matrices, vectors, or scalars of appropriate sizes.

"""
function matrixchainmultiply(A...)
  p = mcm_makep(A...)
  m, s, n = mcm_cost_matrices(p)
  mcm_compute(s, 1, n, A...)
end
"""
  `matrixchainmultiply(fn_name::String, A...)`

Print out the optimal matrix order as an expresion and its
cost. Also use this to...

Construc a function for repeated use (Julia v0.6-?)
---------------------------------------------------
Sometimes, you will do a chain multiplication many times with matrices of the same size. Another use case for this is to do the analysis on the CPU but use the generated function on the GPU. In either case, you only need to find the optimal order and JIT compile that function once. To save this function to the global scope:

```julia
a = rand(Float32, 10000,2000)
b = rand(Float32, 2000,100)
c = rand(Float32, 100)

eval(matrixchainmultiply("mcm_abc", a,b,c))

mcm_abc(a,b,c)

# multiple-dispatch is nice; taking advantage below
using ArrayFire
aa = AFArray(a)
ab = AFArray(b)
ac = AFArray(c)

mcm_abc(aa,ab,ac)  # no new analysis; ke
```
"""
function matrixchainmultiply(fn_name::String, A...)
  matrixchainmultiply_fn(fn_name, A...)
end
function matrixchainmultiply_fn(fn_name::String, A...)
  p = mcm_makep(A...)
  m, s, n = mcm_cost_matrices(p)
  ex = parse(fn_name * "(A...) = " * mcm_print(s, 1, n))
  cost = m[1,n]
  println("Operation: $ex")
  println("Cost: $cost")
  ex
end

"Cost information tells which order to do."
function mcm_cost_matrices(p::Vector{Int})
  ∞ = typemax(Int)  # watch for overflow?
  n = length(p) - 1  # number of matrices
  @assert (n > 1) "TODO: allow to use for short lists"
  # stores costs
  m = zeros(Int, (n, n))
  # records corresponding indices
  s = zeros(Int, (n-1, n-1))  # Cormen has 2nd index going 2..n
  # zero on diagonals
  for i ∈ 1:n
    m[i,i] = 0
  end
  # l is chain length
  for l ∈ 2:n
    for i ∈ 1:(n-l+1)
      j = i + l - 1
      m[i,j] = ∞  # this should be replaced?!
      for k ∈ i:(j-1)
        q1 = m[i,k]
        q2 = m[k+1,j]
        q3 = p[i] * p[k+1] * p[j+1]  # Cormen p vec indexed from 0
        q = +(q1, q2, q3)
        if q < m[i,j]
          m[i,j] = q
          s[i,j-1] = k
        end
      end
    end
  end
  m, s, n
end

"Use this if you want to see what the operation looks like."
function mcm_print(s::Matrix, i::Int, j::Int)
  str = ""
  if i!=j
    a = mcm_print(s, i, s[i, j-1])
    b = mcm_print(s, s[i, j-1] + 1, j)
    return "($a * $b)"
  else
    return "A[$i]"
  end
end

"""
The function describing the optimal order for the input of a specific size.
"""
function mcm_optimalorder(s, i, j, mats...)
  if i!=j
    m1 = mcm_compute(s, i, s[i, j-1], mats...)
    m2 = mcm_compute(s, s[i, j-1] + 1, j, mats...)
    return :(m1 * m2)
  else
    return mats[i]
  end
end

"Do the actual matrix chain multiplication computation."
function mcm_compute(s, i, j, mats...)
  if i!=j
    m1 = mcm_compute(s, i, s[i, j-1], mats...)
    m2 = mcm_compute(s, s[i, j-1] + 1, j, mats...)
    return (m1 * m2)
  else
    return mats[i]
  end
end

"The matrix sizes."
function mcm_makep(mats...)
  # Made some assumptions below that disallow these
  lenmats = length(mats)
  # convert to matrix-style sizes
  msizes = map(msize, mats)
  # produce unique number list
  p = zeros(Int, lenmats + 1)
  for i in 1:lenmats
    firstiter::Bool = (i==1)
    if firstiter
      p[1] = msizes[1][1]
    else
      check1 = msizes[i-1][2]
      check2 = msizes[i][1]
      @assert (check1 == check2) "Matrix sizes don't match ($check1 != $check2)!"
    end
    p[i+1] = msizes[i][2]
  end
  p
end

###

"""
Matrix-style sizes. Return type looks like `(i,j)` where `i` and `j` are `Ints`.

If other packages with matrix types (e.g., OpenCl or ArrayFire) inherit from
AbstractArray, there shouldn't be issues with other packages using this one.
"""
msize(m::AbstractMatrix) = size(m)
msize(v::AbstractVector) = (size(v,1), 1)
msize(s::Number) = (1, 1)

end # module


text:
cd(@__DIR__)
using Pkg
Pkg.activate(".")

using StochasticPrograms
using GLPK
using Pipe
using Random
using Plots
using JuMP
using DataFrames
using CSV

##
function plot_results_bb(sp, pv, w, d; scen = nothing)
    od = optimal_decision(sp)
    t = Int((length(od)-3)/2+3)
    gci = od[4:t]
    gco = od[t+1:end]
    gc = [gci[i]==0 ? -gco[i] : gci[i] for i in 1:length(gci)]
    plt = plot(gc, label = "grid connection")
    if od[1] != 0
        plot!(plt, pv.*od[1], label = "pv")
    end
    if od[2] != 0
        plot!(plt, w.*od[2], label = "wind")
    end
    plot!(plt, -d, label = "demand", xlabel = "t (h)", ylabel = "P (kW)")
    if !isnothing(scen)
    println("Objective value in scenario $scen: $(objective_value(sp, scen))")
    #t = Int((length(od)-3)/2+3)
        ord = optimal_recourse_decision(sp, scen)
        stor_flow = ord[1:t]
        stor = [sum(stor_flow[1:i]) for i in 1:t]
        stor.+=0.5*od[3]
        buyback = ord[t:end]
        plot!(plt, buyback, label = "buy-back")
        plt_st = plot(stor, label = "storage charge")
        plt_general = plot(plt, plt_st, layout = (2, 1))
        display(plt_general)
    else display(plt)
    end
end
##
# Global system parameters
c_i = .03
c_o = .01
timesteps = 1:48
time_scale = 10. *365*24/length(timesteps)
c_pv = 1000.
c_wind = 1000.
c_storage = 400.
c_flex = .5
F = 100.
storage_scale = 1.
##
# Define weather and demand data
pv = CSV.read("../data/pv_Halle18.csv", DataFrame)[timesteps, 1]
wind = CSV.read("../data/wind_Karholz.csv", DataFrame)[timesteps, 1]
demand = CSV.read("../data/demand_Industriepark.csv", DataFrame)[timesteps, 1]
#=pv = CSV.read("../data/basic_example.csv", DataFrame)[timesteps, 3]
wind = CSV.read("../data/basic_example.csv", DataFrame)[timesteps, 4]
demand = CSV.read("../data/basic_example.csv", DataFrame)[timesteps, 2]
=#
##
function scenarios(n, timesteps)
    return [@scenario t_xi = rand(1:length(timesteps)-2) s_xi = rand([-1]) probability = 1. /n for i in 1:n]
end
##
#= Note on signs: s_xi<0 means that energy is requested, s_xi>0 means an additional consumption
To simplify the problem and avoid using inequalities in constraints, let's first study only s_xi<0
Then buy-back B>0, therefore it's associated with price c_i
=#
@stochastic_model em begin
    @stage 1 begin
        @decision(em, u_pv>=0)
        @decision(em, u_wind>=0)
        @decision(em, u_storage>=0)
        @decision(em, gci[t in timesteps]>=0)
        @decision(em, gco[t in timesteps]>=0)
        @objective(em, Min, u_pv*c_pv/time_scale+u_wind*c_wind/time_scale+u_storage*c_storage/time_scale+c_i*sum(gci)-c_o*sum(gco))
    end
    @stage 2 begin
        @uncertain t_xi s_xi #t_xi the time of flexibility demand, s_xi - sign (±1 or 0)
        @known(em, u_pv)
        @known(em, u_wind)
        @known(em, u_storage)
        @known(em, gci)
        @known(em, gco)
        @recourse(em, storage[t in timesteps])
        @recourse(em, B[t in timesteps])
        @constraint(em, [t in timesteps], s_xi*B[t]<=0)
        @constraint(em, [t in 1:t_xi], B[t]==0)
        @constraint(em, [t in 1:(t_xi-1)],
        gci[t]-gco[t]+u_pv*pv[t]+u_wind*wind[t]-demand[t]+storage[t]==0)
        @constraint(em,
        gci[t_xi]-gco[t_xi]+u_pv*pv[t_xi]+u_wind*wind[t_xi]-demand[t_xi]+storage[t_xi]+F*s_xi==0)
        @constraint(em, [t in (t_xi+1):length(timesteps)],
        gci[t]-gco[t]+u_pv*pv[t]+u_wind*wind[t]-demand[t]+storage[t]+B[t]==0)
        @objective(em, Min, c_i*sum(B))
        @constraint(em, [t in timesteps], -0.5*u_storage*storage_scale <= sum(storage[1:t]))
        @constraint(em, [t in timesteps], sum(storage[1:t]) <= 0.5*u_storage*storage_scale)
        @constraint(em, sum(storage) == 0)
        @constraint(em, sum(B) == -s_xi* F)
    end
end
##
#xi = scenarios(10, timesteps)
xi_1 = @scenario t_xi = 2 s_xi = -1 probability = 0.5
xi_2 = @scenario t_xi = 1 s_xi = -1 probability = 0.5
xi = scenarios(10, timesteps)
sp = instantiate(em, xi, optimizer = GLPK.Optimizer)
##
optimize!(sp)
objective_value(sp)

od = optimal_decision(sp)
##
plot_results_bb(sp, pv, wind, demand)

##
# Main result
println("Termination status: $(termination_status(sp))")
println("Objective value: $(objective_value(sp))")
println("Optimal decision: $(optimal_decision(sp))")

##
u_pv = sp[1, :u_pv]
u_wind = sp[1, :u_wind]
u_storage = sp[1, :u_storage]
# First stage
println("value(u_pv) = $(value(u_pv))")
println("value(u_wind) = $(value(u_wind))")
println("value(u_storage) = $(value(u_storage))")

# Second stage
for s in 1:3
    #println("Objective value in scenario $s: $(objective_value(sp, s))")
    println("Optimal recourse in scenario $s: $(optimal_recourse_decision(sp, s))")
    plot_results_bb(sp, pv, wind, demand, scen = s)
end

##

plot_results_bb(sp, pv, wind, demand, scen = 1)

max_stars_repo_path: src/prototype buyback.jl

max_issues_repo_path: src/prototype buyback.jl

max_forks_repo_path: src/prototype buyback.jl


text:
struct Heisenberg{D} <: AbstractHeisenberg{D}
    size::NTuple{D, Int}
    periodic::Bool
    Heisenberg(size::Int...; periodic::Bool) = new{length(size)}(size, periodic)
end

Base.size(model::Heisenberg) = model.size

heisenberg_ij(nbit::Int, i::Int, j::Int=i+1) = put(nbit, i=>X)*put(nbit, j=>X) + put(nbit, i=>Y)*put(nbit, j=>Y) + put(nbit, i=>Z)*put(nbit, j=>Z)
const heisenberg_term = repeat(2, X, 1:2) + repeat(2, Y, 1:2) + repeat(2, Z, 1:2)

function hamiltonian(model::Heisenberg)
    nbit = nspin(model)
    sum(x->heisenberg_ij(nbit, x[1], x[2]), get_bonds(model))*0.25
end

function get_bonds(model::Heisenberg{2})
    m, n = model.size
    cis = LinearIndices(model.size)
    bonds = Tuple{Int, Int, Float64}[]
    for i=1:m, j=1:n
        (i!=m || model.periodic) && push!(bonds, (cis[i,j], cis[i%m+1,j], 1.0))
        (j!=n || model.periodic) && push!(bonds, (cis[i,j], cis[i,j%n+1], 1.0))
    end
    bonds
end

function get_bonds(model::AbstractModel{1})
    nbit, = model.size
    [(i, i%nbit+1, 1.0) for i in 1:(model.periodic ? nbit : nbit-1)]
end

"""
    energy(chem::QuantumMPS, model::AbstractHeisenberg) -> Float64

Ground state energy by sampling Quantum MPS.
The hamiltonian is limited to Heisenberg and J1J2 Type.
"""
function energy(chem::QuantumMPS, model::AbstractHeisenberg)
    energy(chem, X, model) + energy(chem, Y, model) + energy(chem, Z, model)
end

function energy(chem::QuantumMPS, pauli::PauliGate, model::AbstractHeisenberg)
    res = gensample(chem, pauli)
    local eng = 0.0
    for bond in get_bonds(model)
        eng += bond[3]*mean(res[:,bond[1]].*res[:,bond[2]])
    end
    eng/4
end

max_stars_repo_name: GiggleLiu/QuantumMPS

max_issues_repo_name: GiggleLiu/QuantumMPS

max_forks_repo_name: GiggleLiu/QuantumMPS


text:
module Test12

import Automa
import Automa.RegExp: @re_str
const re = Automa.RegExp
using Test

@testset "Test12" begin
    a = re"a*"
    a.actions[:all] = [:a]
    machine = Automa.compile(a)

    ctx = Automa.CodeGenContext()
    @eval function validate(data)
        logger = Symbol[]
        $(Automa.generate_init_code(ctx, machine))
        p_end = p_eof = sizeof(data)
        $(Automa.generate_exec_code(ctx, machine, :debug))
        return logger, cs == 0 ? :ok : cs < 0 ? :error : :incomplete
    end
    @test validate(b"") == ([], :ok)
    @test validate(b"a") == ([:a], :ok)
    @test validate(b"aa") == ([:a, :a], :ok)
    @test validate(b"aaa") == ([:a, :a, :a], :ok)
    @test validate(b"aaab") == ([:a, :a, :a], :error)
end

end


text:
struct Togglers
    options::Vector{Option}
end

function jsrender(session::Session, togglers::Togglers)
    options = togglers.options
    toggles = map(options) do entry
        selected = entry.selected
        isoriginal = entry.value.isoriginal
        reset = Observable(true)
        JSServe.register_resource!(session, reset)
        on(session, reset) do _
            reset!(entry.value)
        end
        modified = DOM.span(
            class="float-right p-4 inline-block hover:text-red-300",
            "⬤",
            style=isoriginal[] ? "display:none" : "display:inline",
            onclick=string(js"JSServe.update_obs($(reset), true)")
        )
        onjs(session, isoriginal, js"""
            function (val) {
                $(modified).style.display = val ? "none" : "inline";
            }
        """
        )
        content = DOM.div(style="display:none", class="p-4 bg-white rounded-b", jsrender(session, entry.value))
        button = DOM.button(
            class="text-blue-800 text-xl font-semibold border-b-2 border-gray-200 hover:bg-gray-200 w-full text-left",
            onclick=string(js"""
                if (!$(modified).isEqualNode(event.target)) {
                    JSServe.update_obs($selected, !(JSServe.get_observable($selected)))
                }
            """), # FIXME: report to JSServe that this requires string to work on reload
            DOM.span(class="pl-4 py-4 inline-block", entry.key),
            modified
        )
        onjs(
            session,
            selected,
            js"""
                function (val) {
                    $(content).style.display = val ? "block" : "none";
                    $(button).classList.toggle("border-b-2")
                }
            """
        )
        return DOM.div(button, content)
    end
    return DOM.div(toggles...)
end


text:
struct HalfEdge
    twin::UInt32
    next::UInt32
end

twin(he::HalfEdge) = he.twin
next(he::HalfEdge) = he.next

struct SimplicialSurface{P <: AbstractPoint}
    vertices::Array{P, 1}
    #halfedges::Array{HalfEdge, 1} # not implemented yet
    edges::Array{SVector{2, UInt32}, 1}
    faces::Array{SVector{3, UInt32}, 1}
    edgetoindex::Dict{SVector{2, UInt32}, UInt32}
    facetoindex::Dict{SVector{3, UInt32}, UInt32}
    edgevertexadjacency::SparseMatrixCSC{UInt32, UInt32}
    faceedgeadjacency::SparseMatrixCSC{UInt32, UInt32}
end

vertices(simplicsurf::SimplicialSurface{<: AbstractPoint}) = simplicsurf.vertices
#halfedges(simplicsurf::SimplicialSurface{<: AbstractPoint}) = simplicsurf.halfedges
edges(simplicsurf::SimplicialSurface{<: AbstractPoint}) = simplicsurf.edges
faces(simplicsurf::SimplicialSurface{<: AbstractPoint}) = simplicsurf.faces
edgetoindex(simplicsurf::SimplicialSurface{<: AbstractPoint}) = simplicsurf.edgetoindex
facetoindex(simplicsurf::SimplicialSurface{<: AbstractPoint}) = simplicsurf.facetoindex
edgevertexadjacency(simplicsurf::SimplicialSurface{<: AbstractPoint}) = simplicsurf.edgevertexadjacency
faceedgeadjacency(simplicsurf::SimplicialSurface{<: AbstractPoint}) = simplicsurf.faceedgeadjacency

sv(x::Vararg{Integer}) = SVector{length(x), UInt32}(x...)
sv(x::NTuple{N, Integer}) where N = sv(x...)

sortsv(x::Vararg{Integer}) = sort(sv(x...))
sortsv(x::NTuple{N, Integer}) where N = sortsv(x...)

struct SurfaceSubset
    vertices::Array{UInt32, 1}
    edges::Array{UInt32, 1}
    faces::Array{UInt32, 1}
end

vertices(subset::SurfaceSubset) = subset.vertices
edges(subset::SurfaceSubset) = subset.edges
faces(subset::SurfaceSubset) = subset.faces

function buildsubsetvector(subset::SurfaceSubset, surface::SimplicialSurface, subsetfunc)::Array{UInt32, 1}
    cursubset = subsetfunc(subset)
    numsubsettotal = length(subsetfunc(surface))
    subsetvector = zeros(UInt32, numsubsettotal)
    for i in 1:length(cursubset)
        subsetvector[cursubset[i]] = one(UInt32)
    end
    return subsetvector
end

buildvertexvector(subset::SurfaceSubset, surface::SimplicialSurface) = buildsubsetvector(subset, surface, vertices)
buildedgevector(subset::SurfaceSubset, surface::SimplicialSurface) = buildsubsetvector(subset, surface, edges)
buildfacevector(subset::SurfaceSubset, surface::SimplicialSurface) = buildsubsetvector(subset, surface, faces)
buildallsubsetvectors(subset::SurfaceSubset, surface::SimplicialSurface) = map(x->buildsubsetvector(subset, surface, x), (vertices, edges, faces))

function star(subset::SurfaceSubset, surface::SimplicialSurface)::SurfaceSubset
    originalvertices, originaledges, originalfaces = buildallsubsetvectors(subset, surface)
    starvertices = originalvertices
    staredges = subsetnormalize(originaledges + edgevertexadjacency(surface) * starvertices)
    starfaces = subsetnormalize(originalfaces + faceedgeadjacency(surface) * staredges)
    starsubset = SurfaceSubset(map(tononzeroindices, (starvertices, staredges, starfaces))...)
    return starsubset
end

function closure(subset::SurfaceSubset, surface::SimplicialSurface)::SurfaceSubset
    originalvertices, originaledges, originalfaces = buildallsubsetvectors(subset, surface)
    closurefaces = originalfaces
    closureedges = subsetnormalize(originaledges + faceedgeadjacency(surface)' * closurefaces)
    closurevertices = subsetnormalize(originalvertices + edgevertexadjacency(surface)' * closureedges)
    closuresubset = SurfaceSubset(map(tononzeroindices, (closurevertices, closureedges, closurefaces))...)
    return closuresubset
end

function link(subset::SurfaceSubset, surface::SimplicialSurface)::SurfaceSubset
    starofclosure = star(closure(subset, surface), surface)
    closureofstar = closure(star(subset, surface), surface)
    return subsetdifference(closureofstar, starofclosure, surface)
end

function subsetdifference(subset1::SurfaceSubset, subset2::SurfaceSubset, surface::SimplicialSurface)
    subsetvectors1 = buildallsubsetvectors(subset1, surface)
    subsetvectors2 = buildallsubsetvectors(subset2, surface)

    return SurfaceSubset(map(tononzeroindices ∘ subsetdifferencehelper, zip(subsetvectors1, subsetvectors2))...)
end

subsetnormalize(x::Array{UInt32, 1})::Array{UInt32, 1} = map(a->min(a, one(UInt32)), x)

function subsetdifferencehelper(arrs::Tuple{Array{UInt32, 1}, Array{UInt32, 1}})::Array{UInt32, 1}
    (x, y) = arrs
    lambdacompare((a, b)) = UInt32((a == one(UInt32)) && (b == zero(UInt32)))
    return map(lambdacompare, zip(x, y))
end

function tononzeroindices(x::Array{UInt32, 1})::Array{UInt32, 1}
    nonzeroindices = Array{UInt32, 1}(undef, sum(x))
    curindex = 1
    for i in 1:length(x)
        if x[i] == one(UInt32)
            nonzeroindices[curindex] = i
            curindex += 1
        end
    end
    return nonzeroindices
end

function constructsimplicialsurface(vertices::AbstractArray{P, 1}, tris::AbstractArray{NgonFace{3, T}})::SimplicialSurface{P} where {P <: AbstractPoint, T}
    normedtris = map(NgonFace{3, UInt32}, tris)
    numtris = length(normedtris)

    #@show P
    faces = Array{SVector{3, UInt32}, 1}(undef, numtris)
    facetoindex = Dict{SVector{3, UInt32}, UInt32}()
    edges = Array{SVector{2, UInt32}, 1}(undef, 0)
    edgetoindex = Dict{SVector{2, UInt32}, UInt32}()
    numedges::UInt32 = zero(UInt32)
    edgevertexadjacencyCOO = Array{Tuple{UInt32, UInt32, UInt32}, 1}(undef, 0)
    faceedgeadjacencyCOO = Array{Tuple{UInt32, UInt32, UInt32}, 1}(undef, 0)
    #halfedgenext = Dict{UInt32, UInt32}()

    indiceshelper = ((1, 2), (2, 3), (1, 3))

    # Iterate over all faces, constructing indices and adjacency data as we go
    for i in 1:numtris
        # Construct face data
        triindcart = sortsv(normedtris[i]...)
        faces[i] = triindcart
        facetoindex[triindcart] = i

        # Iterate over edges in this face
        for j in 1:3
            # Construct edge data
            ind1 = triindcart[indiceshelper[j][1]]
            ind2 = triindcart[indiceshelper[j][2]]
            edgeindcart = sv(ind1, ind2) # don't need sortsv since already sorted from above
            #@show edgeindcart

            # New edge
            if !haskey(edgetoindex, edgeindcart)
                numedges += one(UInt32)
                #@show numedges
                push!(edges, edgeindcart)
                edgetoindex[edgeindcart] = numedges

                # Construct edge vertex adjacency matrix data
                push!(edgevertexadjacencyCOO, (numedges, ind1, one(UInt32)))
                push!(edgevertexadjacencyCOO, (numedges, ind2, one(UInt32)))
            end

            # Construct face edge adjacency matrix data
            edgeindlin = edgetoindex[edgeindcart]
            push!(faceedgeadjacencyCOO, (i, edgeindlin, one(UInt32)))
        end
    end

    #@show edgevertexadjacencyCOO
    #@show faceedgeadjacencyCOO

    edgevertexadjacency = sparse(unzip(edgevertexadjacencyCOO)...)
    #@show edgevertexadjacency
    faceedgeadjacency = sparse(unzip(faceedgeadjacencyCOO)...)
    #@show faceedgeadjacency

    #halfedges = Array{HalfEdge, 1}(undef, numedges * 2)
    #halfedges = Array{HalfEdge, 1}(undef, 0)

    #simplicsurf = SimplicialSurface{P}(vertices, halfedges, edges, faces, edgetoindex, facetoindex, edgevertexadjacency, faceedgeadjacency)
    simplicsurf = SimplicialSurface{P}(vertices, edges, faces, edgetoindex, facetoindex, edgevertexadjacency, faceedgeadjacency)

    return simplicsurf
end

function constructsimplicialsurface(mesh::AbstractMesh)
    return constructsimplicialsurface(coordinates(mesh), faces(mesh))
end

function initfortest()
    mesh = load("meshes/original/tri.ply")
    return mesh
end

unzip(a) = map(x->getfield.(a, x), fieldnames(eltype(a)))

function buildlaplacianoperator(surface::SimplicialSurface{P})::Tuple{SparseMatrixCSC{Float64, Int64}, Array{Float64, 1}} where {P <: AbstractPoint}
    tris = faces(surface)
    numtris = length(tris)
    verts = vertices(surface)
    numverts = length(verts)

    dualvertareas = zeros(Float64, numverts)
    laplacianoperatordict = Dict{Tuple{Int64, Int64}, Float64}()
    for i in 1:numverts
        laplacianoperatordict[i, i] = Float64(1e-8)
    end

    vertexpermutations = ((1, 2, 3), (1, 3, 2), (2, 1, 3), (2, 3, 1), (3, 1, 2), (3, 2, 1))
    for i in 1:numtris
        trivertinds = tris[i]
        #@show trivertinds
        triverts = map(ind->verts[ind], trivertinds)
        #@show triverts, typeof(triverts)
        #@show triverts[1], typeof(triverts[1])

        triarea = norm(cross(triverts[2] - triverts[1], triverts[3] - triverts[1])) / 2
        #@show triarea
        for j in 1:3
            dualvertareas[trivertinds[j]] += Float64(triarea / 3)
        end
        for j in 1:6
            indicesatplay = vertexpermutations[j]
            primaryvertind = trivertinds[indicesatplay[1]]
            secondaryvertind = trivertinds[indicesatplay[2]]
            tertiaryvertind = trivertinds[indicesatplay[3]]
            #@show primaryvertind, secondaryvertind, tertiaryvertind
            primaryvert = triverts[indicesatplay[1]]
            secondaryvert = triverts[indicesatplay[2]]
            tertiaryvert = triverts[indicesatplay[3]]
            #@show primaryvert, secondaryvert, tertiaryvert

            cotalpha = cotan(primaryvert - tertiaryvert, secondaryvert - tertiaryvert) / 2
            #@show cotalpha

            if !haskey(laplacianoperatordict, (primaryvertind, secondaryvertind))
                laplacianoperatordict[(Int64(primaryvertind), Int64(secondaryvertind))] = -cotalpha
            else
                laplacianoperatordict[(Int64(primaryvertind), Int64(secondaryvertind))] -= cotalpha
            end
            laplacianoperatordict[(Int64(primaryvertind), Int64(primaryvertind))] += cotalpha
        end
    end

    #@show dualvertareas
    #@show laplacianoperatordict
    laplacianoperator = sparse(unzip([(row, col, Float64(value)) for ((row, col), value) in laplacianoperatordict])...)
    #@show laplacianoperator

    return laplacianoperator, dualvertareas
end

function cotan(v1::AbstractArray{Float32, 1}, v2::AbstractArray{Float32, 1})::Float32
    cos = sum(v1 .* v2)
    sin = norm(cross(v1, v2))
    return cos / (sin + 1e-10)
end

function solvepoissonproblem(surface::SimplicialSurface{P}, potential::Array{Float64, 1}) where {P <: AbstractPoint}
    lap, dualareas = buildlaplacianoperator(surface)
    rhs = potential .* dualareas
    return lap \ rhs
end

function normalizetozeroone(x::Array)
    fullposx = (x .+ minimum(x))
    return fullposx ./ maximum(fullposx)
end

solvenormpoisson(surface, potential) = normalizetozeroone(solvepoissonproblem(surface, potential))

function solveheatstep(surface::SimplicialSurface{P}, potential::Array{Float64, 1}, stepsize::Float64) where {P <: AbstractPoint}
    lap, dualareas = buildlaplacianoperator(surface)
    numverts = length(dualareas)
    lhs = spdiagm(numverts, numverts, 0 => dualareas) + stepsize .* lap
    rhs = potential .* dualareas
    return lhs \ rhs
end

function solvecurvatureflowstep(surface::SimplicialSurface{P}, stepsize::Float64) where {P <: AbstractPoint}
    lap, dualareas = buildlaplacianoperator(surface)
    numverts = length(dualareas)
    lhs = spdiagm(numverts, numverts, 0 => dualareas) + stepsize .* lap
    rhs = dualareas .* convertpointarraytomatrix(vertices(surface))
    return lhs \ rhs
end

function solvewavestep(surface::SimplicialSurface{P}, value::Array{Float64}, valuevelocity::Array{Float64}, forcing::Array{Float64}, frequency::Float64, stepsize::Float64) where {P <: AbstractPoint}
    lap, dualareas = buildlaplacianoperator(surface)
    numverts = length(dualareas)
    lhsinterior =  frequency .* stepsize .* dualareas .* lap
    lhs = [I -stepsize .* sparse(I, numverts, numverts)
           lhsinterior I]
    rhs = [value .+ forcing
           valuevelocity]

    solution = lhs \ rhs
    valuenext = solution[1:numverts]
    velocitynext = solution[numverts + 1 : 2 * numverts]
    return  valuenext, velocitynext

end

function convertpointarraytomatrix(pointarray::Array{Point3{Float32}, 1})
    numverts = length(pointarray)
    retmat = Array{Float64, 2}(undef, numverts, 3)
    for i in 1:numverts
        retmat[i, :] = Float64.(pointarray[i][:])
    end
    return retmat
end

function convertmatrixtopointarray(matrix::Array{Float64, 2})
    numverts = size(matrix)[1]
    retarray = Array{Point3{Float32}, 1}(undef, numverts)
    for i in 1:numverts
        retarray[i] = Point3{Float32}(Float32.(matrix[i, :]))
    end
    return retarray
end

function updatemeshandsurface!(mesh::AbstractMesh, surface::SimplicialSurface{P}, newverts::Array{Float64, 2}) where {P <: AbstractPoint}
    newvertpoints = convertmatrixtopointarray(newverts)
    meshcoords = coordinates(mesh)
    meshcoords .= newvertpoints
    surfacecoords = vertices(surface)
    surfacecoords .= newvertpoints
end

max_stars_repo_name: zoemcc/geometry

max_issues_repo_name: zoemcc/geometry

max_forks_repo_name: zoemcc/geometry


text:
# Activate project
using Pkg: Pkg
Pkg.activate(".")
Pkg.instantiate()

# Load required packages (sorted alphabetically)
using CSV
using DataFrames
using Dates
using DelimitedFiles
using Distributed
using Formatting
using JLD2
using Plots
using Plots.PlotMeasures
using ProgressMeter
using Random
using RCall
using Shapefile
using SimpleSDMLayers
using Statistics
using StatsBase
using StatsPlots
using ZipFile

# Load custom functions
include(joinpath("lib", "analysis.jl"))
include(joinpath("lib", "betadiv.jl"))
include(joinpath("lib", "csvdata.jl"))
include(joinpath("lib", "landcover.jl"))
include(joinpath("lib", "plotting.jl"))
include(joinpath("lib", "presence-absence.jl"))
include(joinpath("lib", "shapefiles.jl"))

# Manual version-control for oversized files
include(joinpath("lib", "version-control.jl"))
verify_raw_files()
verify_proc_files()
# verify_raw_files(; touch_placeholders=true)
# verify_proc_files(; touch_placeholders=true)

max_stars_repo_name: gabrieldansereau/betadiversity-hotspots

max_issues_repo_name: gabrieldansereau/betadiversity-hotspots

max_forks_repo_name: gabrieldansereau/betadiversity-hotspots


text:
module Reversi

using Parameters

struct State
    p1_placed::UInt64
    p2_placed::UInt64
    player::Int
    round::Int
end
State() = State(0, 0, 1, 0)

encode(r::Integer, c::Integer) = UInt64(1) << (8*(r-1) + (c-1))
encode(rc::Tuple{Integer,Integer}) = encode(rc[1], rc[2])
const MARKERS = Dict(0 => "   ", 2 => " \u25cf ", 1 => " \u25cb ")
const CENTERS_BITS = sum(encode(i) for i in [(4, 4), (5, 5), (5, 4), (4, 5)])

function show_board(io::IO, s::State)
    @unpack p1_placed, p2_placed = s
    # corners
    top_left = "\u250C"
    top_right = "\u2510"
    bottom_left = "\u2514"
    bottom_right = "\u2518"

    # lines
    hori = "\u2500"
    vert = "\u2502"

    # line + center piece
    hori_down = "\u252C"
    hori_up = "\u2534"

    print(io, " ", join(map(i->string("  ", i, " "), 1:8)), "\n")
    print(io, " ", top_left, hori, repeat(string(hori^2, hori_down, hori), 7), hori^2, top_right)

    for r in UInt64(1):8
        println(io)
        print(io, r)
        print(io, vert)
        for c in UInt64(1):8
            v = encode(r, c)
            if v & p1_placed > 0
                print(io, MARKERS[1])
            elseif v & p2_placed > 0
                print(io, MARKERS[2])
            else
                print(io, MARKERS[0])
            end
            print(io, vert)
        end
    end

    println(io, "\n ", bottom_left, hori, repeat(string(hori^2, hori_up, hori), 7), hori^2, bottom_right)
end

function Base.show(io::IO, s::State)
    @unpack player, round = s
    if player == 1 || player == 2
        msg = "It is Player $(player)'s ($(MARKERS[player])) turn"
    else
        msg = "It is player $(player)'s turn"
    end

    println(io, msg, " on round $round")
    println(io, "The score is $(score(s)) and the board looks like:\n")
    show_board(io, s)
end

# -------------- #
# bitboard stuff #
# -------------- #

function _one_directon_right_shift(g::UInt64, p::UInt64, n::Integer)
    g |= p & (g >> (n * 1))
    p &= (p >> (n * 1))
    g |= p & (g >> (n * 2))
    p &= (p >> (n * 2))
    g |= p & (g >> (n * 4))
    g
end

function _one_directon_left_shift(g::UInt64, p::UInt64, n::Integer)
    g |= p & (g << (n * 1))
    p &= (p << (n * 1))
    g |= p & (g << (n * 2))
    p &= (p << (n * 2))
    g |= p & (g << (n * 4))
    g
end

function _get_g_directions(p1::UInt64, p2::UInt64)
    mask_a = 0xfefefefefefefefe
    mask_h = 0x7f7f7f7f7f7f7f7f

    gN = _one_directon_right_shift(p1, p2, 8)
    gS = _one_directon_left_shift(p1, p2, 8)
    gE = _one_directon_left_shift(p1, p2 & mask_a, 1)
    gW = _one_directon_right_shift(p1, p2 & mask_h, 1)
    gNE = _one_directon_right_shift(p1, p2 & mask_a, 7)
    gNW = _one_directon_right_shift(p1, p2 & mask_h, 9)
    gSE = _one_directon_left_shift(p1, p2 & mask_a, 9)
    gSW = _one_directon_left_shift(p1, p2 & mask_h, 7)

    return gN, gS, gE, gW, gNE, gNW, gSE, gSW
end

function legal_actions_bits(mover::UInt64, opp::UInt64)
    occupied = mover | opp
    if (CENTERS_BITS & occupied) != CENTERS_BITS
        return CENTERS_BITS - occupied
    end

    mask_a = 0xfefefefefefefefe
    mask_h = 0x7f7f7f7f7f7f7f7f
    empty = 0xffffffffffffffff ⊻ occupied

    gN, gS, gE, gW, gNE, gNW, gSE, gSW = _get_g_directions(mover, opp)

    legal = UInt64(0)
    legal |= ((gN & ~mover) >> 8) & empty
    legal |= ((gS & ~mover) << 8) & empty
    legal |= ((gE & ~mover & mask_h) << 1) & empty
    legal |= ((gW & ~mover & mask_a) >> 1) & empty
    legal |= ((gNE & ~mover & mask_h) >> 7) & empty
    legal |= ((gNW & ~mover & mask_a) >> 9) & empty
    legal |= ((gSE & ~mover & mask_h) << 9) & empty
    legal |= ((gSW & ~mover & mask_a) << 7) & empty
    legal
end

function next_state_bits(action::UInt64, mover::UInt64, opponent::UInt64)
    flips = UInt64(0)
    gN, gS, gE, gW, gNE, gNW, gSE, gSW = _get_g_directions(
        action, opponent
    )

    mask_a = 0xfefefefefefefefe
    mask_h = 0x7f7f7f7f7f7f7f7f
    mover += action

    ((gN >> 8) & mover > 0)           && (flips |= gN)
    ((gS << 8) & mover > 0)           && (flips |= gS)
    ((gE << 1) & mask_a & mover > 0)  && (flips |= gE)
    ((gW >> 1) & mask_h & mover > 0)  && (flips |= gW)
    ((gNE >> 7) & mask_a & mover > 0) && (flips |= gNE)
    ((gNW >> 9) & mask_h & mover > 0) && (flips |= gNW)
    ((gSE << 9) & mask_a & mover > 0) && (flips |= gSE)
    ((gSW << 7) & mask_h & mover > 0) && (flips |= gSW)

    new_mover = mover | flips
    new_opponent = opponent & (~flips)

    new_mover, new_opponent
end

score(x::UInt64) = count_ones(x)
function bits_to_tuples(x::UInt64)
    out = Tuple{Int,Int}[]
    for r in 1:8
        for c in 1:8
            v = encode(r, c)
            if v & x > 0
                push!(out, (r, c))
            end
        end
    end
    out
end

legal_actions(p1::UInt64, p2::UInt64) = bits_to_tuples(legal_actions_bits(p1, p2))

function game_over(p1::UInt64, p2::UInt64)::Tuple{Bool,Int,Int}
    if p1 == 0
        return true, 0, score(p2)
    end

    if p2 == 0
        return true, score(p1), 0
    end

    occupied = p1 | p2
    if occupied == (UInt64(1) << (64)) - 1
        return true, score(p1), score(p2)
    end

    if isempty(legal_actions(p1, p2)) && isempty(legal_actions(p2, p1))
        return true, score(p1), score(p2)
    end

    false, 0, 0
end

# --------------------------- #
# bitboard api for State type #
# --------------------------- #

function legal_actions_bits(s::State)
    @unpack p1_placed, p2_placed, player = s
    if player == 1
        return legal_actions_bits(p1_placed, p2_placed)
    else
        return legal_actions_bits(p2_placed, p1_placed)
    end
end

score(s::State) = score(s.p1_placed), score(s.p2_placed)
legal_actions(s::State) = bits_to_tuples(legal_actions_bits(s))
function game_over(s::State)
    if s.round  < 4
        return false, 0, 0
    end
    game_over(s.p1_placed, s.p2_placed)
end

function do_action(s::State, action::UInt64)::State
    @unpack p1_placed, p2_placed, player = s
    if player == 1
        new_p1, new_p2 = next_state_bits(action, p1_placed, p2_placed)
        new_player = 2
    else
        new_p2, new_p1 = next_state_bits(action, p2_placed, p1_placed)
        new_player = 1
    end
    State(new_p1, new_p2, new_player, s.round+1)
end

do_action(s::State, a::Tuple{Int,Int}) = do_action(s, encode(a))
do_pass(s::State) = State(s.p1_placed, s.p2_placed, (s.player == 1) + 1, s.round+1)

opponent(i::Integer) = (i == 1) + 1
opponent(s::State) = opponent(state.player)

# ---------- #
# Simulation #
# ---------- #
abstract type Player end

"Using players `p1` and `p2`, simulate the game starting from state `s`"
function play_game(p1::Player, p2::Player, s::State=State(); disp::Bool=true)
    players = (p1, p2)
    history = State[]
    while !(game_over(s)[1])
        push!(history, s)
        disp && show(s)
        actions_bits = legal_actions_bits(s)
        if actions_bits == 0
            s = do_pass(s)
            continue
        end
        action = select_action(players[s.player], s)
        s = do_action(s, action)
    end
    score(s), history
end

# ------- #
# Players #
# ------- #

select_action(p::Player, s::State) = select_action(p, s, legal_actions(s))
select_action(p::Player, s::State, g::Base.Generator) = select_action(p, s, collect(g))
select_action(p::Player, s::State, a::UInt64) = select_action(p, s, bits_to_tuples(a))

# ------------ #
# RandomPlayer #
# ------------ #

struct RandomPlayer <: Player end
select_action(p::RandomPlayer, s::State, actions::AbstractVector) = rand(actions)

# ---------- #
# Evaluators #
# ---------- #
# TODO: optimize over these weights
const weights = Float64[
        120 -20  20   5   5  20 -20 120
        -20 -40  -5  -5  -5  -5 -40 -20
         20  -5  15   3   3  15  -5  20
          5  -5   3   3   3   3  -5   5
          5  -5   3   3   3   3  -5   5
         20  -5  15   3   3  15  -5  20
        -20 -40  -5  -5  -5  -5 -40 -20
        120 -20  20   5   5  20 -20 120
    ]
function weighted_evaluator(s::State)
    p1_score = 0.0
    for (r, c) in bits_to_tuples(s.p1_placed)
        p1_score += weights[r, c]
    end

    p2_score = 0.0
    for (r, c) in bits_to_tuples(s.p1_placed)
        p1_score += weights[r, c]
    end

    p1_score, p2_score
end

# ------- #
# Minimax #
# ------- #

struct Minimax{Teval} <: Player
    evaluator::Teval
    depth::Int
end

function minimax(p::Minimax, s::State, depth::Int)
    if depth == 0 || game_over(s)[1]
        vals = p.evaluator(s)
        val = vals[1] - vals[2]
        if s.player == 1
            return val, (0, 0)
        else
            return -val, (0, 0)
        end
    end

    value(a_state::State) = -minimax(p, a_state, depth-1)[1]

    actions = legal_actions(s)

    if length(actions) == 0  # this guy can't go. Let him pass
        return value(do_pass(s)), (0, 0)
    end

    return maximum((value(do_action(s, a)), a) for a in actions)
end

select_action(p::Minimax, s::State) = minimax(p, s, p.depth)[2]

# ------------------ #
# Stochastic Minimax #
# ------------------ #

struct StochasticMinimax{Teval} <: Player
    evaluator::Teval
    depth::Int
end

function minimax(p::StochasticMinimax, s::State, depth::Int)::Tuple{Float64,Tuple{Int,Int}}
    if depth == 0 || game_over(s)[1]
        vals = p.evaluator(s)
        val = vals[1] - vals[2]
        if s.player == 1
            return val, (0, 0)
        else
            return -val, (0, 0)
        end
    end

    actions = legal_actions(s)

    if length(actions) == 0  # this guy can't go. Let him pass
        return -minimax(p, do_pass(s), depth-1)[1], (0 ,0)
    end

    values = [-minimax(p, do_action(s, a), depth-1)[1] for a in actions]

    probs = values + abs(minimum(values)) + 1  # shift up so they all > 0
    probs = cumsum(probs ./ sum(probs))        # turn into cdf
    index = searchsortedfirst(probs, rand())   # sample from the cdf

    return values[index], actions[index]
end

select_action(p::StochasticMinimax, s::State) = minimax(p, s, p.depth)[2]

# ------------------ #
# Alpha beta pruning #
# ------------------ #

mutable struct AlphaBetaPruning{Teval} <: Player
    evaluator::Teval
    depth::Int
    searched::Int
end

AlphaBetaPruning(evaluator, depth::Int) = AlphaBetaPruning(evaluator, depth, 0)

function alphabeta(p::AlphaBetaPruning, s::State, depth::Int, alpha, beta)
    if depth == 0 || game_over(s)[1]
        vals = p.evaluator(s)
        val = vals[1] - vals[2]
        if s.player == 1
            return val, (0, 0)
        else
            return -val, (0, 0)
        end
    end

    actions = legal_actions(s)

    if length(actions) == 0  # this guy can't go. Let him pass
        return -alphabeta(p, do_pass(s), depth-1, -beta, -alpha)[1], (0, 0)
    end

    best_action = actions[1]
    for action in actions
        p.searched += 1
        if alpha >= beta
            break
        end
        val = -alphabeta(p, do_action(s, action), depth-1, -beta, -alpha)[1]
        if val > alpha
            alpha = val
            best_action = action
        end
    end
    return alpha, best_action
end

function select_action(p::AlphaBetaPruning, s::State)
    start_time = time()
    p.searched = 0
    out = alphabeta(p, s, p.depth, -Inf, Inf)[2]
    end_time = time()
    print("From state ($(s.p1_placed), $(s.p2_placed)) I searched $(p.searched) actions")
    println(" in $(end_time - start_time) seconds")
    out
end

# ----------------- #
# Tournament Server #
# ----------------- #
mutable struct TournamentPlayer{Tp<:Player}
    player::Tp
    me::Int
    them::Int
    sock::TCPSocket
    time1::Float64
    time2::Float64
end

function TournamentPlayer(p::Player, me::Int)
    TournamentPlayer(p, me, (me == 1) + 1, TCPSocket(), NaN, NaN)
end

function read_message(tp::TournamentPlayer)
    sock = tp.sock
    turn = parse(Int, readline(sock))
    round = parse(Int, readline(sock))
    tp.time1 = parse(Float64, readline(sock))
    tp.time2 = parse(Float64, readline(sock))

    # read state from server
    p1_placed = UInt64(0)
    p2_placed = UInt64(0)
    for i in 1:8, j in 1:8
        owned = parse(Int, readline(sock))
        if owned == 1
            p1_placed += encode(9-i, j)
        elseif owned == 2
            p2_placed += encode(9-i, j)
        end
    end
    # terminate message
    readline(sock)

    # return state
    State(p1_placed, p2_placed, turn, round)
end

play_move(tp::TournamentPlayer, rc::Tuple{Int, Int}) = play_move(tp, rc[1], rc[2])
play_move(tp::TournamentPlayer, row, col) = write(tp.sock, string(8-row, '\n', col-1, '\n'))

function update_msg(tp::TournamentPlayer, s::State)
    @unpack me, them = tp
    the_score = score(s)
    if the_score[me] > the_score[them]
        println("Winning $(the_score[me]) to $(the_score[them])")
    else
        println("Losing $(the_score[me]) to $(the_score[them])")
    end
    println("Here's what the board looks like:")
    show_board(STDOUT, s)
end

function play_game(tp::TournamentPlayer)
    @unpack me, them = tp
    tp.sock = connect(3333 + me)
    @show readline(tp.sock)

    while true
        s = read_message(tp)
        if s.player == -999
            final_score = score(s)
            if final_score[me] > final_score[them]
                msg = "Congrats, you won "
            else
                msg = "Sorry, you lost "
            end
            println(msg, "by a score of ", final_score)
            break
        end
        @show tp.time1, tp.time2
        if s.player == me
            update_msg(tp, s)
            play_move(tp, select_action(tp.player, s))
        end
    end
end

#=
usage:

include("reversi.jl")
player = Reversi.AlphaBetaPruning(Reversi.weighted_evaluator, 9)

# assuming we are player 2
tp = Reversi.TournamentPlayer(player, 2)
Reversi.play_game(tp)
=#

#=
TODO:

- Make depth depend on time remaining
- Make weights depend on tile positions
- Do rollouts in parallel
- Memoize?
- Do training to optimize weights. Strategy would be:
    - Play a few thousand games with my static weights to gen training data
        - Could use my alpha-beta with softmax instead of max for move
          selection
        - Would just need to store the History::Vector{State}
    - Train a unique set of weights for each round of the game
    - Start training at latest rounds of the game. I can do a full evaluation
      at these stages because depth isn't too high
    - Loss function could be...
        - cross entropy between softmax dist implied by current weights and ???
        - MSE between greedy policy implied by the weights and "optimal" policy
    - More info here: https://www.tjhsst.edu/~rlatimer/techlab10/Per5/FourthQuarter/ChenPaperQ4-10.pdf
- Use my alpha-beta with low depth as rollout policy for MCTS in a deepmind
  style algo
- timeouts
- Iterative deepening: do depth 2 search, then 3, then 4, and so on until time
  per move runs out
=#

end # module


text:
# LSL.jl: Julia interface for Lab Streaming Layer
# Copyright (C) 2019 Samuel Powell

# test/resolve.jl: test stream resolution functions

@testset "Resolvers" begin

  info_q0 = resolve_streams()
  if length(info_q0) > 0
    @warn "Pre-existing streams exist, tests may fail"
  end

  # Create some unique outlets
  info1 = StreamInfo(name = "TestResolver1",
                     type = "TestResolverType1",
                     channel_count = 1,
                     channel_format = Float64,
                     source_id = "TestResolverSource1")

  outlet1 = StreamOutlet(info1)

  info2 = StreamInfo(name = "TestResolver2",
                     type = "TestResolverType2",
                     channel_count = 2,
                     channel_format = Float64,
                     source_id = "TestResolverSource2")

  outlet2 = StreamOutlet(info2)

  info3 = StreamInfo(name = "TestResolver3",
                     type = "TestResolverType3",
                     channel_count = 3,
                     channel_format = Float64,
                     source_id = "TestResolverSource3")

  outlet3 = StreamOutlet(info3)

  # Add a non-unique (save for source_id)
  info4 = StreamInfo(name = "TestResolver3",
                     type = "TestResolverType3",
                     channel_count = 3,
                     channel_format = Float64,
                     source_id = "TestResolverSource4")

  outlet3 = StreamOutlet(info4)

  # Wait a bit
  sleep(1.0)

  # Resolve all streams
  info_q1 = resolve_streams(timeout = 2.0)
  @test length(info_q1) == 4
  @test length(unique([source_id(i) for i in info_q1])) == 4

  # Resolve each stream by property
  info_q2 = resolve_byprop("name", "TestResolver1", timeout = 2.0)
  info_q3 = resolve_byprop("name", "TestResolver2", timeout = 2.0)
  info_q4 = resolve_byprop("name", "TestResolver3", timeout = 2.0)
  @test length(info_q2) == 1
  @test length(info_q3) == 1
  #@test_broken length(info_q4) == 2

  info_q5 = resolve_byprop("source_id", "TestResolverSource1", timeout = 2.0)
  info_q6 = resolve_byprop("source_id", "TestResolverSource2", timeout = 2.0)
  info_q7 = resolve_byprop("source_id", "TestResolverSource3", timeout = 2.0)
  info_q8 = resolve_byprop("source_id", "TestResolverSource4", timeout = 2.0)
  @test length(info_q5) == 1
  @test length(info_q6) == 1
  #@test_broken length(info_q7) == 1
  @test length(info_q8) == 1

  info_q9  = resolve_byprop("type", "TestResolverType1", timeout = 2.0)
  info_q10 = resolve_byprop("type", "TestResolverType2", timeout = 2.0)
  info_q11 = resolve_byprop("type", "TestResolverType3", timeout = 2.0)
  @test length(info_q9)  == 1
  @test length(info_q10) == 1
  #@test_broken length(info_q11) == 2

  # Resolve each stream by predicate
  info_q2 = resolve_bypred("name='TestResolver1'", timeout = 2.0)
  info_q3 = resolve_bypred("name='TestResolver2'", timeout = 2.0)
  info_q4 = resolve_bypred("name='TestResolver3'", timeout = 2.0)
  @test length(info_q2) == 1
  @test length(info_q3) == 1
  #@test_broken length(info_q4) == 2

  info_q5 = resolve_bypred("source_id='TestResolverSource1'", timeout = 2.0)
  info_q6 = resolve_bypred("source_id='TestResolverSource2'", timeout = 2.0)
  info_q7 = resolve_bypred("source_id='TestResolverSource3'", timeout = 2.0)
  info_q8 = resolve_bypred("source_id='TestResolverSource4'", timeout = 2.0)
  @test length(info_q5) == 1
  @test length(info_q6) == 1
  #@test_broken length(info_q7) == 1
  @test length(info_q8) == 1

  info_q9  = resolve_bypred("type='TestResolverType1'", timeout = 2.0)
  info_q10 = resolve_bypred("type='TestResolverType2'", timeout = 2.0)
  info_q11 = resolve_bypred("type='TestResolverType3'", timeout = 2.0)
  @test length(info_q9)  == 1
  @test length(info_q10) == 1
  #@test_broken length(info_q11) == 2

  # Force destruction of outlet and info before continuing
  LSL._destroy(outlet1)
  LSL._destroy(info1)
  LSL._destroy(outlet2)
  LSL._destroy(info2)
  LSL._destroy(outlet3)
  LSL._destroy(info3)
  LSL._destroy(outlet4)
  LSL._destroy(info4)
  sleep(0.1)

end

# Force cleanup
GC.gc()


text:
@testset "utils" begin
    v1 = [1, 2, 7, 11]
    v2 = [1, 3, 4, 5.1]
    @test AlgebraOfGraphics.mergesorted(v1, v2) == [1, 2, 3, 4, 5.1, 7, 11]
    @test AlgebraOfGraphics.mergesorted(v2, v1) == [1, 2, 3, 4, 5.1, 7, 11]
    @test_throws ArgumentError AlgebraOfGraphics.mergesorted([2, 1], [1, 3])
    @test_throws ArgumentError AlgebraOfGraphics.mergesorted([1, 2], [10, 3])

    e1 = (-3, 11)
    e2 = (-5, 10)
    @test AlgebraOfGraphics.extend_extrema(e1, e2) == (-5, 11)

    @test AlgebraOfGraphics.midpoints(1:10) == 1.5:9.5

    s = Any[1, 2]
    t = NamedArguments((a=3, b=4))
    u = NamedArguments((b=5, c=6))
    v = Any[]
    w = AlgebraOfGraphics.concatenate_values(s, t, u)
    @test w == Any[1, 2, 3, 4, 5, 6]
end

@testset "arguments" begin
    s = Arguments([10, 20, 30])
    @test s == Any[10, 20, 30]
    t = map(keys(s), s) do k, v
        return "key $k and value $v"
    end
    @test t == ["key 1 and value 10", "key 2 and value 20", "key 3 and value 30"]

    s = NamedArguments([:a, :b, :c], [10, 20, 30])
    @test s == NamedArguments((a=10, b=20, c=30))
    t = map(keys(s), s) do k, v
        return "key $k and value $v"
    end
    @test t == NamedArguments(
        [:a, :b, :c],
        ["key a and value 10", "key b and value 20", "key c and value 30"]
    )

    s = NamedArguments([:a, :b, :c], [1, 2, 3])
    odd, even = separate(isodd, s)
    @test odd == NamedArguments((a=1, c=3))
    @test even == NamedArguments((b=2,))

    t = AlgebraOfGraphics.set(s, :d => 5, :a => 3)
    @test t == NamedArguments([:a, :b, :c, :d], [3, 2, 3, 5])

    s = NamedArguments([:a, :b, :c], [1, 2, 3])
    t = AlgebraOfGraphics.set(s)
    @test s == t
    @test keys(s) !== keys(t)

    s = NamedArguments([:a, :b, :c], [1, 2, 3])
    t = AlgebraOfGraphics.filterkeys(!in((:a, :b)), s)
    @test t == NamedArguments([:c], [3])
end

@testset "layout utils" begin
    fig = Figure()
    pos = fig[1, 1]
    subpos = pos[1, 2:3]
    ax = Axis(pos[2, 1])
    @test get_layout(fig) == get_layout(pos) == get_layout(subpos) == get_layout(ax) == fig.layout

    fig = Figure()
    pos = fig[1, 2] = GridLayout()
    @test get_layout(fig) != get_layout(pos)
    @test get_layout(pos) == pos
    @test get_layout(pos).parent == get_layout(fig)
end


text:
let
    function ls(df,x,s,xtmp,g,lsr,c,mayterminate)
        LineSearches.hagerzhang!(df,x,s,xtmp,g,lsr,c,mayterminate,
                                 0.1,0.9,
                                 Inf,5.,1e-6,0.66,
                                 2)
    end
    for optimizer in (ConjugateGradient, GradientDescent, LBFGS, BFGS, Newton, AcceleratedGradientDescent, MomentumGradientDescent)
        println("Testing $(string(optimizer))")
        prob = Optim.UnconstrainedProblems.examples["Exponential"]
        optimize(prob.f, prob.initial_x, optimizer(linesearch = ls),
                 Optim.Options(autodiff = true))
        # TODO: How can I verify that the call to optimize gives a warning?
    end
end


text:
using Documenter, Xpress

makedocs(
    format = :html,
    sitename = "Xpress",
    pages = [
        "Introduction" => "index.md"
    ]
)

deploydocs(
    repo   = "github.com/JuliaOpt/Xpress.jl.git",
    target = "build",
    osname = "linux",
    julia  = "0.6",
    deps   = nothing,
    make   = nothing
)


text:
using Aqua
using DiscriminativeCircuits
using Test

@testset "Aqua tests" begin
    Aqua.test_all(DiscriminativeCircuits,
                    ambiguities = false)
end

max_stars_repo_name: Juice-jl/DiscriminativeCircuits

max_issues_repo_name: Juice-jl/DiscriminativeCircuits

max_forks_repo_name: Juice-jl/DiscriminativeCircuits


text:
import JSON
import Mustache

mutable struct NotebookOutput
end

mutable struct MarkdownOutput
end

mutable struct NowebOutput
end

mutable struct ScriptOutput
end

const output_formats = Dict{String, Any}(
  "noweb" => NowebOutput(),
  "notebook" => NotebookOutput(),
  "markdown" => MarkdownOutput(),
  "script" => ScriptOutput()
)

"""Autodetect format for converter"""
function detect_outformat(outfile::String)
  ext = lowercase(splitext(outfile)[2])

  ext == ".jl" && return "script"
  ext == ".jmd" && return "markdown"
  ext == ".ipynb" && return "notebook"
  return "noweb"
end

"""
    convert_doc(infile::AbstractString, outfile::AbstractString; format::Union{Nothing,AbstractString} = nothing)

Convert Weave documents between different formats

- `infile`: Path of the input document
- `outfile`: Path of the output document
- `format = nothing`: Output document format (optional). It will be detected automatically from the `outfile` extension. You can also specify either of `"script"`, `"markdown"`, `"notebook"`, or `"noweb"`
"""
function convert_doc(infile::AbstractString, outfile::AbstractString; format::Union{Nothing,AbstractString} = nothing)
  doc = read_doc(infile)

  if format == nothing
    format = detect_outformat(outfile)
  end

  converted = convert_doc(doc, output_formats[format])

  open(outfile, "w") do f
    write(f, converted)
  end

  return nothing
end

"""Convert Weave document to Jupyter notebook format"""
function convert_doc(doc::WeaveDoc, format::NotebookOutput)
    nb = Dict()
    nb["nbformat"] =  4
    nb["nbformat_minor"] = 2
    metadata = Dict()
    kernelspec = Dict()
    kernelspec["language"] =  "julia"
    kernelspec["name"] =  "julia-$(VERSION.major).$(VERSION.minor)"
    kernelspec["display_name"] = "Julia $(VERSION.major).$(VERSION.minor).$(VERSION.patch)"
    metadata["kernelspec"] = kernelspec
    language_info = Dict()
    language_info["file_extension"] = ".jl"
    language_info["mimetype"] = "application/julia"
    language_info["name"]=  "julia"
    language_info["version"] = "$(VERSION.major).$(VERSION.minor).$(VERSION.patch)"
    metadata["language_info"] = language_info
    cells = []
    ex_count = 1

    # Handle header
    head_tpl = """
    {{#:title}}# {{:title}}{{/:title}}
    {{#:author}}### {{{:author}}}{{/:author}}
    {{#:date}}### {{{:date}}}{{/:date}}
    """

    if isa(doc.chunks[1], DocChunk)
        doc.chunks[1] = strip_header(doc.chunks[1])
        doc.chunks[1].content[1].content = Mustache.render(head_tpl;
              [Pair(Symbol(k), v) for (k,v) in doc.header]...) * doc.chunks[1].content[1].content
    end

    for chunk in doc.chunks

        if isa(chunk, DocChunk)
            push!(cells,
              Dict("cell_type" => "markdown",
                 "metadata" => Dict(),
              "source" => [strip(join([repr(c) for c in chunk.content], ""))])
                 )
        elseif haskey(chunk.options, :skip) && chunk.options[:skip] == "notebook"
            continue
        else
            push!(cells,
            Dict("cell_type" => "code",
               "metadata" => Dict(),
               "source" => [strip(chunk.content)],
               "execution_count" => nothing,
               "outputs" => []
               ))
       end
   end

    nb["cells"] = cells
    nb["metadata"] = metadata

    json_nb = JSON.json(nb, 2)
    return json_nb
end

"""Convert Weave document to Jupyter notebook format"""
function convert_doc(doc::WeaveDoc, format::MarkdownOutput)
  output = ""
  for chunk in doc.chunks
    if isa(chunk, DocChunk)
      output *= join([repr(c) for c in chunk.content], "")
    else
      output *= "\n" * "```julia"
      isempty(chunk.optionstring) || (output *= ";" * chunk.optionstring)
      output *= "\n" * lstrip(chunk.content)
      output *= "```\n"
    end
  end

  return output
end

"""Convert Weave document to noweb format"""
function convert_doc(doc::WeaveDoc, format::NowebOutput)
  output = ""
  for chunk in doc.chunks
    if isa(chunk, DocChunk)
      output *= join([repr(c) for c in chunk.content], "")
    else
      output *= "\n" * "<<"
      isempty(chunk.optionstring) || (output *= strip(chunk.optionstring))
      output *= ">>="
      output *= "\n" * lstrip(chunk.content)
      output *= "@\n"
    end
  end

  return output
end

"""Convert Weave document to script format"""
function convert_doc(doc::WeaveDoc, format::ScriptOutput)
  output = ""
  for chunk in doc.chunks
    if typeof(chunk) == Weave.DocChunk
      content = join([repr(c) for c in chunk.content], "")
      output *= join(["#' " * s for s in split(content, "\n")], "\n")
    else
      output *= "\n#+ "
      isempty(chunk.optionstring) || (output *= strip(chunk.optionstring))
      output *= "\n\n" * lstrip(chunk.content)
      output *= "\n"
    end
  end

  return output
end

function Base.repr(c::InlineText)
  return c.content
end

function Base.repr(c::InlineCode)
  return "`j $(c.content)`"
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct VerificationIPFlowParameters <: SwaggerModel
    targetResourceId::Any # spec type: Union{ Nothing, String } # spec name: targetResourceId
    direction::Any # spec type: Union{ Nothing, Direction } # spec name: direction
    protocol::Any # spec type: Union{ Nothing, String } # spec name: protocol
    localPort::Any # spec type: Union{ Nothing, String } # spec name: localPort
    remotePort::Any # spec type: Union{ Nothing, String } # spec name: remotePort
    localIPAddress::Any # spec type: Union{ Nothing, String } # spec name: localIPAddress
    remoteIPAddress::Any # spec type: Union{ Nothing, String } # spec name: remoteIPAddress
    targetNicResourceId::Any # spec type: Union{ Nothing, String } # spec name: targetNicResourceId

    function VerificationIPFlowParameters(;targetResourceId=nothing, direction=nothing, protocol=nothing, localPort=nothing, remotePort=nothing, localIPAddress=nothing, remoteIPAddress=nothing, targetNicResourceId=nothing)
        o = new()
        validate_property(VerificationIPFlowParameters, Symbol("targetResourceId"), targetResourceId)
        setfield!(o, Symbol("targetResourceId"), targetResourceId)
        validate_property(VerificationIPFlowParameters, Symbol("direction"), direction)
        setfield!(o, Symbol("direction"), direction)
        validate_property(VerificationIPFlowParameters, Symbol("protocol"), protocol)
        setfield!(o, Symbol("protocol"), protocol)
        validate_property(VerificationIPFlowParameters, Symbol("localPort"), localPort)
        setfield!(o, Symbol("localPort"), localPort)
        validate_property(VerificationIPFlowParameters, Symbol("remotePort"), remotePort)
        setfield!(o, Symbol("remotePort"), remotePort)
        validate_property(VerificationIPFlowParameters, Symbol("localIPAddress"), localIPAddress)
        setfield!(o, Symbol("localIPAddress"), localIPAddress)
        validate_property(VerificationIPFlowParameters, Symbol("remoteIPAddress"), remoteIPAddress)
        setfield!(o, Symbol("remoteIPAddress"), remoteIPAddress)
        validate_property(VerificationIPFlowParameters, Symbol("targetNicResourceId"), targetNicResourceId)
        setfield!(o, Symbol("targetNicResourceId"), targetNicResourceId)
        o
    end
end # type VerificationIPFlowParameters

const _property_map_VerificationIPFlowParameters = Dict{Symbol,Symbol}(Symbol("targetResourceId")=>Symbol("targetResourceId"), Symbol("direction")=>Symbol("direction"), Symbol("protocol")=>Symbol("protocol"), Symbol("localPort")=>Symbol("localPort"), Symbol("remotePort")=>Symbol("remotePort"), Symbol("localIPAddress")=>Symbol("localIPAddress"), Symbol("remoteIPAddress")=>Symbol("remoteIPAddress"), Symbol("targetNicResourceId")=>Symbol("targetNicResourceId"))
const _property_types_VerificationIPFlowParameters = Dict{Symbol,String}(Symbol("targetResourceId")=>"String", Symbol("direction")=>"Direction", Symbol("protocol")=>"String", Symbol("localPort")=>"String", Symbol("remotePort")=>"String", Symbol("localIPAddress")=>"String", Symbol("remoteIPAddress")=>"String", Symbol("targetNicResourceId")=>"String")
Base.propertynames(::Type{ VerificationIPFlowParameters }) = collect(keys(_property_map_VerificationIPFlowParameters))
Swagger.property_type(::Type{ VerificationIPFlowParameters }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_VerificationIPFlowParameters[name]))}
Swagger.field_name(::Type{ VerificationIPFlowParameters }, property_name::Symbol) =  _property_map_VerificationIPFlowParameters[property_name]

const _allowed_VerificationIPFlowParameters_protocol = ["TCP", "UDP"]

function check_required(o::VerificationIPFlowParameters)
    (getproperty(o, Symbol("targetResourceId")) === nothing) && (return false)
    (getproperty(o, Symbol("direction")) === nothing) && (return false)
    (getproperty(o, Symbol("localPort")) === nothing) && (return false)
    (getproperty(o, Symbol("remotePort")) === nothing) && (return false)
    (getproperty(o, Symbol("localIPAddress")) === nothing) && (return false)
    (getproperty(o, Symbol("remoteIPAddress")) === nothing) && (return false)
    true
end

function validate_property(::Type{ VerificationIPFlowParameters }, name::Symbol, val)
    if name === Symbol("protocol")
        Swagger.validate_param(name, "VerificationIPFlowParameters", :enum, val, _allowed_VerificationIPFlowParameters_protocol)
    end
end


text:
# dev test oplus and ominus for SE3

using TransformUtils
using Test

xi = SE3(0)

Dx = SE3([1.0;0.0;0.0], Euler(0.0,0.0,0.0))

xj = xi ⊕ Dx
@test compare(xj, Dx)
@test compare(xi, xj ⊖ Dx)

println("Compare operations on two different SE3 transforms")

Dx1 = SE3(randn(3), convert(SO3, so3(0.25*randn(3))))
Dx2 = SE3(randn(3), convert(SO3, so3(0.25*randn(3))))

xi = SE3(0) ⊕ Dx1
xj = SE3(0) ⊕ Dx2
Dx = SE3(0) ⊖ xi ⊕ xj # inverse(xi) ⊕ xj

@show xi
@show xj
@show Dx
@test compare( xi, xj ⊖ Dx )
@test compare( xj, xi ⊕ Dx )
@test compare( Dx, inverse(xi) ⊕ xj)
@test compare( Dx, xi\xj )

#


text:
using DelimitedFiles
using Distributed
using LinearAlgebra
using ProgressMeter
using QuadGK
using SpecialFunctions
using Statistics

## Parameters
const ħ = 1
const ryd = 13.606              # Rydberg constant in eV
const a0 = 0.529                # Bohr radius in angstroms

const lattice_constant = 7.714  # In Bohr radii

const d1 = [lattice_constant / 2, lattice_constant * √(3) / 2]
const d2 = [-lattice_constant / 2, lattice_constant * √(3) / 2]

# Effective mass in the armchair(x) and zigzag(y) direction in m_e
const mx = 0.15
const my = 1.18

# Area of the unit cell in Bohr radii squared
const UC_area = d1[1] * d2[2] - d1[2] * d2[1];

# Upper limit of propagator integral
const C = 4 * π / UC_area / sqrt(mx * my) * ryd / 3;

const Rot_Mat = [cos(π / 3) sin(π / 3); -sin(π / 3) cos(π / 3)]
const M = [0, π / lattice_constant]

## Integration
const ν = 1e-4;         # Relative tolerance for integration
const η = 1e-6;         # Small number for imaginary parts

const NumEvals = 1e5;   # Max number of integrals in quadgk

## Location struct to label bulk unit cells.
struct Location
    v1::Int64   # Coefficient of d1
    v2::Int64   # Coefficient of d2
end

# Introduce addition and subtraction for the Location struct
# so we can use arithmetic operations on Locations

function Base.:+(l1::Location, l2::Location)
    return Location(l1.v1 + l2.v1, l1.v2 + l2.v2)
end

function Base.:-(l1::Location, l2::Location)
    return Location(l1.v1 - l2.v1, l1.v2 - l2.v2)
end

## Propagator
function P(r, z::ComplexF64)
    md = √(mx * r[1]^2 + my * r[2]^2)
    if md == 0.0
        return (-1 / C) * (log(Complex(1 - (C / z))))
    else
        return (-2 / C) * besselk(0, (md * √(-z / ryd)))
    end
end

function Ξ(r::Location, z::ComplexF64)
    position = r.v1 * d1 + r.v2 * d2
    res =
        P(position, z) * cos(dot(position, M)) +
        P(Rot_Mat' * position, z) * cos(dot(position, Rot_Mat * M)) +
        P(Rot_Mat * position, z) * cos(dot(position, Rot_Mat' * M))
    return res
end

function propagator_matrix(z, Coords::Vector{Location})
    len_coords = length(Coords)
    out = zeros(ComplexF64, len_coords, len_coords)
    for ii = 1:len_coords
        @inbounds for jj = ii:len_coords
            out[ii, jj] = Ξ(Coords[ii] - Coords[jj], z)
            out[jj, ii] = out[ii, jj]
        end
    end
    return out
end


text:
import Base: /
using ..YAXTools: PickAxisArray
struct YAXSlice{T,N,P}
    c::P
    sliceaxes
    otheraxes
end
function YAXSlice(p,dims)
    N = ndims(p)
    ax = caxes(p)
    isliceax = (findAxis.(dims,Ref(ax))...,)
    any(isnothing,isliceax) && error("Axis not in cube")
    iotherax = (filter(i->!in(i,isliceax),1:N)...,)
    YAXSlice{Array{eltype(p),length(dims)},N-length(dims),typeof(p)}(p,(isliceax,getindex.(Ref(ax),isliceax)), (iotherax,getindex.(Ref(ax),iotherax)))
end
/(c::YAXArray, s::Union{String,Symbol}) = YAXSlice(c,(s,))
/(c::YAXArray, s) = YAXSlice(c,s)

dimvals(x::YAXSlice, i) = x.c.axes[x.otheraxes[1][i]].values

function dimname(x::YAXSlice, i)
  axsym(x.otheraxes[2][i])
end

getattributes(x::YAXSlice) = x.c.properties

iscontdim(x::YAXSlice, i) = isa(x.c.axes[x.otheraxes[1][i]], RangeAxis)

function getdata(x::YAXSlice)
    m = fill!(Array{Union{Colon,Bool}}(undef,ndims(x.c)),true)
    m[collect(x.sliceaxes[1])] .= Colon()
    PickAxisArray(x.c.data, m)
end

getCubeDes(s::YAXSlice) = string(join(axname.(s.sliceaxes[2]), " x "), " slices over an ", getCubeDes(s.c))
cubesize(s::YAXSlice) = cubesize(s.c)
Base.ndims(s::YAXSlice{<:Any,N}) where N = N

Base.show(io::IO,s::YAXSlice) = YAXArrays.Cubes.show_yax(io,s)


text:
if !isdefined(Makie, :FigurePosition)
    const FigurePosition = Makie.GridPosition
    get_parent(x::Makie.GridPosition) = x.layout.parent
else
    using Makie: FigurePosition
    get_parent(x::FigurePosition) = x.fig
end

# We can't rely on inference to always give us fully typed
# Vector{<: Number} so we add `{T} where T` to the the mix
# This makes the number like type a bit absurd, but is still nice for
# documentation purposes!
const NumberLike = Union{Number,Missing,Nothing,T} where {T}
const NumberVector = AbstractVector{<:NumberLike}
const NumberMatrix = AbstractMatrix{<:NumberLike}

#####
##### Helpers for theming and color generation...May want to move them to Colors.jl / Makie.jl
#####
using Makie.Colors: LCHab, distinguishable_colors, RGB, Colorant

function get_theme(scene, key1::Symbol, key2::Symbol; defaults...)
    return get_theme(Makie.get_scene(scene), key1, key2; defaults...)
end

function get_theme(fig::FigurePosition, key1::Symbol, key2::Symbol; defaults...)
    return get_theme(get_parent(fig), key1, key2; defaults...)
end

# This function helps us to get the theme from a scene, that we can apply to our plotting functions
function get_theme(scene::Scene, key1::Symbol, key2::Symbol; defaults...)
    scene_theme = theme(scene)
    sub_theme = get(scene_theme, key1, Theme())
    # The priority is key1.key2 > defaults > key2
    # this way defaults are overwritten by theme options specifically set for our theme.
    # Consider Kappas.Axis for key1/2, what we want is, that if there are defaults in Kappas.Axis
    # they should overwrite our generic lighthouse defaults. But anything not specified in Kappas.Axis/defaults,
    # should fall back to scene_theme.Axis
    return merge(get(sub_theme, key2, Theme()), Theme(; defaults...),
                 get(scene_theme, key2, Theme()))
end

function high_contrast(background_color::Colorant, target_color::Colorant;
                       # chose from whole lightness spectrum
                       lchoices=range(0; stop=100, length=15))
    target = LCHab(target_color)
    color = distinguishable_colors(1, [RGB(background_color)]; dropseed=true, lchoices=lchoices,
                                   cchoices=[target.c], hchoices=[target.h])
    return RGBAf0(color[1], Makie.Colors.alpha(target_color))
end

"""
    Tuple{<:NumberVector, <: NumberVector}

Tuple of X, Y coordinates
"""
const XYVector = Tuple{<:NumberVector,<:NumberVector}

"""
    Union{XYVector, AbstractVector{<: XYVector}}

A series of XYVectors, or a single xyvector.
"""
const SeriesCurves = Union{XYVector,AbstractVector{<:XYVector}}

function series_plot!(subfig::FigurePosition, per_class_pr_curves::SeriesCurves,
                      class_labels::Union{Nothing,AbstractVector{String}}; legend=:lt, title="No title",
                      xlabel="x label", ylabel="y label", solid_color=nothing, color=nothing,
                      linewidth=nothing, scatter=NamedTuple())

    axis_theme = get_theme(subfig, :SeriesPlot, :Axis; title=title, titlealign=:left, xlabel=xlabel,
                           ylabel=ylabel, aspect=AxisAspect(1),
                           xticks=0:0.2:1, yticks=0.2:0.2:1)

    ax = Axis(subfig; axis_theme...)
    # Not the most elegant, but this way we can let the theming to the Series theme, or
    # pass it through explicitely
    series_theme = get_theme(subfig, :SeriesPlot, :Series; scatter...)
    isnothing(solid_color) || (series_theme[:solid_color] = solid_color)
    isnothing(color) || (series_theme[:color] = color)
    isnothing(linewidth) || (series_theme[:linewidth] = linewidth)
    series_theme = merge(series_theme, Attributes(;scatter...))
    hidedecorations!(ax; label=false, ticklabels=false, grid=false)
    limits!(ax, 0, 1, 0, 1)
    Makie.series!(ax, per_class_pr_curves; labels=class_labels, series_theme...)
    if !isnothing(legend)
        axislegend(ax; position=legend)
    end
    return ax
end

function plot_pr_curves!(subfig::FigurePosition, per_class_pr_curves::SeriesCurves,
                         class_labels::Union{Nothing,AbstractVector{String}}; legend=:lt, title="PR curves",
                         xlabel="True positive rate", ylabel="Precision", scatter=NamedTuple(),
                         solid_color=nothing)
    return series_plot!(subfig, per_class_pr_curves, class_labels; legend=legend, title=title, xlabel=xlabel,
                        ylabel=ylabel, scatter=scatter, solid_color=solid_color)
end

function plot_roc_curves!(subfig::FigurePosition, per_class_roc_curves::SeriesCurves,
                          per_class_roc_aucs::NumberVector, class_labels::AbstractVector{<:String};
                          legend=:rb, title="ROC curves", xlabel="False positive rate",
                          ylabel="True positive rate")
    auc_labels = [@sprintf("%s (AUC: %.3f)", class, per_class_roc_aucs[i])
                  for (i, class) in enumerate(class_labels)]

    return series_plot!(subfig, per_class_roc_curves, auc_labels; legend=legend, title=title, xlabel=xlabel,
                        ylabel=ylabel)
end

function plot_reliability_calibration_curves!(subfig::FigurePosition,
                                              per_class_reliability_calibration_curves::SeriesCurves,
                                              per_class_reliability_calibration_scores::NumberVector,
                                              class_labels::AbstractVector{String}; legend=:rb)
    calibration_score_labels = map(enumerate(class_labels)) do (i, class)
        @sprintf("%s (MSE: %.3f)", class, per_class_reliability_calibration_scores[i])
    end

    scatter_theme = get_theme(subfig, :ReliabilityCalibrationCurves, :Scatter; marker=Circle,
                              markersize=5, strokewidth=0)
    ideal_theme = get_theme(subfig, :ReliabilityCalibrationCurves, :Ideal; color=(:black, 0.5),
                            linestyle=:dash, linewidth=2)

    ax = series_plot!(subfig, per_class_reliability_calibration_curves, calibration_score_labels;
                      legend=legend, title="Prediction reliability calibration",
                      xlabel="Predicted probability bin", ylabel="Fraction of positives",
                      scatter=scatter_theme)
    #TODO: mean predicted value histogram underneath?? Maybe important...
    # https://scikit-learn.org/stable/modules/calibration.html
    linesegments!(ax, [0, 1], [0, 1]; ideal_theme...)
    return ax
end

function set_from_kw!(theme, key, kw, default)
    if haskey(kw, key)
        theme[key] = getproperty(kw, key)
    elseif !haskey(theme, key)
        theme[key] = default
    end
    return
end

function plot_binary_discrimination_calibration_curves!(subfig::FigurePosition,
                                                        calibration_curve::SeriesCurves, calibration_score,
                                                        per_expert_calibration_curves::SeriesCurves,
                                                        per_expert_calibration_scores, optimal_threshold,
                                                        discrimination_class::AbstractString; kw...)
    kw = values(kw)
    scatter_theme = get_theme(subfig, :BinaryDiscriminationCalibrationCurves, :Scatter; strokewidth=0)
    # Hayaah, this theme merging is getting out of hand
    # but we want kw > BinaryDiscriminationCalibrationCurves > Scatter, so we need to somehow set things
    # after the theme merging above, especially, since we also pass those to series!,
    # which then again tries to merge the kw args with the theme.
    set_from_kw!(scatter_theme, :makersize, kw, 5)
    set_from_kw!(scatter_theme, :marker, kw, :rect)

    per_expert = get_theme(subfig, :BinaryDiscriminationCalibrationCurves, :PerExpert; solid_color=:darkgrey,
                           color=nothing)
    set_from_kw!(per_expert, :linewidth, kw, 2)
    ax = series_plot!(subfig, per_expert_calibration_curves, nothing; legend=nothing,
                      title="Detection calibration", xlabel="Expert agreement rate",
                      ylabel="Predicted positive probability", scatter=scatter_theme, per_expert...)

    calibration = get_theme(subfig, :BinaryDiscriminationCalibrationCurves, :CalibrationCurve;
                            solid_color=:navyblue, markerstrokewidth=0)

    set_from_kw!(calibration, :markersize, kw, 5)
    set_from_kw!(calibration, :marker, kw, :rect)
    set_from_kw!(calibration, :linewidth, kw, 2)

    Makie.series!(ax, calibration_curve; calibration...)

    ideal_theme = get_theme(subfig, :BinaryDiscriminationCalibrationCurves, :Ideal; color=(:black, 0.5),
                            linestyle=:dash)
    set_from_kw!(ideal_theme, :linewidth, kw, 2)
    linesegments!(ax, [0, 1], [0, 1]; label="Ideal", ideal_theme...)
    #TODO: expert agreement histogram underneath?? Maybe important...
    # https://scikit-learn.org/stable/modules/calibration.html
    return ax
end

function plot_confusion_matrix!(subfig::FigurePosition, confusion::NumberMatrix,
                                class_labels::AbstractVector{String}, normalize_by::Symbol;
                                annotation_text_size=20, colormap=:Blues)
    normdim = get((Row=2, Column=1), normalize_by) do
        return error("normalize_by must be either :Row or :Column, found: $(normalize_by)")
    end

    nclasses = length(class_labels)
    if size(confusion) != (nclasses, nclasses)
        error("Labels must match size of square confusion matrix. Found $(nclasses) labels for an $(size(confusion)) matrix")
    end
    confusion = round.(confusion ./ sum(confusion; dims=normdim); digits=3)
    class_indices = 1:nclasses
    text_theme = get_theme(subfig, :ConfusionMatrix, :Text; textsize=annotation_text_size)
    heatmap_theme = get_theme(subfig, :ConfusionMatrix, :Heatmap; nan_color=(:black, 0.0))
    axis_theme = get_theme(subfig, :ConfusionMatrix, :Axis; xticklabelrotation=pi / 4,
                           titlealign=:left, title="$(string(normalize_by))-Normalized Confusion",
                           xlabel="Elected Class", ylabel="Predicted Class",
                           xticks=(class_indices, class_labels),
                           yticks=(class_indices, class_labels),
                           aspect=AxisAspect(1))

    ax = Axis(subfig; axis_theme...)

    hidedecorations!(ax; label=false, ticklabels=false, grid=false)
    ylims!(ax, nclasses, 0)
    tightlimits!(ax)
    plot_bg_color = to_color(ax.backgroundcolor[])
    crange = (0.0, 1.0)
    nan_color = to_color(heatmap_theme.nan_color[])
    cmap = to_colormap(to_value(pop!(heatmap_theme, :colormap, colormap)))
    heatmap!(ax, confusion'; colorrange=crange, colormap=cmap, nan_color=nan_color, heatmap_theme...)
    text_color = to_color(to_value(pop!(text_theme, :color, :black)))
    function label_color(i, j)
        c = confusion[i, j]
        bg_color = if isnan(c) || ismissing(c)
            Makie.Colors.alpha(nan_color) <= 0.0 ? plot_bg_color : nan_color
        else
            Makie.interpolated_getindex(cmap, c, crange)
        end
        return high_contrast(bg_color, text_color)
    end
    annos = vec([(string(confusion[i, j]), Point2f0(j, i)) for i in class_indices, j in class_indices])
    colors = vec([label_color(i, j) for i in class_indices, j in class_indices])
    text!(ax, annos; align=(:center, :center), color=colors, textsize=annotation_text_size, text_theme...)
    return ax
end

function text_attributes(values, groups, bar_colors, bg_color, text_color)
    aligns = NTuple{2,Symbol}[]
    offsets = NTuple{2,Int}[]
    text_colors = RGBAf0[]
    for (i, k) in enumerate(values)
        group = groups isa AbstractVector ? groups[i] : groups
        bar_color = bar_colors[group]
        # Plot text inside bar
        if k > 0.85
            push!(aligns, (:right, :center))
            push!(offsets, (-10, 0))
            push!(text_colors, high_contrast(bar_color, text_color))
        else
            # plot text next to bar
            push!(aligns, (:left, :center))
            push!(offsets, (10, 0))
            push!(text_colors, high_contrast(bg_color, text_color))
        end
    end
    return aligns, offsets, text_colors
end

function plot_kappas!(subfig::FigurePosition, per_class_kappas::NumberVector,
                      class_labels::AbstractVector{String}, per_class_IRA_kappas=nothing;
                      color=[:lightgrey, :lightblue], annotation_text_size=20)
    nclasses = length(class_labels)
    axis_theme = get_theme(subfig, :Kappas, :Axis; aspect=AxisAspect(1), titlealign=:left,
                           xlabel="Cohen's kappa", xticks=[0, 1], yreversed=true,
                           yticks=(1:nclasses, class_labels))

    text_theme = get_theme(subfig, :Kappas, :Text; textsize=annotation_text_size)
    text_color = to_color(to_value(pop!(text_theme, :color, to_color(:black))))
    bars_theme = get_theme(subfig, :Kappas, :BarPlot; color=color)
    bar_colors = to_color.(bars_theme.color[])

    ax = Axis(subfig[1, 1]; axis_theme...)
    bg_color = to_color(ax.backgroundcolor[])
    xlims!(ax, 0, 1)
    if isnothing(per_class_IRA_kappas)
        ax.title = "Algorithm-expert agreement"
        annotations = map(enumerate(per_class_kappas)) do (i, k)
            return (string(round(k; digits=3)), Point2f0(max(0, k), i))
        end
        aligns, offsets, text_colors = text_attributes(per_class_kappas, 2, bar_colors, bg_color, text_color)
        barplot!(ax, per_class_kappas; direction=:x, color=bar_colors[2])
        text!(ax, annotations; align=aligns, offset=offsets, color=text_colors, text_theme...)
    else
        ax.title = "Inter-rater reliability"
        values = vcat(per_class_kappas, per_class_IRA_kappas)
        groups = vcat(fill(2, nclasses), fill(1, nclasses))
        xvals = vcat(1:nclasses, 1:nclasses)
        cmap = bar_colors
        bars = barplot!(ax, xvals, max.(0, values); dodge=groups, color=groups, direction=:x, colormap=cmap)
        # This is a bit hacky, but for now the easiest way to figure out the exact, dodged positions
        rectangles = bars.plots[][1][]
        dodged_y = last.(minimum.(rectangles)) .+ (last.(widths.(rectangles)) ./ 2)
        textpos = Point2f0.(max.(0, values), dodged_y)

        labels = string.(round.(values; digits=3))
        aligns, offsets, text_colors = text_attributes(values, groups, bar_colors, bg_color, text_color)
        text!(ax, labels; position=textpos, align=aligns, offset=offsets, color=text_colors, text_theme...)
        labels = ["Expert-vs-expert IRA", "Algorithm-vs-expert"]
        entries = map(c -> PolyElement(; color=c, strokewidth=0, strokecolor=:white), cmap)
        legend = Legend(subfig[1, 1, Bottom()], entries, labels; tellwidth=false, tellheight=true,
                        labelsize=12, padding=(0, 0, 0, 0), framevisible=false, patchsize=(10, 10),
                        patchlabelgap=6)
        legend.margin = (0, 0, 0, 60)
    end
    return ax
end

"""
    evaluation_metrics_plot(data::Dict; resolution=(1000, 1000), textsize=12)

Plots all evaluation metrics Lighthouse has to offer.
"""
function evaluation_metrics_plot(data::Dict; resolution=(1000, 1000), textsize=12)
    fig = Figure(; resolution=resolution, Axis=(titlesize=17,))

    # Confusion
    plot_confusion_matrix!(fig[1, 1], data["confusion_matrix"], data["class_labels"], :Column;
                           annotation_text_size=textsize)
    plot_confusion_matrix!(fig[1, 2], data["confusion_matrix"], data["class_labels"], :Row;
                           annotation_text_size=textsize)
    # Kappas
    IRA_kappa_data = nothing
    multiclass = length(data["class_labels"]) > 2
    labels = multiclass ? vcat("Multiclass", data["class_labels"]) : data["class_labels"]
    kappa_data = multiclass ? vcat(data["multiclass_kappa"], data["per_class_kappas"]) :
                 data["per_class_kappas"]

    if issubset(["multiclass_IRA_kappas", "per_class_IRA_kappas"], keys(data))
        IRA_kappa_data = multiclass ? vcat(data["multiclass_IRA_kappas"], data["per_class_IRA_kappas"]) :
                         data["per_class_IRA_kappas"]
    end

    plot_kappas!(fig[1, 3], kappa_data, labels, IRA_kappa_data; annotation_text_size=textsize)

    # Curves
    ax = plot_roc_curves!(fig[2, 1], data["per_class_roc_curves"], data["per_class_roc_aucs"],
                          data["class_labels"]; legend=nothing)

    plot_pr_curves!(fig[2, 2], data["per_class_pr_curves"], data["class_labels"]; legend=nothing)

    plot_reliability_calibration_curves!(fig[3, 1], data["per_class_reliability_calibration_curves"],
                                         data["per_class_reliability_calibration_scores"],
                                         data["class_labels"]; legend=nothing)

    legend_pos = 2:3
    if haskey(data, "discrimination_calibration_curve")
        legend_pos = 3
        plot_binary_discrimination_calibration_curves!(fig[3, 2], data["discrimination_calibration_curve"],
                                                       data["discrimination_calibration_score"],
                                                       data["per_expert_discrimination_calibration_curves"],
                                                       data["per_expert_discrimination_calibration_scores"],
                                                       data["optimal_threshold"],
                                                       data["class_labels"][data["optimal_threshold_class"]])
    end
    legend_plots = filter(Makie.MakieLayout.get_plots(ax)) do plot
        return haskey(plot, :label)
    end
    elements = map(legend_plots) do elem
        return [PolyElement(; color=elem.color, strokecolor=:transparent)]
    end

    function label_str(i)
        auc = round(data["per_class_roc_aucs"][i]; digits=2)
        mse = round(data["per_class_reliability_calibration_scores"][i]; digits=2)
        return ["""ROC AUC  $auc
                   Cal. MSE    $mse
                   """]
    end
    classes = data["class_labels"]
    nclasses = length(classes)
    class_labels = label_str.(1:nclasses)
    Legend(fig[3, legend_pos], elements, class_labels, classes; nbanks=2, tellwidth=false, tellheight=false,
           labelsize=11, titlesize=14, titlegap=5, groupgap=6, labelhalign=:left, labelvalign=:center)
    colgap!(fig.layout, 2)
    rowgap!(fig.layout, 4)
    return fig
end

# Helper to more easily define the non mutating versions
function axisplot(func, args; resolution=(800, 600), plot_kw...)
    fig = Figure(; resolution=resolution)
    ax = func(fig[1, 1], args...; plot_kw...)
    # ax.plots[1] is not really that great, but there isn't a FigureAxis object right now
    # this will need to wait for when we figure out a better recipe integration
    return Makie.FigureAxisPlot(fig, ax, ax.scene.plots[1])
end

"""
    plot_confusion_matrix!(subfig::FigurePosition, args...; kw...)

    plot_confusion_matrix(confusion::AbstractMatrix{<: Number}, class_labels::AbstractVector{String}, normalize_by::Symbol;
                          resolution=(800,600),
                          annotation_text_size=20)

Lighthouse plots confusion matrices, which are simple tables
showing the empirical distribution of predicted class (the rows)
versus the elected class (the columns). These come in two variants:

* row-normalized: this means each row has been normalized to sum to 1. Thus, the row-normalized confusion matrix shows the empirical distribution of elected classes for a given predicted class. E.g. the first row of the row-normalized confusion matrix shows the empirical probabilities of the elected classes for a sample which was predicted to be in the first class.
* column-normalized: this means each column has been normalized to sum to 1. Thus, the column-normalized confusion matrix shows the empirical distribution of predicted classes for a given elected class. E.g. the first column of the column-normalized confusion matrix shows the empirical probabilities of the predicted classes for a sample which was elected to be in the first class.

```
fig, ax, p = plot_confusion_matrix(rand(2, 2), ["1", "2"], :Row)
fig = Figure()
ax = plot_confusion_matrix!(fig[1, 1], rand(2, 2), ["1", "2"], :Column)
```
"""
plot_confusion_matrix(args...; kw...) = axisplot(plot_confusion_matrix!, args; kw...)

"""
    plot_reliability_calibration_curves!(fig::SubFigure, args...; kw...)

    plot_reliability_calibration_curves(per_class_reliability_calibration_curves::SeriesCurves,
                                        per_class_reliability_calibration_scores::NumberVector,
                                        class_labels::AbstractVector{String};
                                        legend=:rb, resolution=(800, 600))

"""
function plot_reliability_calibration_curves(args...; kw...)
    return axisplot(plot_reliability_calibration_curves!, args; kw...)
end

"""
    plot_binary_discrimination_calibration_curves!(fig::SubFigure, args...; kw...)

    plot_binary_discrimination_calibration_curves!(calibration_curve::SeriesCurves, calibration_score,
                                                   per_expert_calibration_curves::SeriesCurves,
                                                   per_expert_calibration_scores, optimal_threshold,
                                                   discrimination_class::AbstractString;
                                                   marker=:rect, markersize=5, linewidth=2)

"""
function plot_binary_discrimination_calibration_curves(args...; kw...)
    return axisplot(plot_binary_discrimination_calibration_curves!, args; kw...)
end

"""
    plot_pr_curves!(subfig::FigurePosition, args...; kw...)

    plot_pr_curves(per_class_pr_curves::SeriesCurves,
                class_labels::AbstractVector{<: String};
                resolution=(800, 600),
                legend=:lt, title="PR curves",
                xlabel="True positive rate", ylabel="Precision",
                linewidth=2, scatter=NamedTuple(), color=:darktest)

- `scatter::Union{Nothing, NamedTuple}`: can be set to a named tuples of attributes that are forwarded to the scatter call (e.g. markersize). If nothing, no scatter is added.

"""
plot_pr_curves(args...; kw...) = axisplot(plot_pr_curves!, args; kw...)

"""
    plot_roc_curves!(subfig::FigurePosition, args...; kw...)

    plot_roc_curves(per_class_roc_curves::SeriesCurves,
                    per_class_roc_aucs::NumberVector,
                    class_labels::AbstractVector{<: String};
                    resolution=(800, 600),
                    legend=:lt,
                    title="ROC curves",
                    xlabel="False positive rate",
                    ylabel="True positive rate",
                    linewidth=2, scatter=NamedTuple(), color=:darktest)

- `scatter::Union{Nothing, NamedTuple}`: can be set to a named tuples of attributes that are forwarded to the scatter call (e.g. markersize). If nothing, no scatter is added.

"""
plot_roc_curves(args...; kw...) = axisplot(plot_roc_curves!, args; kw...)

"""
    plot_kappas!(subfig::FigurePosition, args...; kw...)

    plot_kappas(per_class_kappas::NumberVector,
                class_labels::AbstractVector{String},
                per_class_IRA_kappas=nothing;
                resolution=(800, 600),
                annotation_text_size=20)
"""
plot_kappas(args...; kw...) = axisplot(plot_kappas!, args; kw...)


text:
NumTypes = Union{Number, Vector{Float32}, Vector{Float64}} # TODO use where T<>:Number in 0.6

#TODO add to Erdos api
function test_networks_eq(g::ANetOrDiNet, h::ANetOrDiNet)
    @test h == g

    @test sort(gprop_names(h)) == sort(gprop_names(g))
    @test sort(vprop_names(h)) == sort(vprop_names(g))
    @test sort(eprop_names(h)) == sort(eprop_names(g))

    for (pname, ph) in gprop(h)
        @test has_gprop(g, pname)
        # @test typeof(ph) == typeof(gprop(g, pname))
        if  typeof(ph) <: NumTypes
            @test ph ≈ gprop(g, pname)
        else
            @test ph == gprop(g, pname)
        end
    end

    for (pname, ph) in vprop(h)
        @test has_vprop(g, pname)
        pg = vprop(g, pname)
        # @test valtype(pg) == valtype(ph)
        for i=1:nv(g)
            @test  haskey(ph, i) == haskey(pg, i)
            haskey(ph, i) || continue
            if valtype(ph) <: NumTypes
                @test ph[i] ≈ pg[i]
            else
                @test ph[i] == pg[i]
            end
        end
    end

    for (pname, ph) in eprop(h)
        @test has_eprop(g, pname)
        pg = eprop(g, pname)
        # @test valtype(pg) == valtype(ph)
        for e in edges(h)
            i, j = src(e), dst(e)
            @test  haskey(ph, e) == haskey(pg, i, j)
            haskey(ph, e) || continue
            if valtype(ph) <: NumTypes
                @test ph[e] ≈ pg[i,j]
            else
                @test ph[e] == pg[i,j]
            end
        end
    end
end


text:
# FIXME temporary copy/paste from:
# https://github.com/FluxML/MacroTools.jl/blob/master/src/utils.jl#L55-L87

isexpr(x::Expr) = true
isexpr(x) = false
isexpr(x::Expr, ts...) = x.head in ts
isexpr(x, ts...) = any(T->isa(T, Type) && isa(x, T), ts)

isline(ex) = isexpr(ex, :line) || isa(ex, LineNumberNode)

rmlines(x) = x
function rmlines(x::Expr)
  # Do not strip the first argument to a macrocall, which is
  # required.
  if x.head == :macrocall && length(x.args) >= 2
    Expr(x.head, x.args[1], nothing, filter(x->!isline(x), x.args[3:end])...)
  else
    Expr(x.head, filter(x->!isline(x), x.args)...)
  end
end

striplines(ex) = prewalk(rmlines, ex)

walk(x, inner, outer) = outer(x)
walk(x::Expr, inner, outer) = outer(Expr(x.head, map(inner, x.args)...))
prewalk(f, x)  = walk(f(x), x -> prewalk(f, x), identity)