text:
### ADDING  ###

"""
`add!(H::SimpleHypergraph{T}, v::T)` to add a vertex to `H`. **Example**:
`add!(H,3)`.

`add!(H::SimpleHypergraph{T}, e::Set{T})` to add a hyperedge. This
also works if `e` is a one-dimensional array of `T`s.
**Example**: `add!(H,[1,2,3])`.
"""
function add!(H::SimpleHypergraph{T}, v::T)::Bool where {T}
    if in(v, H.V)
        return false
    end
    push!(H.V, v)
    H.VE[v] = Set{Set{T}}()
    return true
end

function add!(H::SimpleHypergraph{T}, e::Set{T})::Bool where {T}
    if in(e, H.E)
        return false  # already have this edge
    end

    # make sure we have all the vertices
    for v in e
        add!(H, v)
    end

    # add the edge
    push!(H.E, e)

    # update the VE dictionary
    for v in e
        push!(H.VE[v], e)
    end

    return true
end

function add!(H::SimpleHypergraph{T}, e::Vector{T})::Bool where {T}
    add!(H, Set(e))
end

function add!(H::SimpleHypergraph{T}, e::T...) where {T}
    add!(H, Set(e))
end

### DELETING ###

"""
`delete!(H::SimpleHypergraph,v)` deletes vertex `v` from `H` as well as
all edges that contain `v`.

`delete!(H::SimpleHypergraph,e)` deletes edge `e` from `H`. Here,
`e` is either a `Set` or a `Vector` of vertices.
"""
function delete!(H::SimpleHypergraph{T}, e::Set{T})::Bool where {T}
    if !has(H, e)   # nothing to do if this edge isn't in H
        return false
    end

    # Delete e from H.E
    delete!(H.E, e)

    # for every v in e, delete e from VE[v]
    for v in e
        delete!(H.VE[v], e)
    end

    return true
end

function delete!(H::SimpleHypergraph{T}, e::Vector{T})::Bool where {T}
    delete!(H, Set(e))
end

function delete!(H::SimpleHypergraph{T}, e::T...) where {T}
    delete!(H, Set(e))
end

function delete!(H::SimpleHypergraph{T}, v::T)::Bool where {T}
    if !has(H, v)
        return false
    end

    # delete all edges that contain v
    for e in H[v]
        delete!(H, e)
    end

    # delete v from H.V
    delete!(H.V, v)

    return true
end


text:
struct RootedTree{T<:Integer}
    ancestors::Dict{T,T}
    root::T
end

root(T::RootedTree) = T.root
ancestor(T::RootedTree{S},i::S) where S<:Integer = T.ancestors[i]
vertices(T::RootedTree) = keys(T.ancestors)∪[root(T)]

function ancestormap(T::RootedTree{S}) where S<:Integer
    D = Dict{S,Set{S}}(k=>Set() for k in vertices(T))
    for (child,ancestor) in T.ancestors
        push!(D[ancestor],child)
    end
    return D
end

function treesum(T::RootedTree{S},f::Dict{S,U}) where {S<:Integer, U<:Real}
    D = Dict{S,U}(k=>zero(U) for k in vertices(T))
    A = ancestormap(T)
    R = root(T)
    function updatebranch!(k)
        D[k] += f[k]
        if k ≠ R
            D[k] += D[ancestor(T,k)]
        end
        foreach(updatebranch!,A[k])
    end
    updatebranch!(R)
    return D
end

function backsum(T::RootedTree{S},f::Dict{S,U}) where {S<:Integer, U<:Real}
    D = Dict{S,U}()
    A = ancestormap(T)
    function updatevertex!(k)
        if k in keys(D)
            return D[k]
        else
            D[k] = f[k] + reduce(+, [updatevertex!(v) for v in A[k]], init = 0)
        end
        D[k]
    end
    updatevertex!(root(T))
    D
end

import Base: +, -
+(D::Dict{<:Any,<:Real},t::Real) = Dict(k=>v+t for (k,v) in D)
-(D::Dict{<:Any,<:Real},t::Real) = D + (-t)

struct WoodedTriangulation{T<:Integer}
    P::PlanarMap{T}
    bluetree::RootedTree{T}
    redtree::RootedTree{T}
    greentree::RootedTree{T}
end

PlanarMap(W::WoodedTriangulation) = W.P

length(WT::WoodedTriangulation) = length(WT.P)

function show(io::IO,WT::WoodedTriangulation)
    print(io,"WoodedTriangulation(")
    print(io,string(length(WT.P)))
    print(io,")")
end

outerface(WT::WoodedTriangulation) =
    map(root,(WT.bluetree,WT.redtree,WT.greentree))

interiorvertices(WT) = setdiff(1:length(WT.P),outerface(WT))

blueroot(WT::WoodedTriangulation) = root(WT.bluetree)
redroot(WT::WoodedTriangulation) = root(WT.redtree)
greenroot(WT::WoodedTriangulation) = root(WT.greentree)

isblue(WT::WoodedTriangulation{T},i::T,j::T) where T<:Integer =
    i ∉ outerface(WT) && ancestor(WT.bluetree,i) == j
isred(WT::WoodedTriangulation{T},i::T,j::T) where T<:Integer =
    i ∉ outerface(WT) && ancestor(WT.redtree,i) == j
isgreen(WT::WoodedTriangulation{T},i::T,j::T) where T<:Integer =
    i ∉ outerface(WT) && ancestor(WT.greentree,i) == j

struct ColoredEdge
    tail::Int64
    head::Int64
    color::Symbol
end

function edges(WT::WoodedTriangulation)
    interior_edges = ColoredEdge[]
    for (d,c) in zip((WT.bluetree,WT.redtree,WT.greentree),(:blue,:red,:green))
        append!(interior_edges,vcat([[ColoredEdge(k,d[k],c),
                                 ColoredEdge(d[k],k,c)] for k in keys(d)]...))
    end
    roots = [f(WT) for f in (blueroot,redroot,greenroot)]
    vcat(interior_edges,
    [ColoredEdge(a,b,:gray) for a in roots for b in roots if a ≠ b])
end

function faces(WT::WoodedTriangulation{T}) where T<:Integer
    n = length(WT.P)-3
    all_edges = edges(WT)
    colors = Dict((E.tail,E.head)=>E.color for E in all_edges);
    interior_edges = filter(E->E.color ≠ :gray, all_edges)
    found_edges = Set{NTuple{2,T}}()
    all_faces = Face{T}[]
    u = blueroot(WT)
    v = cw(WT.P,u,greenroot(WT))
    while length(found_edges) < length(interior_edges)
        if colors[(u,v)] == :blue
            (u,v) = v,cw(WT.P,v,u)
        else
            (u,v) = u,cw(WT.P,u,v)
        end
        if (u,v) in found_edges
            continue
        end
        f = face(WT.P,u,v)
        push!(found_edges,edges(f)...)
        push!(all_faces,f)
    end
    return all_faces
end

"""
    schnyderwood(P::PlanarMap,
                 outface=face(P,1,neighbors(P,1)[1]))

Find a Schnyder wood associated with the plane
triangulation P, with given outer face
"""
function schnyderwood(RNG::Random.AbstractRNG,
                      P::PlanarMap{T};
                      outface::Face=face(P,1,neighbors(P,1)[1])) where T<: Integer
    PWT = PartialWoodedTriangulation(P;outface=outface)
    while !allfound(PWT)
        v = rand(RNG,PWT.shellcandidates)
        shell!(PWT,v)
    end
    return WoodedTriangulation(PWT)
end

function schnyderwood(P::PlanarMap{T};
                      outface::Face=face(P,1,neighbors(P,1)[1])) where T<: Integer
      schnyderwood(Random.GLOBAL_RNG,P,outface=outface)
end

struct PartialWoodedTriangulation{T<:Integer}
    P::PlanarMap{T}
    bluetree::Dict{T,T}
    redtree::Dict{T,T}
    greentree::Dict{T,T}
    blueroot::T
    redroot::T
    greenroot::T
    discovered::Set{T}
    shellcandidates::Set{T}
    fneighborcount::Vector{Int64}
end

function PartialWoodedTriangulation(P::PlanarMap{T};
                        outface::Face=face(P,1,neighbors(P,1)[1])) where T<:Integer
    blueroot, greenroot, redroot = outface # bgr order not a mistake
    bluetree = Dict{T,T}()
    redtree = Dict{T,T}()
    greentree = Dict{T,T}()
    discovered = Set([blueroot])
    shellcandidates = Set([blueroot])
    affectedvertices = Set{T}()
    # number of frontier neighbors of each vertex in P
    fneighborcount = zeros(Int64,length(P))
    for v in outface
        for w in neighbors(P,v)
            fneighborcount[w] += 1
        end
    end
    PartialWoodedTriangulation(
    P,bluetree,redtree,greentree,blueroot,redroot,greenroot,
    discovered,shellcandidates,fneighborcount
    )
end

function WoodedTriangulation(PWT::PartialWoodedTriangulation)
    if !allfound(PWT)
        error("Partial wooded triangulation not complete")
    end
    WoodedTriangulation(PWT.P,
            map(RootedTree,(PWT.bluetree, PWT.redtree, PWT.greentree),
                           (PWT.blueroot, PWT.redroot, PWT.greenroot))...)
end

function allfound(PWT::PartialWoodedTriangulation)
    length(PWT.P) == length(PWT.discovered) + 2
end

function shell!(PWT::PartialWoodedTriangulation{T},v::T) where T<:Integer
    push!(PWT.discovered,v)
    affectedvertices = Set{T}()
    if v == PWT.blueroot
        N = collect(rotate(neighbors(PWT.P,v),PWT.redroot))
    else
        N = downwardneighbors(neighbors(PWT.P,v),PWT.bluetree[v])
        filter!(x->x ∉ PWT.discovered,N)
        PWT.redtree[v] = N[1]
        PWT.greentree[v] = N[end]
    end
    for w in N[2:end-1]
        PWT.bluetree[w] = v
    end
    # update frontier neighbor count
    for w in neighbors(PWT.P,v)
        PWT.fneighborcount[w] -= 1
        push!(affectedvertices,w)
    end
    for u in N[2:end-1]
        for w in neighbors(PWT.P,u)
            PWT.fneighborcount[w] += 1
            push!(affectedvertices,w)
        end
    end
    # update active vertex list
    pop!(PWT.shellcandidates,v)
    for w in N
        if PWT.fneighborcount[w] == 2 && w ≠ PWT.redroot && w ≠ PWT.greenroot
            push!(PWT.shellcandidates,w)
        end
    end
    for w in affectedvertices
        if PWT.fneighborcount[w] ≠ 2 && w in PWT.shellcandidates
            pop!(PWT.shellcandidates,w)
        end
    end
    if 0 in PWT.shellcandidates
        error("hey!")
    end
end

function downwardneighbors(N::NeighborCycle{T},v::T) where T<:Integer
    p = position(N,v)
    N[to(p+1,p-1)]
end

"""
    schnyder_red(WT::WoodedTriangulation)

Walk the green tree to determine the red coordinate
(the number of faces between the green and blue flow lines
from each vertex)
"""
function schnyder_red(WT::WoodedTriangulation{T}) where T<:Integer
    # RD = red descendants, GA = green ancestors, BA = blue ancestors
    RD = backsum(WT.redtree,Dict(k=>1 for k in vertices(WT.redtree)))-1
    GA = treesum(WT.greentree,Dict(k=>1 for k in vertices(WT.greentree)))-1
    BA = treesum(WT.bluetree,Dict(k=>1 for k in vertices(WT.bluetree)))-1
    # BR = sum of red descendants with respect to blue tree
    # GR = sum of red descendants with respect to
    RD[blueroot(WT)] = 0; RD[greenroot(WT)] = 0
    BR = treesum(WT.bluetree,RD)
    GR = treesum(WT.greentree,RD)
    Dict(k=>2(BR[k] + GR[k] - RD[k]) + GA[k] + BA[k] - 1 for k in interiorvertices(WT))
end

function schnyder_green(WT::WoodedTriangulation{T}) where T<:Integer
    schnyder_red(WoodedTriangulation(WT.P,WT.redtree,WT.greentree,WT.bluetree))
end

function schnyder_blue(WT::WoodedTriangulation{T}) where T<:Integer
    schnyder_red(WoodedTriangulation(WT.P,WT.greentree,WT.bluetree,WT.redtree))
end

function schnyder_coords(WT::WoodedTriangulation)
    redcoords = schnyder_red(WT)
    greencoords = schnyder_green(WT)
    D = Dict(k=>(redcoords[k],greencoords[k]) for k in keys(redcoords))
    n = length(WT)
    merge!(D,Dict(blueroot(WT)=>(0,0),redroot(WT)=>(2n-5,0),greenroot(WT)=>(0,2n-5)))
    [D[k] for k=1:length(D)]
end

"""
    schnyderaverage(P::PlanarMap,
                    outface::Face,
                    n::Integer=100)
"""
function schnyderaverage(RNG::Random.AbstractRNG,
                         P::PlanarMap{T};
                         outface::Face=face(P,1,neighbors(P,1)[1]),
                         n::Integer=100) where T<: Integer
    N = length(P)
    A = zeros(N)
    for i=1:n
        D = schnyder_coords(schnyderwood(RNG,P,outface=outface))
        A = [D[i] .+ A[i] for i=1:N]
    end
    [a./n for a in A]
end

function schnyderaverage(P::PlanarMap;kwargs...)
    schnyderaverage(Random.GLOBAL_RNG,P;kwargs...)
end

"""
Non-mutating version of shell!
"""
function shell(P::PartialWoodedTriangulation{T},v::T) where T<:Integer
    Q = deepcopy(P)
    shell!(Q,v)
    Q
end

"""
    allwoods(PWT::PartialWoodedTriangulation,C::Channel)

Send to `C` a sequence of all `WoodedTriangulation`s on `P`
obtained by further shelling `PWT`
"""
function allwoods(PWT::PartialWoodedTriangulation,C::Channel)
    if allfound(PWT)
        put!(C,WoodedTriangulation(PWT))
    else
        for v in copy(PWT.shellcandidates)
            allwoods(shell(PWT,v),C)
        end
    end
end

"""
    allwoods(P::PlanarMap,C::Channel)

Send a sequence of all `WoodedTriangulation`s on `P`
to the channel `C`

# Examples:

```julia
# Count the number of Schnyder woods on a
# a randomly sampled wooded triangulation of size 10
sum(1 for W in Channel(C->allwoods(PlanarMap(UWT(10)),C)))
```
"""
function allwoods(P::PlanarMap,C::Channel;
                  outface::Face=face(P,1,neighbors(P,1)[1]))
    allwoods(PartialWoodedTriangulation(P;outface=outface),C)
end

ext: jl


text:
using ModelingToolkit, MethodOfLines, DomainSets, Test, Symbolics, SymbolicUtils

@testset "Test 01: discretization of variables, center aligned grid" begin
    # Test centered order

    @parameters x, y, t
    @variables u(..) v(..)

    t_min= 0.
    t_max = 2.0
    x_min = 0.
    x_max = 2.
    y_min = 0.
    y_max = 2.
    dx = 0.1; dy = 0.2
    order = 2

    domains = [t ∈ Interval(t_min, t_max), x ∈ Interval(x_min, x_max), y ∈ Interval(y_min,y_max)]

    pde = u(t, x, y) ~ v(t, x, y) + Differential(x)(u(t, x, y)) + Differential(y)(v(t, x, y))

    bcs = [u(t_min, x, y) ~ 0, u(t_max, x, y) ~ 0, v(t_min, x, y) ~ 0, v(t_max, x, y) ~ 0]

    @named pdesys = PDESystem(pde,bcs,domains,[t,x,y],[u(t,x,y),v(t,x,y)])

    depvar_ops = map(x->operation(x.val),pdesys.depvars)

    depvars_lhs = MethodOfLines.get_depvars(pde.lhs, depvar_ops)
    depvars_rhs = MethodOfLines.get_depvars(pde.rhs, depvar_ops)
    depvars = collect(depvars_lhs ∪ depvars_rhs)
    # Read the independent variables,
    # ignore if the only argument is [t]
    indvars = first(Set(filter(xs->!isequal(xs, [t]), map(arguments, depvars))))
    x̄ = first(Set(filter(!isempty, map(u->filter(x-> t === nothing || !isequal(x, t.val), arguments(u)), depvars))))

    disc = MOLFiniteDifference([x=>dx, y=>dy], t; approx_order=order)

    s = MethodOfLines.DiscreteSpace(domains, depvars, x̄, disc)

    discx = x_min:dx:x_max
    discy = y_min:dy:y_max

    grid = Dict([x=>discx, y=>discy])

    @test s.grid[x] == discx
    @test s.grid[y] == discy

    @test s.axies[x] == s.grid[x]
    @test s.axies[y] == s.grid[y]

    #@test all([all(I[i] .∈ (collect(s.Igrid),)) for I in values(s.Iedges), i in [1,2]])

    @test all(s.Iaxies[u] == s.Igrid[u] for u in depvars)

end

@testset "Test 02: discretization of variables, edge aligned grid" begin
    # Test centered order
    @parameters x, y, t
    @variables u(..) v(..)

    t_min= 0.
    t_max = 2.0
    x_min = 0.
    x_max = 2.
    y_min = 0.
    y_max = 2.
    dx = 0.1; dy = 0.2
    order = 2

    domains = [t ∈ Interval(t_min, t_max), x ∈ Interval(x_min, x_max), y ∈ Interval(y_min,y_max)]

    pde = u(t, x, y) ~ v(t, x, y) + Differential(x)(u(t, x, y)) + Differential(y)(v(t, x, y))
    bcs = [u(t_min, x, y) ~ 0, u(t_max, x, y) ~ 0, v(t_min, x, y) ~ 0, v(t_max, x, y) ~ 0]

    @named pdesys = PDESystem(pde,bcs,domains,[t,x,y],[u(t,x,y),v(t,x,y)])

    depvar_ops = map(x->operation(x.val),pdesys.depvars)

    depvars_lhs = MethodOfLines.get_depvars(pde.lhs, depvar_ops)
    depvars_rhs = MethodOfLines.get_depvars(pde.rhs, depvar_ops)
    depvars = collect(depvars_lhs ∪ depvars_rhs)
    # Read the independent variables,
    # ignore if the only argument is [t]
    indvars = first(Set(filter(xs->!isequal(xs, [t]), map(arguments, depvars))))
    x̄ = first(Set(filter(!isempty, map(u->filter(x-> t === nothing || !isequal(x, t.val), arguments(u)), depvars))))

    disc = MOLFiniteDifference([x=>dx, y=>dy], t; approx_order=order, grid_align=edge_align)

    s = MethodOfLines.DiscreteSpace(domains, depvars, x̄, disc)

    discx = (x_min-dx/2): dx : (x_max+dx/2)
    discy = (y_min-dy/2): dy : (y_max+dy/2)

    grid = Dict([x=>discx, y=>discy])

    @test s.grid[x] == discx
    @test s.grid[y] == discy

    @test s.axies[x] != s.grid[x]
    @test s.axies[y] != s.grid[y]

    #@test all([all(I[i] .∈ (collect(s.Igrid),)) for I in values(s.Iedges), i in [1,2]])

    @test all(s.Iaxies[u] != s.Igrid[u] for u in depvars)

end


text:
import AutoryBroadcastMacros
AutoryBroadcastMacros.runtests()


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

mutable struct VirtualNetworkRule <: SwaggerModel
    id::Any # spec type: Union{ Nothing, String } # spec name: id
    action::Any # spec type: Union{ Nothing, String } # spec name: action
    state::Any # spec type: Union{ Nothing, String } # spec name: state

    function VirtualNetworkRule(;id=nothing, action="Allow", state=nothing)
        o = new()
        validate_property(VirtualNetworkRule, Symbol("id"), id)
        setfield!(o, Symbol("id"), id)
        validate_property(VirtualNetworkRule, Symbol("action"), action)
        setfield!(o, Symbol("action"), action)
        validate_property(VirtualNetworkRule, Symbol("state"), state)
        setfield!(o, Symbol("state"), state)
        o
    end
end # type VirtualNetworkRule

const _property_map_VirtualNetworkRule = Dict{Symbol,Symbol}(Symbol("id")=>Symbol("id"), Symbol("action")=>Symbol("action"), Symbol("state")=>Symbol("state"))
const _property_types_VirtualNetworkRule = Dict{Symbol,String}(Symbol("id")=>"String", Symbol("action")=>"String", Symbol("state")=>"String")
Base.propertynames(::Type{ VirtualNetworkRule }) = collect(keys(_property_map_VirtualNetworkRule))
Swagger.property_type(::Type{ VirtualNetworkRule }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_VirtualNetworkRule[name]))}
Swagger.field_name(::Type{ VirtualNetworkRule }, property_name::Symbol) =  _property_map_VirtualNetworkRule[property_name]

const _allowed_VirtualNetworkRule_action = ["Allow"]

const _allowed_VirtualNetworkRule_state = ["provisioning", "deprovisioning", "succeeded", "failed", "networkSourceDeleted"]

function check_required(o::VirtualNetworkRule)
    (getproperty(o, Symbol("id")) === nothing) && (return false)
    true
end

function validate_property(::Type{ VirtualNetworkRule }, name::Symbol, val)
    if name === Symbol("action")
        Swagger.validate_param(name, "VirtualNetworkRule", :enum, val, _allowed_VirtualNetworkRule_action)
    end
    if name === Symbol("state")
        Swagger.validate_param(name, "VirtualNetworkRule", :enum, val, _allowed_VirtualNetworkRule_state)
    end
end


text:
module TestConversion
using StaticArrays
using Test
using PhaseSpaceIO
using PhaseSpaceIO.Testing
using PhaseSpaceIO: position
using PhaseSpaceIO.Conversion: to_egs, to_iaea, phsp_convert
using Accessors

@testset "to_egs, to_iaea" begin
    for _ in 1:100
        Ni, Nf = rand(0:3,2)
        p_iaea = arbitrary(IAEAParticle{Nf, Ni})
        @reset p_iaea.typ = rand([electron, positron, photon])
        p_egs = to_egs(p_iaea)
        p_iaea2 = to_iaea(p_egs,z=p_iaea.z)
        @test position(p_egs) == position(p_iaea)[1:2]
        @test direction(p_egs) ≈ direction(p_iaea)
        @test particle_type(p_egs) == particle_type(p_iaea)
        @test p_egs.E ≈ p_iaea.E
        @test p_egs.new_history == p_iaea.new_history
        @test p_egs.weight == p_iaea.weight

        @test position(p_iaea) ≈ position(p_iaea2)
        @test direction(p_iaea) ≈ direction(p_iaea2)
        @test particle_type(p_iaea2) == particle_type(p_iaea)
        @test p_iaea2.E ≈ p_iaea.E
        @test p_iaea2.new_history == p_iaea.new_history
        @test p_egs.weight == p_iaea.weight

        p_iaea_bad = @set p_iaea.typ = rand([proton, neutron])
        @test_throws ArgumentError to_egs(p_iaea_bad)
    end

    # latch and ZLAST
    p_egs = arbitrary(EGSParticle{Float32})
    p_iaea = to_iaea(p_egs, z=0)
    @test p_iaea.extra_floats === (p_egs.zlast,)
    ilatch = Int32(UInt32(p_egs.latch))
    @test p_iaea.extra_ints === (ilatch, )
end

@testset "phsp_convert" begin
    src = testdatapath("photon_electron_positron.egsphsp")
    dst = IAEAPath(tempname())
    z = rand(Float32)

    phsp_convert(src, dst, z=z)
    ps_egs = PhspVector(src)
    ps_iaea = PhspVector(dst)
    ps_iaea2 = map(ps_egs) do p
        to_iaea(p,z=z)
    end

    @test ps_iaea == ps_iaea2
    orig_egs = ps_egs.header.originalcount
    orig_iaea = parse(Float32, ps_iaea.header.attributes[:ORIG_HISTORIES])
    @test orig_egs == orig_iaea
end

end#module


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

## dummy stub for https://github.com/JuliaBinaryWrappers/CompilerSupportLibraries_jll.jl

baremodule CompilerSupportLibraries_jll
using Base, Libdl, Base.BinaryPlatforms
Base.Experimental.@compiler_options compile=min optimize=0 infer=false

const PATH_list = String[]
const LIBPATH_list = String[]

export libgfortran, libstdcxx, libgomp

# These get calculated in __init__()
PATH = Ref("")
LIBPATH = Ref("")
artifact_dir = ""
libgfortran_handle = C_NULL
libgfortran_path = ""
libstdcxx_handle = C_NULL
libstdcxx_path = ""
libgomp_handle = C_NULL
libgomp_path = ""

if Sys.iswindows()
    if arch(HostPlatform()) == "x86_64"
        const libgcc_s = "libgcc_s_seh-1.dll"
    else
        const libgcc_s = "libgcc_s_sjlj-1.dll"
    end
    const libgfortran = string("libgfortran-", libgfortran_version(HostPlatform()).major, ".dll")
    const libstdcxx = "libstdc++-6.dll"
    const libgomp = "libgomp-1.dll"
elseif Sys.isapple()
    if arch(HostPlatform()) == "aarch64"
        const libgcc_s = "@rpath/libgcc_s.2.dylib"
    else
        const libgcc_s = "@rpath/libgcc_s.1.dylib"
    end
    const libgfortran = string("@rpath/", "libgfortran.", libgfortran_version(HostPlatform()).major, ".dylib")
    const libstdcxx = "@rpath/libstdc++.6.dylib"
    const libgomp = "@rpath/libgomp.1.dylib"
else
    const libgcc_s = "libgcc_s.so.1"
    const libgfortran = string("libgfortran.so.", libgfortran_version(HostPlatform()).major)
    const libstdcxx = "libstdc++.so.6"
    const libgomp = "libgomp.so.1"
end

function __init__()
    global artifact_dir = dirname(Sys.BINDIR)
    global LIBPATH[] = joinpath(Sys.BINDIR, Base.LIBDIR, "julia")
    global libgcc_s_handle = dlopen(libgcc_s)
    global libgcc_s_path = dlpath(libgcc_s_handle)
    global libgfortran_handle = dlopen(libgfortran)
    global libgfortran_path = dlpath(libgfortran_handle)
    global libstdcxx_handle = dlopen(libstdcxx)
    global libstdcxx_path = dlpath(libstdcxx_handle)
    global libgomp_handle = dlopen(libgomp)
    global libgomp_path = dlpath(libgomp_handle)
end

# JLLWrappers API compatibility shims.  Note that not all of these will really make sense.
# For instance, `find_artifact_dir()` won't actually be the artifact directory, because
# there isn't one.  It instead returns the overall Julia prefix.
is_available() = true
find_artifact_dir() = artifact_dir
dev_jll() = error("stdlib JLLs cannot be dev'ed")
best_wrapper = nothing
get_libgfortran_path() = libgfortran_path
get_libstdcxx_path() = libstdcxx_path
get_libgomp_path() = libgomp_path

end  # module CompilerSupportLibraries_jll

max_stars_repo_name: dmageeLANL/julia

max_issues_repo_name: dmageeLANL/julia

max_forks_repo_name: dmageeLANL/julia


text:
using Distributions
using PSFModels
using StableRNGs
using StaticArrays
using Test

rng = StableRNG(112358)

modname(::typeof(gaussian)) = "gaussian"
modname(::typeof(airydisk)) = "airydisk"
modname(::typeof(moffat)) = "moffat"

function test_model_interface(K)
    # test defaults
    m = K(x=0, y=0, fwhm=10)
    val_dir = @inferred K(0, 0; x=0, y=0, fwhm=10)
    val_call = @inferred m(0, 0)
    val_tup = @inferred m((0, 0))
    val_vec = @inferred m(SA[0, 0])
    @test val_dir ≈ val_call ≈ val_tup ≈ val_vec ≈ 1

    # test new position
    m = K(x=12, y=13, fwhm=10)
    val_dir = @inferred K(12, 13; x=12, y=13, fwhm=10)
    @test m(12, 13) ≈ val_dir ≈ 1

    # test diagonal fwhm
    m = K(x=0, y=0, fwhm=(10, 9))
    @test m(3, 5) ≈ K(3, 5; x=0, y=0, fwhm=(10, 9)) ≈ K(3, 5; x=0, y=0, fwhm=SA[10, 9])

    # test rotation
    m = K(x=0, y=0, fwhm=(3, 5))
    m90 = K(x=0, y=0, fwhm=(5, 3), theta=90)
    @test m(7, 8) ≈ m90(7, 8)

    mwarn = K(x=0, y=0, fwhm=10, theta=30)
    name = modname(K)
    expected_log = (:warn, "isotropic $name is not affected by non-zero rotation angle 30")
    @test_logs expected_log mwarn(0, 0)

    # test different amplitude
    m = K(x=0, y=0, amp=2, fwhm=9)
    val_dir = @inferred K(0, 0; x=0, y=0, amp=2.0, fwhm=9)
    @test m(0, 0) ≈ val_dir ≈ 2.0

    # test background
    m = K(x=0, y=0, fwhm=9, bkg=10)
    val_dir = @inferred K(0, 0; x=0, y=0, fwhm=9, bkg=10)
    @test m(0, 0) ≈ val_dir ≈ 11

    # test different type
    m = K(BigFloat; x=0, y=0, fwhm=10)
    val_dir = @inferred K(BigFloat, 0, 0; x=0, y=0, fwhm=10)
    @test val_dir isa BigFloat
    @test m(0, 0) ≈ val_dir ≈ BigFloat(1)
end

@testset "Model Interface - $K" for K in (gaussian, airydisk, moffat)
    test_model_interface(K)
end

@testset "gaussian" begin
    m = gaussian(x=0, y=0, fwhm=10)
    expected = exp(-4 * log(2) * sum(abs2, SA[1, 2]) / 100)
    @test m(1, 2) ≈ expected

    m = gaussian(x=0, y=0, fwhm=(10, 9))
    wdist = (1/10)^2 + (2/9)^2
    expected = exp(-4 * log(2) * wdist)
    @test m(1, 2) ≈ expected

    # test Normal alias
    @test normal(0, 0; x=0, y=0, fwhm=10) === gaussian(0, 0; x=0, y=0, fwhm=10)
end

@testset "airydisk" begin
    fwhm = 10
    m = airydisk(x=0, y=0, fwhm=fwhm)
    ld = fwhm / 1.028
    radius = 1.22 * ld
    # first radius is 0
    @test m(radius, 0) ≈ 0 atol=1e-5
    @test m(-radius, 0) ≈ 0 atol=1e-5
    @test m(0, radius) ≈ 0 atol=1e-5
    @test m(0, -radius) ≈ 0 atol=1e-5

    # second radius is 0
    radius2 = 2.23 * ld
    @test m(radius2, 0) ≈ 0 atol=1e-5
    @test m(-radius2, 0) ≈ 0 atol=1e-5
    @test m(0, radius2) ≈ 0 atol=1e-5
    @test m(0, -radius2) ≈ 0 atol=1e-5

    fwhm = (10, 6)
    m = airydisk(x=0, y=0, fwhm=fwhm)
    r1 = fwhm[1] * 1.18677
    r2 = fwhm[2] * 1.18677
    # first radius is 0
    @test m(r1, 0) ≈ 0 atol=1e-5
    @test m(-r1, 0) ≈ 0 atol=1e-5
    @test m(0, r2) ≈ 0 atol=1e-5
    @test m(0, -r2) ≈ 0 atol=1e-5

    # test with ratio
    mratio = airydisk(x=0, y=0, fwhm=10, ratio=sqrt(0.5))
    # test attenuation
    @test mratio(0, 0) ≈ 4
    @test mratio(radius, 0) > m(radius, 0)
    @test mratio(-radius, 0) > m(-radius, 0)
    @test mratio(0, radius) > m(0, radius)
    @test mratio(0, -radius) > m(0, -radius)

    mratio = airydisk(x=0, y=0, fwhm=(10, 6), ratio=sqrt(0.5))
    r1 = fwhm[1] * 1.18677
    r2 = fwhm[2] * 1.18677
    # test attenuation
    @test mratio(0, 0) ≈ 4
    @test mratio(r1, 0) > m(r1, 0)
    @test mratio(-r1, 0) > m(-r1, 0)
    @test mratio(0, r2) > m(0, r2)
    @test mratio(0, -r2) > m(0, -r2)
end

@testset "moffat" begin
    m = moffat(x=0, y=0, fwhm=10)
     expected = inv(1 + sum(abs2, SA[1, 2]) / 25)
    @test m(1, 2) ≈ expected

    m = moffat(x=0, y=0, fwhm=(10, 9))
    wdist = (1/5)^2 + (2/4.5)^2
    expected = inv(1 + wdist)
    @test m(1, 2) ≈ expected

    # different alpha
    m = moffat(x=0, y=0, fwhm=10, alpha=2)
    expected = inv(1 + sum(abs2, SA[1, 2]) / 25)^2
    @test m(1, 2) ≈ expected
end

include("plotting.jl")
include("fitting.jl")


text:
@testset "testing catobs & getindex operations missing values" begin
    a = BagNode(ArrayNode(rand(3,4)), [1:4], nothing)
    e = BagNode(missing, AlignedBags([0:-1]), nothing)

    x = reduce(catobs, [a, e])
    @test x.data.data == a.data.data
    @test x.bags.bags == [1:4, 0:-1]

    x = reduce(catobs, [e, a])
    @test x.data.data == a.data.data
    @test x.bags.bags == [0:-1, 1:4]

    x = reduce(catobs, [e, e])
    @test ismissing(x.data)
    @test x.bags.bags == [0:-1, 0:-1]

    x = reduce(catobs, [a, e])
    @test  isnothing(x[2].metadata)
    Mill.emptyismissing!(true)
    @test  ismissing(x[2].data)
    Mill.emptyismissing!(false)
    @test  nobs(x[2].data) == 0
    @test  x[2].bags.bags == [0:-1]

    @test  x[1].data.data == a.data.data
    @test  x[1].bags.bags == [1:4]

    @test catobs(a, e).data.data == a.data.data
    @test catobs(a, e).bags.bags == [1:4, 0:-1]

    a = BagNode(a, [1:1], missing)
    @test catobs(a, e).bags.bags == [1:1, 0:-1]
    @test catobs(e, a).bags.bags == [0:-1, 1:1]
    @test catobs(e, a).data.data.data == a.data.data.data

    @test catobs(a,e,a,e).data.data.data == catobs(a,a).data.data.data
    @test catobs(a,e,a,e).bags.bags == [1:1, 0:-1, 2:2, 0:-1]

    @test_throws ErrorException BagNode(missing, AlignedBags([1:3]), nothing)
    @test_throws ErrorException BagNode(missing, Mill.ScatteredBags([[1,2,3]]), nothing)
end

@testset "testing catobs & getindex operations missing values for weighted" begin
    a = WeightedBagNode(ArrayNode(rand(3,4)), [1:4], [1.0, 0.0, 1.0, 0.5], nothing)
    e = WeightedBagNode(missing, AlignedBags([0:-1]), [], nothing)

    x = reduce(catobs, [a, e])
    @test x.data.data == a.data.data
    @test x.bags.bags == [1:4, 0:-1]
    @test x.weights == [1.0, 0.0, 1.0, 0.5]

    x = reduce(catobs, [e, a])
    @test x.data.data == a.data.data
    @test x.bags.bags == [0:-1, 1:4]
    @test x.weights == [1.0, 0.0, 1.0, 0.5]

    x = reduce(catobs, [e, e])
    @test ismissing(x.data)
    @test x.bags.bags == [0:-1, 0:-1]
    @test x.weights == []

    x = reduce(catobs, [a, e])
    @test  isnothing(x[2].metadata)
    Mill.emptyismissing!(true)
    @test  ismissing(x[2].data)
    Mill.emptyismissing!(false)
    @test  nobs(x[2].data) == 0
    @test  x[2].bags.bags == [0:-1]
    @test x[1].weights == [1.0, 0.0, 1.0, 0.5]
    @test x[2].weights == []
    @test  x[1].data.data == a.data.data
    @test  x[1].bags.bags == [1:4]
end

@testset "testing model operations missing values" begin
    a = BagNode(ArrayNode(rand(3, 4)), [1:4], nothing)
    e = BagNode(missing, AlignedBags([0:-1]), nothing)
    m = BagModel(ArrayModel(Dense(3, 2)), SegmentedMean(2), ArrayModel(Dense(2, 2)))

    @testset "BagNode" begin
        x = reduce(catobs, [a, e])
        @test m(x).data[:, 1] ≈ m(a).data
        @test m(x).data[:, 2] ≈ m(e).data
        @test m(x).data ≈ hcat(m(a).data, m(e).data)
    end

    @testset "ProductNode" begin
        t1 = ProductNode((a, a))
        t2 = ProductNode((a, e))
        t3 = ProductNode((e, a))
        t4 = ProductNode((e, e))
        tt = [t1, t2, t3, t4]
        x  = reduce(catobs, tt)
        tm = ProductModel((m, m))
        o = tm(x).data
        for i in 1:length(tt)
            @test o[:,i] ≈ tm(tt[i]).data
        end
        for i in 1:length(tt)
            @test o[:,i] ≈ tm(x[i]).data
        end
    end
end

@testset "testing reduction with missing values" begin
    a = ProductNode((x=ArrayNode(rand(3,4)), y=ArrayNode(rand(3,4))))
    b = ProductNode((x=ArrayNode(rand(3,4)), y=ArrayNode(rand(3,4))))

    c = ProductNode((
        x=BagNode(missing, AlignedBags([0:-1]), nothing), y=ArrayNode(rand(3,1)),
    ))
    d = ProductNode((
        x=BagNode(a, [1:4], nothing), y=ArrayNode(rand(3,1)),
    ))

    ba = BagNode(a, [1:4], nothing)
    bb = BagNode(b, [1:4], nothing)
    bc = BagNode(missing, AlignedBags([0:-1]), nothing)

    bba = BagNode(c, [1:1], nothing)
    bbb = BagNode(d, [1:1], nothing)
    bbc = BagNode(missing, AlignedBags([0:-1]), nothing)

    x = reduce(catobs, [c, d])
    y = reduce(catobs, [ba, bb, bc])
    z = reduce(catobs, [bba, bbb, bbc])

    mx = reflectinmodel(x, d -> Dense(d, 2, relu))
    my = reflectinmodel(y, d -> Dense(d, 2, relu))
    mz = reflectinmodel(z, d -> Dense(d, 2, relu))

    @test mx(x).data[:, 1] ≈ mx(c).data
    @test mx(x).data[:, 2] ≈ mx(d).data

    @test my(y).data[:, 1] ≈ my(ba).data
    @test my(y).data[:, 2] ≈ my(bb).data
    @test my(y).data[:, 3] ≈ my(bc).data

    @test mz(z).data[:, 1] ≈ mz(bba).data
    @test mz(z).data[:, 2] ≈ mz(bbb).data
    @test mz(z).data[:, 3] ≈ mz(bbc).data
end


text:
function submit(cur_player::AbstractPlayer, cur_name::AbstractString)

  cur_round = cur_player.round

  cur_request = Dict(
    "score" => Dict(
      "name" => cur_name,
      "round" => Dict(
        "id" => cur_player.round.server_id,
        "password" => cur_player.round.client_password
      )
    )
  )

  did_save_score = false

  cur_response = nothing

  cur_error = nothing

  try
    cur_response = HTTP.post(
      scores_endpoint,
      json_headers,
      JSON.json(cur_request)
    )

    did_save_score = ( cur_response.status == 201 )

    did_save_score || ( cur_error = cur_response )
  catch cur_error
    any(
      cur_error_type -> isa(cur_error, cur_error_type),
      api_errors
    ) || rethrow(cur_error)
  end

  if !did_save_score
    println(cur_error)
    return
  end

  update_scores(cur_player)

  set_default_values(cur_player)

  tetris_js(
    cur_player.game.scope,
    JSString("""
      \$(".cs-scores-screen").removeClass("hidden");
      \$(".js-new-score").addClass("hidden");
      \$(".js-new-score input").val("");
    """)
  )

end


text:
# Cu_Lorentz-Drude_0.207-12.4µm (Johnson & Christy)
# wavelength (micron) Re(n) Im(n)
const cu_data = [0.20664 1.1018 1.8654
0.21094 1.1579 1.8851
0.21532 1.2145 1.8982
0.2198 1.2707 1.9041
0.22437 1.325 1.9029
0.22903 1.3761 1.8947
0.23379 1.4225 1.88
0.23865 1.4632 1.8596
0.24361 1.4969 1.8346
0.24868 1.5229 1.8062
0.25385 1.5405 1.7759
0.25912 1.5496 1.7452
0.26451 1.5501 1.7155
0.27001 1.5424 1.6882
0.27562 1.527 1.6647
0.28135 1.5045 1.6461
0.2872 1.4758 1.6335
0.29317 1.4419 1.6279
0.29926 1.4041 1.63
0.30548 1.3636 1.6405
0.31183 1.3219 1.6598
0.31832 1.2807 1.688
0.32493 1.2418 1.725
0.33169 1.2068 1.7703
0.33858 1.1776 1.8227
0.34562 1.1556 1.8809
0.35281 1.1423 1.943
0.36014 1.1386 2.0067
0.36763 1.1448 2.0693
0.37527 1.1608 2.1277
0.38307 1.1854 2.1787
0.39103 1.2163 2.219
0.39916 1.2496 2.2459
0.40746 1.2805 2.258
0.41593 1.3033 2.2558
0.42458 1.3123 2.2423
0.4334 1.3036 2.2223
0.44241 1.2752 2.2022
0.45161 1.2279 2.1883
0.461 1.1644 2.1861
0.47058 1.0892 2.1995
0.48036 1.0073 2.2301
0.49035 0.92356 2.2781
0.50054 0.84204 2.3418
0.51094 0.7656 2.4188
0.52157 0.6959 2.5066
0.53241 0.63358 2.6025
0.54348 0.57861 2.7044
0.55477 0.53052 2.8106
0.56631 0.48865 2.9199
0.57808 0.45227 3.0314
0.5901 0.42068 3.1446
0.60236 0.39326 3.259
0.61488 0.36943 3.3746
0.62767 0.34872 3.491
0.64071 0.33072 3.6084
0.65403 0.31507 3.7266
0.66763 0.30149 3.8457
0.68151 0.28971 3.9658
0.69567 0.27955 4.0869
0.71014 0.27081 4.2092
0.7249 0.26336 4.3326
0.73997 0.25706 4.4573
0.75535 0.2518 4.5834
0.77105 0.24751 4.7109
0.78708 0.24409 4.8399
0.80344 0.24149 4.9706
0.82014 0.23965 5.1029
0.83719 0.23853 5.2371
0.85459 0.23807 5.373
0.87236 0.23825 5.511
0.89049 0.23904 5.6509
0.909 0.24042 5.793
0.9279 0.24237 5.9372
0.94719 0.24487 6.0836
0.96688 0.24791 6.2324
0.98698 0.25148 6.3836
1.0075 0.25557 6.5372
1.0284 0.26019 6.6934
1.0498 0.26532 6.8523
1.0716 0.27096 7.0138
1.0939 0.27713 7.1781
1.1167 0.28382 7.3452
1.1399 0.29104 7.5152
1.1636 0.2988 7.6883
1.1878 0.30709 7.8643
1.2124 0.31595 8.0436
1.2376 0.32537 8.226
1.2634 0.33536 8.4118
1.2896 0.34595 8.6009
1.3164 0.35714 8.7934
1.3438 0.36896 8.9895
1.3717 0.38142 9.1891
1.4003 0.39454 9.3925
1.4294 0.40834 9.5995
1.4591 0.42284 9.8104
1.4894 0.43806 10.025
1.5204 0.45403 10.244
1.552 0.47077 10.467
1.5842 0.4883 10.694
1.6172 0.50666 10.925
1.6508 0.52587 11.161
1.6851 0.54596 11.4
1.7201 0.56696 11.645
1.7559 0.5889 11.894
1.7924 0.6118 12.147
1.8297 0.6357 12.405
1.8677 0.66064 12.668
1.9065 0.68664 12.936
1.9462 0.71373 13.209
1.9866 0.74195 13.486
2.0279 0.77134 13.769
2.0701 0.80191 14.057
2.1131 0.8337 14.35
2.157 0.86674 14.649
2.2019 0.90106 14.952
2.2476 0.93669 15.262
2.2944 0.97364 15.576
2.3421 1.0119 15.897
2.3907 1.0516 16.223
2.4404 1.0926 16.554
2.4912 1.1351 16.892
2.543 1.1789 17.235
2.5958 1.2241 17.584
2.6498 1.2707 17.939
2.7049 1.3186 18.3
2.7611 1.3679 18.668
2.8185 1.4185 19.041
2.8771 1.4704 19.421
2.9369 1.5235 19.808
2.9979 1.5778 20.2
3.0602 1.6332 20.6
3.1239 1.6896 21.006
3.1888 1.747 21.419
3.2551 1.8053 21.839
3.3228 1.8643 22.266
3.3918 1.9241 22.7
3.4623 1.9844 23.143
3.5343 2.0452 23.592
3.6078 2.1064 24.05
3.6828 2.1679 24.517
3.7593 2.2296 24.991
3.8375 2.2914 25.475
3.9173 2.3534 25.968
3.9987 2.4154 26.471
4.0818 2.4775 26.983
4.1667 2.5396 27.507
4.2533 2.6018 28.04
4.3417 2.6641 28.585
4.4319 2.7267 29.142
4.5241 2.7895 29.711
4.6181 2.8529 30.292
4.7141 2.9168 30.886
4.8121 2.9816 31.494
4.9122 3.0474 32.115
5.0143 3.1145 32.751
5.1185 3.1832 33.401
5.2249 3.2536 34.066
5.3335 3.3261 34.746
5.4444 3.4011 35.442
5.5576 3.4787 36.154
5.6731 3.5593 36.882
5.791 3.6432 37.626
5.9114 3.7308 38.388
6.0343 3.8223 39.166
6.1597 3.9182 39.962
6.2878 4.0186 40.775
6.4185 4.1239 41.606
6.5519 4.2344 42.456
6.6881 4.3504 43.323
6.8271 4.4722 44.209
6.9691 4.6002 45.114
7.1139 4.7346 46.038
7.2618 4.8758 46.981
7.4128 5.0239 47.944
7.5669 5.1794 48.926
7.7242 5.3425 49.928
7.8847 5.5136 50.95
8.0486 5.6929 51.993
8.2159 5.8808 53.056
8.3867 6.0776 54.14
8.5611 6.2835 55.245
8.739 6.499 56.372
8.9207 6.7244 57.52
9.1061 6.96 58.691
9.2954 7.2061 59.883
9.4887 7.4632 61.098
9.6859 7.7316 62.335
9.8873 8.0117 63.596
10.093 8.3039 64.879
10.303 8.6086 66.186
10.517 8.9261 67.517
10.735 9.257 68.872
10.959 9.6017 70.251
11.186 9.9606 71.655
11.419 10.334 73.083
11.656 10.723 74.536
11.899 11.127 76.015
12.146 11.548 77.519
12.398 11.985 79.048]


text:
using BasisFunctions, FrameFun, Test

@testset  "AugmentationPlatform" begin
    plat1 = AugmentationPlatform(platform(Fourier(10)),Fourier(5))
    plat2 = ONB_plus_K(platform(Fourier(10)),sqrt, 5)
    plat3 = AugmentationPlatform(platform(ChebyshevT(10)),ChebyshevT(5))
    plat4 = ONB_plus_K(platform(ChebyshevT(10)),sqrt, 5)
    plat5 = AugmentationPlatform(platform(ChebyshevT(10)),Fourier(5))
    plat6 = AugmentationPlatform(platform(Fourier(10)),ChebyshevT(5))

    for plat in (plat1,plat2,plat3,plat4,plat5,plat6)
        @test SamplingStyle(plat) == OversamplingStyle()
        @test DictionaryStyle(plat) == FrameStyle()
        @test SolverStyle(plat,OversamplingStyle()) == AZStyle()
    end

    for plat in (plat1,plat2,plat6)
        @test measure(plat) == FourierWeight()
    end
    for plat in (plat3,plat4,plat5)
        @test measure(plat) == ChebyshevWeight()
    end

    dict1 = dictionary(plat1,20)
    dict2 = dictionary(plat2,20)
    dict3 = dictionary(plat3,20)
    dict4 = dictionary(plat4,20)
    dict5 = dictionary(plat5,20)
    dict6 = dictionary(plat6,20)

    for dict in (dict1,dict2,dict3,dict4,dict5,dict6)
        @test dict isa MultiDict
        @test dimensions(dict) == [20,5]
    end

    # @test_throws Exception dualdictionary(plat1,20,measure(plat1))
    # @test_throws Exception dualdictionary(plat2,20,measure(plat2))
    # @test_throws Exception dualdictionary(plat3,20,measure(plat3))
    # @test_throws Exception dualdictionary(plat4,20,measure(plat4))
    # @test_throws Exception dualdictionary(plat5,20,measure(plat5))
    # @test_throws Exception dualdictionary(plat6,20,measure(plat6))

    # TODO Product platform tests

end

@testset "WeightedSumPlatforms" begin
    plat = WeightedSumPlatform(platform(Fourier(10)),x->1., sqrt)
    @test SamplingStyle(plat) == OversamplingStyle()
    @test DictionaryStyle(plat) == FrameStyle()
    @test SolverStyle(plat,OversamplingStyle()) == AZStyle()
    @test measure(plat) == FourierWeight()
    dict = dictionary(plat,(10,10))
    @test dict isa MultiDict
    @test dimensions(dict) == [10,10]
    ddict =  dualdictionary(plat,(10,10),FourierWeight())
    @test ddict isa MultiDict
    @test dimensions(ddict) == [10,10]
    @test param_first(plat) == (10,10)

end


text:
using Test

using MathOptInterface
const MOI = MathOptInterface
using JuMP
using ParameterJuMP

using GLPK
factory = GLPK.Optimizer

include("tests.jl")

@testset "ParameterJuMP tests" begin
    test0(factory)
    test1(factory)
    test2(factory)
    test3(factory)
    test4(factory)
    test5(factory)
    test6(factory)
    test7(factory)
    test8(factory)
    test9(factory)
    test10(factory)
    test11(factory)
    test12(factory)
    test13(factory)
    test14(factory)
    test15(factory)
    test16(factory)
end
;


text:
using LazySets, MathematicalPredicates, SparseArrays

function building_specification()
    # initial set:
    # - xᵢ ∈ [0.0002, 0.00025] for i ∈ [1, 10]
    # - xᵢ ∈ [-0.0001, 0.0001] for i = 25
    # - xᵢ = 0 for all other i
    X0 = Hyperrectangle(low=[fill(0.0002, 10); zeros(14); -0.0001; zeros(23)],
                        high=[fill(0.00025, 10); zeros(14); 0.0001; zeros(23)])

    # safety property: x25 ≤ 6e-3
    property = is_contained_in(HalfSpace(sparsevec([25], [1.0], 48), 6e-3))

    # time horizon: 20 time units
    time_horizon = 20.0

    # specification
    O = Dict(:T => time_horizon, :property => property)

    return X0, O
end

max_stars_repo_name: JuliaReach/ReachabilityBenchmarks

max_issues_repo_name: JuliaReach/ReachabilityBenchmarks

max_forks_repo_name: JuliaReach/ReachabilityBenchmarks


text:
"""
    assemble_stiffness_matrix(elements, nodes[, D])

Assemble 3D stiffness matrix using P1 finite elements.
A diffusion tensor `D` may be provided.

This function is based on the Matlab function `stiffness_matrixP1_3D.m` from the Matlab
nodal matrix assembly toolbox by Jan Valdman and Talal Rahman.

https://uk.mathworks.com/matlabcentral/fileexchange/27826-fast-fem-assembly-nodal-elements

Talal Rahman and Jan Valdman: Fast MATLAB assembly of FEM matrices in 2D and 3D: nodal
elements, Applied Mathematics and Computation 219, 7151–7158 (2013).
"""
function assemble_stiffness_matrix(elements, nodes, D = I(3))
    nelement = size(elements, 1)
    nodes = permutedims(nodes[elements, :], [3, 2, 1])

    integration_point = [1 // 4; 1 // 4; 1 // 4]
    ∇φ, detjac, _ = compute_∇φ(nodes, integration_point)

    volumes = abs.(detjac[:]) / 6
    ∇φ = ∇φ[:, :, 1, :]

    z = zeros(4, 4, nelement)
    for ielement = 1:nelement
        z[:, :, ielement] = volumes[ielement] * ∇φ[:, :, ielement]' * D * ∇φ[:, :, ielement]
    end
    y = reshape(repeat(elements, 1, 4)', 4, 4, nelement)
    x = permutedims(y, [2, 1, 3])

    # Enforce symmetry
    sym(x) = (x + x') / 2

    sym(sparse(x[:], y[:], z[:]))
end

"""
    compute_∇φ(nodes, points)

Compute the gradients of the finite element basis functions φ.
"""
function compute_∇φ(nodes, points)
    nelement = size(nodes, 3)
    npoint = size(points, 2)
    dshape = shapeder(points)
    ∇φ = zeros(3, 4, npoint, nelement)
    detjac = zeros(npoint, nelement)
    jac = zeros(3, 3, npoint, nelement)
    for ipoint = 1:npoint
        for ielement = 1:nelement
            j = dshape[:, :, ipoint] * nodes[:, :, ielement]'
            jinv = inv(j)
            detjac[ipoint, ielement] = abs(det(j))
            jac[:, :, ipoint, ielement] = j
            ∇φ[:, :, ipoint, ielement] = jinv * dshape[:, :, ipoint]
        end
    end
    ∇φ, detjac, jac
end

"""
    shapeder(points)

Derivative of shape functions with respect to reference coordinates (`ξ₁`, `ξ₂`, `ξ₃`).
"""
function shapeder(points)
    npoint = size(points, 2)
    dshape = [
        1.0 0.0 0.0 -1.0
        0.0 1.0 0.0 -1.0
        0.0 0.0 1.0 -1.0
    ]
    dshape = repeat(dshape, 1, 1, npoint)
end


text:
# This file contains the Port tool for connecting the tools of Causal

"""
    $TYPEDEF

Abstract type of [`Outport`](@ref) and [`Inport`](@ref). See also: [`Outport`](@ref), [`Inport`](@ref).
"""
abstract type AbstractPort{P} <: AbstractVector{P} end

"""
    $TYPEDEF

Constructs an `Outport` with `numpins` [`Outpin`](@ref).

# Fields

    $TYPEDFIELDS

!!! warning Element type of an `Outport` must be `Outpin`. See also [`Outpin`](@ref)

# Example
```jldoctest
julia> Outport{Int}(2)
2-element Outport{Outpin{Int64}}:
 Outpin(eltype:Int64, isbound:false)
 Outpin(eltype:Int64, isbound:false)

julia> Outport()
1-element Outport{Outpin{Float64}}:
 Outpin(eltype:Float64, isbound:false)
```
"""
struct Outport{P} <: AbstractPort{P}
    "Pins of the output port"
    pins::Vector{P}
    "Unique identifier"
    id::UUID
    Outport(pins::AbstractVector{P}) where {T, P<:Outpin{T}} = new{P}(pins, uuid4())
end
Outport(pin::Outpin) = Outport([pin])
Outport{T}(numpins::Int=1) where T = Outport([Outpin{T}() for i = 1 : numpins])
Outport(numpins::Int=1) = Outport{Float64}(numpins)

show(io::IO, outport::Outport) = print(io, "Outport(numpins:$(length(outport)), eltype:$(eltype(outport)))")

"""
    $TYPEDEF

Constructs an `Inport` with `numpins` [`Inpin`](@ref).

# Fields

    $TYPEDFIELDS

!!! warning Element type of an `Inport` must be `Inpin`. See also [`Inpin`](@ref)

# Example
```jldoctest
julia> Inport{Int}(2)
2-element Inport{Inpin{Int64}}:
 Inpin(eltype:Int64, isbound:false)
 Inpin(eltype:Int64, isbound:false)

julia> Inport()
1-element Inport{Inpin{Float64}}:
 Inpin(eltype:Float64, isbound:false)
```
"""
struct Inport{P} <: AbstractPort{P}
    "Pins of the port"
    pins::Vector{P}
    Inport(pins::AbstractVector{P}) where {T, P<:Inpin{T}} = new{P}(pins)
end
Inport(pin::Inpin) = Inport([pin])
Inport{T}(numpins::Int=1) where T = Inport([Inpin{T}() for i = 1 : numpins])
Inport(numpins::Int=1) = Inport{Float64}(numpins)

show(io::IO, inport::Inport) = print(io, "Inport(numpins:$(length(inport)), eltype:$(eltype(inport)))")

"""
    $SIGNATURES

Returns the data type of `port`.
"""
datatype(port::AbstractPort{<:AbstractPin{T}}) where {T} = T

##### AbstractVector interface
"""
    $SIGNATURES

Retruns size of `port`.
"""
size(port::AbstractPort) = size(port.pins)

"""
    $SIGNATURES

Returns elements from `port` at index `idx`. Same as `port[idx]`.

# Example
```jldoctest
julia> op = Outport(3)
3-element Outport{Outpin{Float64}}:
 Outpin(eltype:Float64, isbound:false)
 Outpin(eltype:Float64, isbound:false)
 Outpin(eltype:Float64, isbound:false)

julia> op[1]
Outpin(eltype:Float64, isbound:false)

julia> op[end]
Outpin(eltype:Float64, isbound:false)

julia> op[:]
3-element Array{Outpin{Float64},1}:
 Outpin(eltype:Float64, isbound:false)
 Outpin(eltype:Float64, isbound:false)
 Outpin(eltype:Float64, isbound:false)
```
"""
getindex(port::AbstractPort, idx::Vararg{Int, N}) where {N} = port.pins[idx...]

"""
    $SIGNATURES

Sets `item` to `port` at index `idx`. Same as `port[idx] = item`.

# Example
```jldoctest
julia> op = Outport(3)
3-element Outport{Outpin{Float64}}:
 Outpin(eltype:Float64, isbound:false)
 Outpin(eltype:Float64, isbound:false)
 Outpin(eltype:Float64, isbound:false)

julia> op[1] = Outpin()
Outpin(eltype:Float64, isbound:false)

julia> op[end] = Outpin()
Outpin(eltype:Float64, isbound:false)

julia> op[1:2] = [Outpin(), Outpin()]
2-element Array{Outpin{Float64},1}:
 Outpin(eltype:Float64, isbound:false)
 Outpin(eltype:Float64, isbound:false)
```
"""
setindex!(port::AbstractPort, item, idx::Vararg{Int, N}) where {N} = port.pins[idx...] = item

##### Reading from and writing into from buses
"""
    $SIGNATURES

Takes an element from `inport`. Each link of the `inport` is a read and a vector containing the results is returned.

!!! warning
    The `inport` must be readable to be read. That is, there must be a runnable tasks bound to links of the `inport`
    that writes data to `inport`.

# Example
```jldoctest
julia> op, ip = Outport(), Inport()
(Outport(numpins:1, eltype:Outpin{Float64}), Inport(numpins:1, eltype:Inpin{Float64}))

julia> ls = connect!(op, ip)
1-element Array{Link{Float64},1}:
 Link(state:open, eltype:Float64, isreadable:false, iswritable:false)

julia> t = @async for val in 1 : 5
       put!(op, [val])
       end;

julia> take!(ip)
1-element Array{Float64,1}:
 1.0

julia> take!(ip)
1-element Array{Float64,1}:
 2.0
```
"""
take!(inport::Inport) = take!.(inport[:])

"""
    $SIGNATURES

Puts `vals` to `outport`. Each item in `vals` is putted to the `links` of the `outport`.

!!! warning
    The `outport` must be writable to be read. That is, there must be a runnable tasks bound to links of the
    `outport` that reads data from `outport`.

# Example
```jldoctest
julia> op, ip = Outport(), Inport()
(Outport(numpins:1, eltype:Outpin{Float64}), Inport(numpins:1, eltype:Inpin{Float64}))

julia> ls = connect!(op, ip)
1-element Array{Link{Float64},1}:
 Link(state:open, eltype:Float64, isreadable:false, iswritable:false)

julia> t = @async while true
       val = take!(ip)
       all(val .=== NaN) && break
       println("Took " * string(val))
       end;

julia> put!(op, [1.])
Took [1.0]
1-element Array{Float64,1}:
 1.0

julia> put!(op, [NaN])
1-element Array{Float64,1}:
 NaN
```
"""
function put!(outport::Outport, vals)
    put!.(outport[:], vals)
    vals
end

##### Interconnection of busses.
"""
    $SIGNATURES

Returns a new port that is similar to `port` with the same element type. The number of links in the new port is `nlinks` and
data buffer length is `ln`.
"""
similar(outport::Outport{P}, numpins::Int=length(outport)) where {P<:Outpin{T}} where {T} = Outport{T}(numpins)
similar(inport::Inport{P}, numpins::Int=length(inport)) where {P<:Inpin{T}} where {T} = Inport{T}(numpins)


text:
## mfccs.jl
## (c) 2013--2014 David A. van Leeuwen

## Recoded from / inspired by melfcc from Dan Ellis's rastamat package.

using SpecialFunctions

## This function accepts a vector of sample values, below we will generalize to arrays,
## i.e., multichannel data
## Recoded from rastamat's "melfcc.m" (c) Dan Ellis.
## Defaults here are HTK parameters, this is contrary to melfcc
function mfcc{T<:AbstractFloat}(x::Vector{T}, sr::Real=16000.0; wintime=0.025, steptime=0.01, numcep=13,
              lifterexp=-22, sumpower=false, preemph=0.97, dither=false, minfreq=0.0, maxfreq=sr/2,
              nbands=20, bwidth=1.0, dcttype=3, fbtype=:htkmel, usecmp=false, modelorder=0)
    if (preemph != 0)
        x = filt(TFFilter([1., -preemph], [1.]), x)
    end
    pspec = powspec(x, sr, wintime=wintime, steptime=steptime, dither=dither)
    aspec = audspec(pspec, sr, nfilts=nbands, fbtype=fbtype, minfreq=minfreq, maxfreq=maxfreq, sumpower=sumpower, bwidth=bwidth)
    if usecmp
        #  PLP-like weighting/compression
        aspec = postaud(aspec, maxfreq, fbtype)
    end
    if modelorder > 0
        if dcttype != 1
            error("Sorry, modelorder>0 and dcttype ≠ 1 is not implemented")
        end
        # LPC analysis
        lpcas = dolpc(aspec, modelorder)
        # cepstra
        cepstra = lpc2cep(lpcas, numcep)
    else
        cepstra = spec2cep(aspec, numcep, dcttype)
    end
    cepstra = lifter(cepstra, lifterexp)'
    meta = Dict("sr" => sr, "wintime" => wintime, "steptime" => steptime, "numcep" => numcep,
            "lifterexp" => lifterexp, "sumpower" => sumpower, "preemph" => preemph,
            "dither" => dither, "minfreq" => minfreq, "maxfreq" => maxfreq, "nbands" => nbands,
            "bwidth" => bwidth, "dcttype" => dcttype, "fbtype" => fbtype,
            "usecmp" => usecmp, "modelorder" => modelorder)
    return (cepstra, pspec', meta)
end
mfcc{T<:AbstractFloat}(x::Array{T}, sr::Real=16000.0; args...) = @parallel (tuple) for i=1:size(x)[2] mfcc(x[:,i], sr; args...) end

## default feature configurations, :rasta, :htk, :spkid_toolkit, :wbspeaker
## With optional extra agrs... you can specify more options
function mfcc{T<:AbstractFloat}(x::Vector{T}, sr::AbstractFloat, defaults::Symbol; args...)
    if defaults == :rasta
        mfcc(x, sr; lifterexp=0.6, sumpower=true, nbands=40, dcttype=2, fbtype=:mel, args...)
    elseif defaults ∈ [:spkid_toolkit, :nbspeaker]
        mfcc(x, sr; lifterexp=0.6, sumpower=true, nbands=30, dcttype=2, fbtype=:mel, minfreq=130., maxfreq=3900., numcep=20, args...)
    elseif defaults == :wbspeaker
        mfcc(x, sr; lifterexp=0.6, sumpower=true, nbands=63, dcttype=2, fbtype=:mel, minfreq=62.5, maxfreq=7937.5, numcep=20, args...)
    elseif defaults == :htk
        mfcc(x, sr; args...)
    else
        error("Unknown set of defaults ", defaults)
    end
end

## our features run down with time, this is essential for the use of DSP.filt()
function deltas{T<:AbstractFloat}(x::Matrix{T}, w::Int=9)
    (nr, nc) = size(x)
    if nr==0 || w <= 1
        return x
    end
    hlen = floor(Int, w/2)
    w = 2hlen+1                 # make w odd
    win = collect(convert(T, hlen):-1:-hlen)
    x1 = reshape(x[1,:], 1, size(x,2)) ## julia v0.4 v0.5 changeover
    xend = reshape(x[end,:], 1, size(x,2))
    xx = vcat(repmat(x1, hlen, 1), x, repmat(xend, hlen, 1)) ## take care of boundaries
    norm = 3 / (hlen * w * (hlen+1))
    return norm * filt(TFFilter(win, [1.]), xx)[2hlen+(1:nr),:]
end

import Base.Sort.sortperm
sortperm(a::Array, dim::Int) = mapslices(sortperm, a, dim)

function warpstats{T<:Real}(x::Matrix{T}, w::Int=399)
    nx, nfea = size(x)
    wl = min(w, nx)
    hw = (wl+1) / 2
    erfinvtab = √2 * erfinv.(collect(1:wl) / hw .- 1)
    rank = similar(x, Int)
    if nx < w
        index = sortperm(x, 1)
        for j in 1:nfea
            rank[index[:,j], j] = collect(1:nx)
        end
    else
        for j in 1:nfea            # operations in outer loop over columns, better for memory cache
            for i in 1:nx
                s = max(1, i - round(Int, hw) + 1)
                e = s + w - 1
                if e > nx
                    d = e - nx
                    e -= d
                    s -= d
                end
                r = 1
                for k in s:e
                    r += x[i, j] > x[k, j]
                end
                rank[i, j] = r
            end
        end
    end
    return rank, erfinvtab
end

function warp{T<:Real}(x::Matrix{T}, w::Int=399)
    rank, erfinvtab = warpstats(x, w)
    return erfinvtab[rank]
end
warp{T<:Real}(x::Vector{T}, w::Int=399) = warp(x'', w)

function WarpedArray(x::Matrix, w::Int)
    rank, erfinvtab = warpstats(x, w)
    WarpedArray(rank, erfinvtab)
end

## mean and variance normalization
znorm(x::Array, dim::Int=1) = znorm!(copy(x), dim)
znorm!(x::Array, dim::Int=1) = broadcast!(/, x, broadcast!(-, x, x, mean(x, dim)), std(x, dim))

## short-term mean and variance normalization
function stmvn(x::Vector, w::Int=399)
    y = similar(x)
    hw = w ÷ 2 ## effectively makes `w` odd...
    nx = length(x)
    nx ≤ 1 && return x
    ## initialize sum x and sum x^2
    sx = (hw+1) * x[1]
    sxx = (hw+1) * x[1] * x[1]
    for i = 2:min(nx, hw+1)
        sx += x[i]
        sxx += x[i] * x[i]
    end
    if hw+1 > nx
        sx += (hw+1 - nx) * x[nx]
        sxx += (hw+1 - nx) * x[nx] * x[nx]
    end
    for i = 1:nx
        μ = sx / w
        σ = sqrt((sxx - μ * sx) / (w - 1))
        y[i] = (x[i] - μ) / σ
        mi = max(i - hw, 1)
        ma = min(i + hw + 1, nx)
        sx += x[ma] - x[mi]
        sxx += x[ma] * x[ma] - x[mi] * x[mi]
    end
    return y
end

stmvn(m::Matrix, w::Int=399, dim::Int=1) = mapslices(x->stmvn(x, w), m, dim)

## Shifted Delta Coefficients by copying
function sdc{T<:AbstractFloat}(x::Matrix{T}, n::Int=7, d::Int=1, p::Int=3, k::Int=7)
    nx, nfea = size(x)
    n <= nfea || error("Not enough features for n")
    hnfill = (k-1) * p / 2
    nfill1, nfill2 = floor(Int, hnfill), ceil(Int, hnfill)
    xx = vcat(zeros(eltype(x), nfill1, n), deltas(x[:,1:n], 2d+1), zeros(eltype(x), nfill2, n))
    y = zeros(eltype(x), nx, n*k)
    for (dt,offset) = zip(0:p:p*k-1, 0:n:n*k-1)
        y[:,offset+(1:n)] = xx[dt+(1:nx),:]
    end
    return y
end


text:
"""
    OpticalFlowAlgorithm

An optical flow algorithm with given parameters.
"""
abstract type OpticalFlowAlgorithm end

"""
    Farneback(Args...)

A method for dense optical flow estimation developed by Gunnar Farneback. It
computes the optical flow for all the points in the frame using the polynomial
representation of the images. The idea of polynomial expansion is to approximate
the neighbourhood of a point in a 2D function with a polynomial. Displacement
fields are estimated from the polynomial coefficients depending on how the
polynomial transforms under translation.

# Options
Various options for the fields of this type are described in more detail
below.

## Choices for `iterations`

Number of iterations the displacement estimation algorithm is run at each point.
If left unspecified a default value of seven iterations is assumed.

## Choices for `estimation_window`

Determines the neighbourhood size over which information will be intergrated
when determining the displacement of a pixel. The total size equals
`2*estimation_window + 1`. If left unspecified a default value of eleven is
assumed.

## Choices for `σ_estimation_window`

Standard deviation of a Gaussian weighting filter used to weigh the contribution
of a pixel's neighbourhood when determining the displacement of a pixel.
If left unspecified a default value of nine is assumed.

## Choices for `expansion_window`

Determines the size of the pixel neighbourhood used to find polynomial expansion
for each pixel; larger values mean that the image will be approximated with
smoother surfaces, yielding more robust algorithm and more blurred motion field.
The total size equals `2*expansion_window + 1`. If left unspecified a default
value of six is assumed.

## Choices for `σ_expansion_window`

Standard deviation of the Gaussian that is used to smooth the image for the purpose
of approximating it with a polynomial expansion. If left unspecified a default
value of five is assumed.

# References

1. Farnebäck G. (2003) Two-Frame Motion Estimation Based on Polynomial Expansion. In: Bigun J., Gustavsson T. (eds) Image Analysis. SCIA 2003. Lecture Notes in Computer Science, vol 2749. Springer, Berlin, Heidelberg
2. Farnebäck, G.: Polynomial Expansion for Orientation and Motion Estimation. PhD thesis, Linköping University, Sweden, SE-581 83 Linköping, Sweden (2002) Dissertation No 790, ISBN 91-7373-475-6.

"""
struct Farneback{F <: Float64, I <: Int} <: OpticalFlowAlgorithm
    iterations::I
    estimation_window::I
    σ_estimation_window::F
    expansion_window::I
    σ_expansion_window::F
end

Farneback(iterations::Int; estimation_window::Int = 11, σ_estimation_window::Real = 9.0, expansion_window::Int =  6, σ_expansion_window::Real = 5.0) = Farneback(iterations, estimation_window, σ_estimation_window, expansion_window, σ_expansion_window)

"""
    LucasKanade(Args...)

A differential method for optical flow estimation developed by Bruce D. Lucas
and Takeo Kanade. It assumes that the flow is essentially constant in a local
neighbourhood of the pixel under consideration, and solves the basic optical flow
equations for all the pixels in that neighbourhood by the least squares criterion.

# Options
Various options for the fields of this type are described in more detail
below.

## Choices for `iterations`

The termination criteria of the iterative search algorithm, that is, the number of
iterations.

## Choices for `window_size`

Size of the search window at each pyramid level; the total size of the
window used is 2*window_size + 1. If left unspecified, a default value of
eleven is assumed.

## Choices for `pyramid_levels`

0-based maximal pyramid level number; if set to 0, pyramids are not used
(single level), if set to 1, two levels are used, and so on. IF left unspecified,
a default value of four is assumed.

## Choices for `eigenvalue_threshold`

The algorithm calculates the minimum eigenvalue of a (2 x 2) normal matrix of
optical flow equations, divided by number of pixels in a window; if this value
is less than `eigenvalue_threshold`, then a corresponding feature is filtered
out and its flow is not processed (default value is 1e-4).

## `ϵ` termination criteria

Minimum required change in displacement at which the iterative algorithm continues to work.
Default value is `1e-2`.

## References

1. B. D. Lucas, & Kanade. "An Interative Image Registration Technique with an Application to Stereo Vision," DARPA Image Understanding Workshop, pp 121-130, 1981.
2. J.-Y. Bouguet, “Pyramidal implementation of the afﬁne lucas-kanade feature tracker description of the algorithm,” Intel Corporation, vol. 5,no. 1-10, p. 4, 2001.
"""
struct LucasKanade <: OpticalFlowAlgorithm
    iterations::Int
    window_size::Int
    pyramid_levels::Int
    eigenvalue_threshold::Float64
    ϵ::Float64
end

LucasKanade(
    iterations::Integer = 20; window_size::Integer = 11, pyramid_levels::Integer = 4,
    eigenvalue_threshold::Real = 1e-4, ϵ::Real = 1e-2,
) = LucasKanade(iterations, window_size,  pyramid_levels, eigenvalue_threshold, ϵ)

"""
    flow = optical_flow(source, target, Farneback(Args...))
    flow = optical_flow(source, target, displacement, Farneback(Args...))

Returns the dense optical flow from the `source` to the `target` image using the `Farneback` algorithm.

# Details

The `source` and `target` images must be `Gray` types.

The `displacement` argument allows you to specify an initial
guess for the optical flow and must be of type `Array{SVector{2, Float64}, 2}`.
The elements of `displacement` should represent the flow required to map the
`(row, column)` of each pixel in the `source` image to the `target` image.

"""
function optical_flow(source::AbstractArray{T, 2}, target::AbstractArray{T,2}, algorithm::Farneback) where T <: Gray
    # Sanity checks.
    @assert size.(axes(source)) == size.(axes(target)) "Images must have the same size"
    optflow(source, target, algorithm)
end

function optical_flow(source::AbstractArray{T, 2}, target::AbstractArray{T,2}, displacement::Array{SVector{2, Float64}, 2}, algorithm::Farneback) where T <: Gray
    # Sanity checks.
    @assert size.(axes(source)) == size.(axes(target)) "Images must have the same size"
    @assert size.(axes(source)) == size.(axes(displacement)) "Optical flow field must match image size"
    optflow!(source, target, copy(displacement), algorithm)
end

"""
    flow, indicator  = optical_flow(source, target, points, LucasKanade(Args...))
    flow, indicator  = optical_flow(source, target, points, displacement, LucasKanade(Args...))

Returns the optical flow from the `source` to the `target` image for the specified `points` using the `LucasKanade` algorithm.

# Details

The `source` and `target` images must be `Gray` types.

The `points` argument is of type `Array{SVector{2, Float64}, 2}` and represents
a set of keypoints in the `source` image for which the optical flow is to be
computed. The coordinates  of a `point` represent the  `(row, column)` of a
pixel in the image. The `displacement` argument allows you to specify an initial
guess for the optical flow.

The function returns `flow` of type `Array{SVector{2, Float64}, 2}` which
matches the length of `points` and represents the displacement needed to map
a point in the `source` image to the `target` image.

The `indicator` is a vector of boolean values (one for each point in `points`)
which signals whether a particular point was successfully tracked or not.

In order to use the `flow` to index the corresponding point in the `target` image
you first need to round it to the nearest integer, and check that it falls within
the bounds of the `target` image dimensions.

"""
function optical_flow(source::AbstractArray{T, 2}, target::AbstractArray{T,2}, points::Array{SVector{2, Float64}, 1}, algorithm::LucasKanade) where T <: Gray
    # Sanity checks.
    @assert size.(axes(source)) == size.(axes(target)) "Images must have the same size"

    optflow(source, target, points, algorithm)
end

function optical_flow(source::AbstractArray{T, 2}, target::AbstractArray{T,2}, points::Array{SVector{2, Float64}, 1},  displacement::Array{SVector{2, Float64}, 1}, algorithm::LucasKanade) where T <: Gray
    # sanity checks
    @assert size.(axes(source)) == size.(axes(target)) "Images must have the same size"
    @assert size.(axes(points)) == size.(axes(displacement)) "Vector of points and vector of displacement must have the same size"
    optflow!(source, target, points, copy(displacement), algorithm)
end

#------------------
# IMPLEMENTATIONS
#------------------

include("lucas_kanade.jl")
include("farneback.jl")


text:
@testset "HTTP.Form for multipart/form-data" begin
    headers = Dict("User-Agent" => "HTTP.jl")
    body = HTTP.Form(Dict())
    uri = "https://httpbin.org/post"
    @testset "Setting of Content-Type" begin
        for r in (HTTP.request("POST", uri, headers, body), HTTP.post(uri, headers, body))
            @test r.status == 200
            json = JSON.parse(IOBuffer(HTTP.payload(r)))
            @test startswith(json["headers"]["Content-Type"], "multipart/form-data; boundary=")
        end
    end
    @testset "Deprecation of HTTP.post without header for body::Form" begin
        proj = normpath(joinpath(pathof(HTTP), "..", "..", "Project.toml"))
        # Extract version = "(...)" from Project.toml
        vers = VersionNumber(match(r"^version\s*=\s*\"(.*?)\"$"m, read(proj, String)).captures[1])
        if vers.minor == 9 && vers.major == 0 # Keep deprecation around for at least the 0.9 release series
            depwarn_flag = Base.JLOptions().depwarn
            if depwarn_flag == 0 # silent
                @test @test_logs HTTP.post(uri, body).status == 200
            elseif depwarn_flag == 1 # warning
                @test @test_logs (:warn, r"deprecated") HTTP.post(uri, body).status == 200
            else # depwarn_flag == 2 # error
                @test_throws ErrorException HTTP.post(uri, body)
            end
        else # Next breaking release
            try
                HTTP.post(uri, body)
            catch e
                if !(e isa MethodError)
                    @warn "Deprecation for HTTP.post(uri, body) should be removed."
                end
            end
        end
    end
    @testset "HTTP.Multipart ensure show() works correctly" begin
        # testing that there is no error in printing when nothing is set for filename
        str = sprint(show, (HTTP.Multipart(nothing, IOBuffer("some data"), "plain/text", "", "testname")))
        @test findfirst("contenttype=\"plain/text\"", str) !== nothing
    end
    @testset "HTTP.Multipart test constructor" begin
        @test_nowarn HTTP.Multipart(nothing, IOBuffer("some data"), "plain/text", "", "testname")
        @test_throws MethodError HTTP.Multipart(nothing, "some data", "plain/text", "", "testname")
    end

    @testset "Boundary" begin
        @test HTTP.Form(Dict()) isa HTTP.Form
        @test HTTP.Form(Dict(); boundary="a") isa HTTP.Form
        @test HTTP.Form(Dict(); boundary=" Aa1'()+,-.:=?") isa HTTP.Form
        @test HTTP.Form(Dict(); boundary='a'^70) isa HTTP.Form
        @test_throws ArgumentError HTTP.Form(Dict(); boundary="")
        @test_throws ArgumentError HTTP.Form(Dict(); boundary='a'^71)
        @test_throws ArgumentError HTTP.Form(Dict(); boundary="a ")
    end
end


text:
using AcuteML, BenchmarkTools

# Type Definition

@aml mutable struct Person "person", check_course
    age::UInt64, "~"
    field, "study-field"
    GPA::Float64 = 4.5, "~", GPAcheck
    courses::Vector{String}, "taken-courses"
    id::Int64, att"~"
end

@aml mutable struct University
    doc"university"
    name, att"university-name"
    people::Vector{Person}, "person"
end

# Value Checking Functions
GPAcheck(x) = x <= 4.5 && x >= 0

function check_course(age, field, GPA, courses, id)

    if field == "Mechanical Engineering"
        relevant = ["Artificial Intelligence", "Robotics", "Machine Design"]
    elseif field == "Computer Engineering"
        relevant = ["Julia", "Algorithms"]
    else
        error("study field is not known")
    end

    return any(in.(courses, Ref(relevant)))
end

################################################################
#Benchmark
function creation()
    P1 = Person(age = 24, field = "Mechanical Engineering", courses = ["Artificial Intelligence", "Robotics"], id = 1)

    P2 = Person(age = 18, field = "Computer Engineering", GPA = 4, courses = ["Julia"], id = 2)

    U = University(name = "Julia University", people = [P1, P2])

    return P1, P2, U
end

function mutability(P1, P2, U)
    P1.age = UInt(22)
    # P1.courses = ["Artificial Intelligence", "Robotics", "Machine Design"]

    P2.GPA = 4.2 # mutability support
    U.name = "MIT"
    return P1, P2, U
end

function extraction()
    xml = parsexml("""
      <?xml version="1.0" encoding="UTF-8"?>
      <university university-name="MIT">
        <person id="1">
          <age>22</age>
          <study-field>Mechanical Engineering</study-field>
          <GPA>4.5</GPA>
          <taken-courses>Artificial Intelligence</taken-courses>    <taken-courses>Robotics</taken-courses>
        </person>
        <person id="2">
          <age>18</age>
          <study-field>Computer Engineering</study-field>
          <GPA>4.2</GPA>
          <taken-courses>Julia</taken-courses>
        </person>
      </university>
      """)
    U = University(xml)
end

@btime creation();

P1, P2, U = creation();

@btime mutability($P1, $P2, $U);

@btime extraction();


text:
"""

    HASLER_AND_SUSSTRUNK_M3 <: NoReferenceIQI
    M =  hasler_and_susstrunk_m3(img)

Calculates the colorfulness of a natural image using method M3 from
[1]. As a guide to interpretation of results, the authors suggest:

| Attribute           | M3 |
|:--------------------|----|
| Not colorful        |  0 |
| slightly colorful   | 15 |
| moderately colorful | 33 |
| averagely colorful  | 45 |
| quite colorful      | 59 |
| highly colorful     | 82 |
| extremely colorful  |109 |

[1] Hasler, D. and Süsstrunk, S.E., 2003, June. Measuring colorfulness
in natural images. In Human vision and electronic imaging VIII
(Vol. 5007, pp. 87-96). International Society for Optics and
Photonics.

"""
struct HASLER_AND_SUSSTRUNK_M3 <: NoReferenceIQI end

# api
(iqi::HASLER_AND_SUSSTRUNK_M3)(img) = hasler_and_susstrunk_m3(img)

@doc (@doc HASLER_AND_SUSSTRUNK_M3)
function hasler_and_susstrunk_m3(img::AbstractArray{<:AbstractRGB})

    R = 255 .* float(red.(img))
    G = 255 .* float(green.(img))
    B = 255 .* float(blue.(img))

    rg = R .- G
    μrg, σrg = mean(rg), std(rg)

    yb = 0.5 .* (R .+ G) .- B
    μyb,  σyb = mean(yb), std(yb)

    μrgyb = sqrt(μrg^2 + μyb^2)
    σrgyb = sqrt(σrg^2 + σyb^2)

    return σrgyb + 0.3 * μrgyb

end

hasler_and_susstrunk_m3(img::GenericGrayImage) = 0

"""
     M =  colorfulness(img)
     M =  colorfulness(HASLER_AND_SUSSTRUNK_M3(), img)

Measures the colorfulness of a natural image. Uses the
[`HASLER_AND_SUSSTRUNK_M3`](@ref) method by default.

See also: [`HASLER_AND_SUSSTRUNK_M3`](@ref).

"""
colorfulness(m::HASLER_AND_SUSSTRUNK_M3, img) = m(img)

colorfulness(img) = HASLER_AND_SUSSTRUNK_M3()(img)


text:
#=
formulates and solves the (dual of the) polynomial envelope problem described in:
D. Papp and S. Yildiz. Sum-of-squares optimization without semidefinite programming
available at https://arxiv.org/abs/1712.01792
=#

using SparseArrays

struct PolyEnvelopeNative{T <: Real} <: ExampleInstanceNative{T}
    n::Int
    rand_halfdeg::Int
    num_polys::Int
    env_halfdeg::Int
    primal_wsos::Bool # use primal formulation, else use dual
end

function build(inst::PolyEnvelopeNative{T}) where {T <: Real}
    (n, num_polys) = (inst.n, inst.num_polys)
    @assert inst.rand_halfdeg <= inst.env_halfdeg
    domain = PolyUtils.BoxDomain{T}(-ones(T, n), ones(T, n))

    # generate interpolation
    (U, pts, Ps, _, w) = PolyUtils.interpolate(domain,
        inst.env_halfdeg, get_quadr = true)

    # generate random data
    L = binomial(n + inst.rand_halfdeg, n)
    c_or_h = vec(Ps[1][:, 1:L] * rand(T(-9):T(9), L, num_polys))

    if inst.primal_wsos
        # WSOS cone in primal
        c = -w
        A = zeros(T, 0, U)
        b = T[]
        G = repeat(sparse(one(T) * I, U, U), outer = (num_polys, 1))
        h = c_or_h
    else
        # WSOS cone in dual
        c = c_or_h
        A = repeat(sparse(one(T) * I, U, U), outer = (1, num_polys))
        b = w
        G = Diagonal(-one(T) * I, num_polys * U) # TODO uniformscaling
        h = zeros(T, num_polys * U)
    end

    cones = Cones.Cone{T}[Cones.WSOSInterpNonnegative{T, T}(U, Ps,
        use_dual = !inst.primal_wsos) for k in 1:num_polys]

    model = Models.Model{T}(c, A, b, G, h, cones)
    return model
end


text:
# example using network from the RSSA Models
# download input files from: https://gist.github.com/isaacsas/7648fe065f1884307a2906a0a48ea12b

using DiffEqBase, Catalyst, OrdinaryDiffEq, Sundials
using ReactionNetworkImporters
using TimerOutputs

# parameters
networkname = "tester"
tf = 10.

# input files, specify path
datadir  = joinpath(@__DIR__,"../data/repressilator")
fname    = joinpath(datadir, "Repressilator.net")

const to = TimerOutput()
reset_timer!(to)

# get the reaction network
@timeit to "netgen" prn = loadrxnetwork(RSSAFile(), networkname, speciesf, rxsf; printrxs = false)
rn = prn.rn; initialpop = prn.u₀
@timeit to "addodes" addodes!(rn; build_jac=false, build_symfuncs=false, build_paramjac=false)
@timeit to "ODEProb" oprob = ODEProblem(rn,convert.(Float64,initialpop),(0.,tf))
show(to)
println()

# note solvers run faster the second time

# I haven't been able to successfully solve the system with Rodas4/Rodas5
#reset_timer!(to); @timeit to "Rodas4" begin sol = solve(oprob,Rodas4(autodiff=false),dense=false,calck=false); end; show(to)
#reset_timer!(to); @timeit to "Rodas4" begin sol = solve(oprob,Rodas4(autodiff=false),dense=false,calck=false); end; show(to)
#reset_timer!(to); @timeit to "Rodas4" begin sol = solve(oprob,Rodas4(autodiff=false),dense=false); end;

# similarly CVODE_BDF with gmres is very slow and I haven't been able to complete a simulation
#reset_timer!(to); @timeit to "CVODE_BDF" begin sol = solve(oprob,CVODE_BDF(linear_solver=:GMRES),dense=false); end; show(to)
#reset_timer!(to); @timeit to "CVODE_BDF" begin sol = solve(oprob,CVODE_BDF(linear_solver=:GMRES),dense=false); end; show(to)

# CVODE_BDF with LU seems to finish, on my machine it takes upwards of 500 seconds
reset_timer!(to); @timeit to "CVODE_BDF" begin sol = solve(oprob,CVODE_BDF(),dense=false); end; show(to)
reset_timer!(to); @timeit to "CVODE_BDF" begin sol = solve(oprob,CVODE_BDF(),dense=false); end; show(to)


text:
#####################################################################################################
#                                         Load Packages and Functions
#####################################################################################################
cd(@__DIR__)
using Pkg
Pkg.activate("..")
using GlobalMatchingModels, StatsPlots, Random
include("MINERVA_Functions.jl")
Random.seed!(356355)
#####################################################################################################
#                                        Specify  Model and Simulate Experiment
#####################################################################################################
# initialize model
model = MINERVA(n_features=20, n_traces=60, L=.50)
# simulate experiment
sim_data = mapreduce(x->simulate(model, 4, 5), hcat, 1:10000)'
#####################################################################################################
#                                      Plot Simulation
#####################################################################################################
pyplot()
density(sim_data, alpha=.9, grid=false, norm=true, legendtitle="Frequency", label=[0:5;]',
    legendtitlefontsize=9, xaxis=font(12), yaxis=font(12), xlabel="Echo Intensity", ylabel="Density",
    linewidth=1.5, bins=50, size=(600,300))
savefig("MINERVA.png")

# function noise_distribution(model, stimuli, probe, n_sim)
#     intensity = fill(0.0, n_sim)
#     for i in 1:n_sim
#         encode!(model, stimuli)
#         intensity[i] = compute_intensity(model, probe)
#     end
#     return intensity
# end

# function signal_distribution(model, stimuli, probe, n_sim)
#     intensity = fill(0.0, n_sim)
#     for i in 1:n_sim
#         encode!(model, stimuli)
#         intensity[i] = compute_intensity(model, probe)
#     end
#     return intensity
# end

# model = MINERVA(h=.75, n_features=30, n_traces=20, L=.50)

# stimuli = generate_stimuli(model)

# target_probe = random_target(stimuli)
# intensity = signal_distribution(model, stimuli, target_probe, 10_000)
# println(mean(intensity), " ", std(intensity))

# distractor_probe = random_distractor(model)
# intensity = signal_distribution(model, stimuli, target_probe, 10_000)
# println(mean(intensity), " ", std(intensity))

# function random_target(stimuli)
#     idx = rand(1:size(stimuli, 1))
#     return stimuli[idx,:]
# end

# function generate_stimuli(model)
#     @unpack n_traces, n_features = model
#     prototype = random_distractor(n_features)
#     stimuli = fill(0, n_traces, n_features)
#     for t in 1:n_traces, f in 1:n_features
#         stimuli[t,f] = stimulus_similarity(model, prototype[f])
#     end
#     return stimuli
# end


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

"""
    TTree <: AbstractVector{Any}

`UpROOT.TTree` is a wrapper around Python objects with mix-in
`uproot.tree.TTreeMethods`. It behaves like a Julia `AbstractVector` and
`Tables.Table` (with column access).

`TTree` is an on-disk data structure, iteration over single elements is *very*
inefficient.

Even though `TTree` is semantically a vector over the tree entries, a two
dimensional `getindex(tree::TTree, entries, branches)` is defined for
convenience. So you may use, e.g.:

```julia
tree[1:5, (:branch1, :branch2)]
tree[:, [:branch1, :branch2]]
tree[[1,4,7,10], ["branch1", "branch2"]]
```

Limitations: Write access is not implemented yet.
"""
struct TTree <: AbstractVector{NamedTuple}
    pyobj::PyObject

    function TTree(x::PyObject)
        if !pybuiltin(:isinstance)(x, uproot.tree.TTreeMethods)
            throw(ArgumentError("Python object $x is not a ROOT TTree"))
        end
        new(x)
    end
end

export TTree

pyobj(tree::TTree) = Base.getfield(tree, :pyobj)

Base.propertynames(tree::TTree) = Symbol.(pyobj(tree).keys())

Base.getproperty(tree::TTree, name::Symbol) = TBranch(get(pyobj(tree), String(name)))

Base.size(tree::TTree) = (pyobj(tree).numentries,)

Base.IndexStyle(tree::TTree) = IndexLinear()

function Base.getindex(tree::TTree, idxs::AbstractUnitRange)
    @boundscheck checkbounds(tree, idxs)
    cols = pyobj(tree).arrays(entrystart = first(idxs) - 1, entrystop = last(idxs))

    colnames = Tuple(Symbol.(keys(cols)))
    colvals_py = values(cols)
    colvals_jl = _branch_content.(py2jl.(colvals_py))

    Table(NamedTuple{colnames}(colvals_jl))
end

Base.getindex(tree::TTree, i::Integer) = first(tree[i:i])

Base.getindex(tree::TTree, ::Colon) = getindex(tree, eachindex(tree))

Base.getindex(tree::TTree, idxs, col::Symbol) = getproperty(tree, col)[idxs]

function Base.getindex(tree::TTree, idxs, cols::Union{Tuple,AbstractArray})
    colnames = map(Symbol, (cols...,))
    colvalues = map(c -> tree[idxs, c], colnames)
    Table(NamedTuple{colnames}(colvalues))
end

Tables.istable(::Type{TTree}) = true

Tables.columnaccess(::Type{TTree}) = true

Tables.columns(tree::TTree) = Dict(name => getproperty(tree, name) for name in Tables.columnnames(tree))

function Tables.schema(tree::TTree)
    props = Base.propertynames(tree)
    Tables.Schema(props, map(_ -> Any, props))
end

TypedTables.Table(tree::TTree) = tree[:]

"""
    TBranch <: AbstractVector{Any}

`UpROOT.TBranch` is a wrapper around Python objects with mix-in
`uproot.tree.TBranchMethods`. It behaves like a Julia `AbstractVector`. Ff
the branch has children, it also behaves like a `Tables.Table` (with column
access).

`TBranch` is an on-disk data structure, iteration over single elements is *very*
inefficient.

Limitations: Write access is not implemented yet.
"""
struct TBranch <: AbstractVector{Any}
    pyobj::PyObject

    function TBranch(x::PyObject)
        if !pybuiltin(:isinstance)(x, uproot.tree.TBranchMethods)
            throw(ArgumentError("Python object $x is not a ROOT TBranch"))
        end
        new(x)
    end
end

export TBranch

Base.propertynames(branch::TBranch) = Symbol.(pyobj(branch).keys())

Base.getproperty(branch::TBranch, name::Symbol) = TBranch(get(pyobj(branch), String(name)))

Base.size(branch::TBranch) = (pyobj(branch).numentries,)

Base.IndexStyle(branch::TBranch) = IndexLinear()

function Base.getindex(branch::TBranch, idxs::AbstractUnitRange)
    @boundscheck checkbounds(branch, idxs)
    content = pyobj(branch).array(entrystart = first(idxs) - 1, entrystop = last(idxs))
    _branch_content(py2jl(content))
end

Base.getindex(branch::TBranch, i::Integer) = first(branch[i:i])

Base.getindex(branch::TBranch, ::Colon) = getindex(branch, eachindex(branch))

Base.copy(branch::TBranch) = Array(branch)

Base.Array(branch::TBranch) = branch[:]

Tables.istable(::Type{TBranch}) = true

Tables.columnaccess(::Type{TBranch}) = true

Tables.columns(branch::TBranch) = branch # Dict(broadcast(item -> Symbol(item[1]) => TBranch(item[2]), pyobj(branch).items()))

function Tables.schema(branch::TBranch)
    props = Base.propertynames(branch)
    Tables.Schema(props, map(_ -> Any, props))
end

function read_ttree(tree::PyObject, branchnames::Array{<:AbstractString})
    tree_data = tree.arrays(branchnames)

    df = DataFrame()
    for bn in branchnames
        df[Symbol(bn)] = py2jl(tree_data[bn])
    end

    df
end

function read_ttree(filename::AbstractString, treename::AbstractString)
    file = uproot.open(filename)
    tree = get(file, treename)
    branchnames = tree.keys()
    read_ttree(tree, branchnames)
end

function read_ttree(filename::AbstractString, treename::AbstractString, branchnames::Array{<:AbstractString})
    file = uproot.open(filename)
    tree = get(file, treename)
    read_ttree(tree, branchnames)
end

function read_ttree(filename::AbstractString, treename::AbstractString, branchnames::Array{Symbol})
    read_ttree(filename, treename, String.(branchnames))
end

_branch_content(A::AbstractVector) = A

# uproot returns some branches containing fix-size n-dim array data as
# (n+1)-dim arrays. Permutation of dimensions in _branch_content changes
# memory layout and preserves order of indices. So a branch containing
# static C arrays of size [2][3] becomes a vector of arrays of size (2, 3) in
# Julia.
function _branch_content(A::AbstractArray)
    dimidxs = ntuple(i -> i, Val(length(size(A))))
    branchdim = first(dimidxs)
    innerdims = Base.tail(dimidxs)
    VectorOfSimilarArrays(permutedims(A, (innerdims..., branchdim)))
end


text:
using MatterPower
using PyCall
using Roots
using Tables, CSV
# %% Call the python wrapper for CLASS, `classy`, via PyCall
classy = pyimport("classy")
# Create an instance of the CLASS wrapper
cosmo = classy.Class()
# Create a dictionary of the cosmological parameters
deg_ncdm = true # 3 massive neutrinos with degenerate mass?
params = Dict(
   "output" => "mPk",
   "P_k_max_h/Mpc" => 100,
   "z_pk" => 5,
   "non linear" => "halofit",
   "A_s" => 2.097e-9,
   "n_s" => 0.9652,
   "k_pivot" => 0.05,
   "h" => 0.6737,
   "omega_b" => 0.02233,
   "omega_cdm" => 0.1198,
   "N_ncdm" => 1,
)
if deg_ncdm
   push!(params, "m_ncdm" => 0.02, "deg_ncdm" => 3, "N_ur" => 0.00641)
   mν = params["m_ncdm"] * params["deg_ncdm"]
else
   push!(params, "m_ncdm" => 0.06, "N_ur" => 2.0328)
   mν = params["m_ncdm"]
end
h0 = params["h"]
Ωcb = (params["omega_b"] + params["omega_cdm"]) / h0^2
# Set the parameters and run the CLASS code
cosmo.set(params)
cosmo.compute()

# %% Compute non-linear mass, Mstar, at various redshifts
redshift = 0:0.1:2
nred = length(redshift)
Rstar = zeros(nred)
Mstar = zeros(nred)
for ired = 1:nred
   z = redshift[ired]
   # Define functions to return a linear baryon+CDM power spectrum
   # Note: The CLASS code takes wavenumbers in units of 1/Mpc (no h) and
   # return power spectra in units of Mpc^3 (no 1/h^3).
   pkcb(kovh) = cosmo.pk_cb_lin(kovh * h0, z) * h0^3
   # %% Compute the non-linear mass, M*, defined by σ(R*,z) = 1.6865
   # and M* = (4π/3)ρm R*^3 = (4π/3)ρc Ωm R*^3
   # ρc is the present-day critical density of the Universe in units of h^2 M⊙/Mpc^3
   f(R) = √sigma2(pkcb, R) - 1.6865
   Rstar[ired] = find_zero(f, (0.2, 3), Bisection())
end
ρc = 2.775e11
Mstar = 4π / 3 * ρc * Ωcb * Rstar .^ 3

#%% Save table data to table1.csv
t = Tables.table([redshift Rstar Mstar])
header = ["z", "Rstar", "Mstar"]
CSV.write("redshift_Rstar_Mstar.csv", t, header = header)

# %% Clean CLASS (the equivalent of the struct_free() in the `main`
# of CLASS. This step is primordial when running in a loop over different
# cosmologies, as you will saturate your memory very fast if you ommit
# it.
cosmo.struct_cleanup()
# If you want to change completely the cosmology, you should also
# clean the arguments, otherwise, if you are simply running on a loop
# of different values for the same parameters, this step is not neededc
cosmo.empty()


text:
using Images

BufferedImage = Union{Void,AbstractArray}
Stream = Union{Void,Base.LibuvStream}

getWidth(image::AbstractArray) = size(image)[1]
getHeight(image::AbstractArray) = size(image)[2]

printStackTrace(ex) = Base.showerror(STDOUT, ex, catch_backtrace())

type Move
  fromX::Int
  fromY::Int
  toX::Int
  toY::Int
  Move(fromX::Int, fromY::Int, toX::Int, toY::Int) = new(fromX, fromY, toX, toY)
end

toString(this::Move) = string(this.fromX, ",", this.fromY, " -> ", this.toX, ",", this.toY)

type PrintWriter
  stream::Stream
end

pw_println(this::PrintWriter, args...) = println(this.stream, args...)
pw_flush(this::PrintWriter) = flush(this.stream)

type Socket
  hostname::String
  port::Int
  stream::Stream
  Socket() = new("",0,nothing)
  function Socket(hostname::String, port::Int)
    this = new(hostname,port,nothing)
    this.stream  = connect(hostname, port)
    this
  end
end

type SocketStream
  socket::Union{Void, Socket}
  SocketStream() = new(nothing)
  SocketStream(x::Socket) = new(x)
end

InputStream = SocketStream
OutputStream = SocketStream

STREAM = []

function st_read(this::InputStream)
  value = 0
  #@async if isopen(this.socket.stream) value = parse(Integer, string(read(this.socket.stream)...)) end # parse(Integer,)  while isopen(this.socket.stream) println(readline(this.socket.stream)) end
  #@async while true push!(value, read(this.socket.stream)) end
  #@async while isopen(this.socket.stream) println(readline(this.socket.stream)) end
  value = read(this.socket.stream, 1)[1]
  println("Value: ", value)
  value
end

function st_write(this::OutputStream, value::Any)
  if isopen(this.socket.stream) write(this.socket.stream,value) end
end

st_flush(this::OutputStream) = flush(this.socket.stream)

getOutputStream(this::Socket) = SocketStream(this)
getInputStream = getOutputStream

SystemCurrentTimeMillis() = Dates.time() * 1000

type NetworkClient
  jdField_a_of_type_JavaNetSocket::Socket
  jdField_a_of_type_JavaIoInputStream::InputStream
  jdField_a_of_type_JavaIoOutputStream::OutputStream
  jdField_a_of_type_Int::Int
  b::Int
  c::Int
  NetworkClient() = new(Socket(),InputStream(),OutputStream(),0,0)
end

function NetworkClient(hostname::String, teamName::String, logo::BufferedImage)
  this = NetworkClient()
  if (logo == nothing) || (getWidth(logo) != 256) || (getHeight(logo) != 256)
    throw(ErrorException("You have to provide a 256x256 image as your team logo")) #RuntimeException
  end
  this.jdField_a_of_type_JavaNetSocket = Socket(hostname, 22135)
  this.jdField_a_of_type_JavaIoOutputStream = getOutputStream(this.jdField_a_of_type_JavaNetSocket)
  this.jdField_a_of_type_JavaIoInputStream = getInputStream(this.jdField_a_of_type_JavaNetSocket)
  st_write(this.jdField_a_of_type_JavaIoOutputStream, 1)
  st_flush(this.jdField_a_of_type_JavaIoOutputStream)
  if Int8(st_read(this.jdField_a_of_type_JavaIoInputStream)) != 1
    throw(ErrorException("Outdated client software - update your client!"))
  end
  pw = PrintWriter(this.jdField_a_of_type_JavaIoOutputStream.socket.stream)
  pw_println(pw, teamName)
  pw_println(pw, encodeImageBase64(logo))
  pw_flush(pw)
  hostname = st_read(this.jdField_a_of_type_JavaIoInputStream)
  this.jdField_a_of_type_Int = (hostname & 0x3)
  this.b = trunc(hostname / 4)
  l = SystemCurrentTimeMillis()
  st_write(this.jdField_a_of_type_JavaIoOutputStream, 0)
  st_flush(this.jdField_a_of_type_JavaIoOutputStream)
  st_read(this.jdField_a_of_type_JavaIoInputStream)
  this.c = trunc(Int(SystemCurrentTimeMillis() - l) / 2)
  println("Expected network latency in milliseconds: ", this.c)
  this
end

function receiveMove(this::NetworkClient)::Move
  try
    i = 0
    if (i = st_read(this.jdField_a_of_type_JavaIoInputStream)) == 88 #no input
      return nothing
    end
    if i == 99 #invalid
      throw(ErrorException("You got kicked because your move was invalid!"))
    end
    j = st_read(this.jdField_a_of_type_JavaIoInputStream)
    return Move(i / 9, i % 9, j / 9, j % 9)
  catch localIOException #IOException
    throw(ErrorException(string("Failed to receive move: ", localIOException)))
  end
end

function sendMove(this::NetworkClient, move::Move)
  try
    st_write(this.jdField_a_of_type_JavaIoOutputStream, move.fromX * 9 + move.fromY)
    st_write(this.jdField_a_of_type_JavaIoOutputStream, move.toX * 9 + move.toY)
    st_flush(this.jdField_a_of_type_JavaIoOutputStream)
  catch localIOException #IOException
    throw(ErrorException(string("Failed to send move: ", localIOException)))
  end
end

getMyPlayerNumber(this::NetworkClient) = this.jdField_a_of_type_Int
getTimeLimitInSeconds(this::NetworkClient) = this.b
getExpectedNetworkLatencyInMilliseconds(this::NetworkClient) = this.c

function encodeImageBase64(image::BufferedImage)
  io = IOBuffer()
  save(Images.Stream(Images.format"PNG", io), image) # color image is still a gray image on server
  base64encode(io.data)
end

#img = Gray.(img)
#imgMatrix = vec(UInt8.(reinterpret.(channelview(img))))
#convert(Array{UInt8},Images.raw(grayImage))
#grayimg = (c -> begin g=N0f8((c.r+c.g+c.b)/3); r=RGBA(g,g,g,N0f8(1.0)); end).(img)
#save("_temp.png", img)
#content=""
#open("_temp.png") do file
#   content=readstring(file)
#end
#data = convert(Vector{UInt8},content))