text:
using CayleyDickson
using Base.Test: @test, @test_throws

@test begin
    a = CayleyDickson.randomBigFloat()
    isreal(Exo3Real(a))
end

@test begin
    a = 1
    !isreal(Exo3Real(a, a, a, a))
end

@test begin
    z = random(Exo3Real{Int})
    z == +(z)
end

@test begin
    a = 1
    real(Exo3Real(a)) == a
end

@test begin
    a = 1
    l = unreal(Exo3Real(0, a, a, a, a, a, a, a))
    r = [a, a, a, a, a, a, a]
    l == r
end

@test begin
    io = IOBuffer()
    show(io, Exo3Real(1.0, 2, 3, 4, 5, 6, 7, 8))
    l = String(take!(io))
    r = "[1: 1.0, a: 2.0, b: 3.0, ab: 4.0, c: 5.0, ac: 6.0, bc: 7.0, (ab)c: 8.0]"
    l == r
end

@test begin
    io = IOBuffer()
    show(io, Exo3Real(1.0, 2, 3, 4, 5, 6, 7))
    l = String(take!(io))
    r = "[1: 1.0, a: 2.0, b: 3.0, ab: 4.0, c: 5.0, ac: 6.0, bc: 7.0, (ab)c: 0.0]"
    l == r
end

@test begin
    io = IOBuffer()
    show(io, Exo3Real(1.0, 2, 3, 4, 5, 6))
    l = String(take!(io))
    r = "[1: 1.0, a: 2.0, b: 3.0, ab: 4.0, c: 5.0, ac: 6.0, bc: 0.0, (ab)c: 0.0]"
    l == r
end

@test begin
    io = IOBuffer()
    show(io, Exo3Real(1.0, 2, 3, 4, 5))
    l = String(take!(io))
    r = "[1: 1.0, a: 2.0, b: 3.0, ab: 4.0, c: 5.0, ac: 0.0, bc: 0.0, (ab)c: 0.0]"
    l == r
end

@test begin
    io = IOBuffer()
    show(io, Exo3Real(1.0, 2, 3, 4))
    l = String(take!(io))
    r = "[1: 1.0, a: 2.0, b: 3.0, ab: 4.0, c: 0.0, ac: 0.0, bc: 0.0, (ab)c: 0.0]"
    l == r
end

@test begin
    a = CayleyDickson.randomBigInt()
    b = CayleyDickson.randomBigFloat()
    c = CayleyDickson.randomBigInt()
    d = CayleyDickson.randomBigFloat()
    l = Exo3Real(Exo2Real(a, b, c, d))
    r = Exo3Real(a, b, c, d)
    l == r
end

@test begin
    a = CayleyDickson.randomBigInt()
    b = CayleyDickson.randomBigFloat()
    l = Exo3Real(Exo1Real(a, b))
    r = Exo3Real(a, b)
    l == r
end

@test begin
    a = CayleyDickson.randomBigInt()
    b = CayleyDickson.randomBigFloat()
    c = CayleyDickson.randomBigInt()
    l = Exo3Real(a, b, c)
    r = Exo3Real(a, b, c, 0)
    l == r
end

@test begin
    x = random(Exo3Real{BigFloat})
    y = random(Exo3Real{BigFloat})
    x != y
end

@test begin
    a = CayleyDickson.randomBigInt()
    z = random(Exo3Real{BigInt})
    l = a + z
    r = z + a
    l == r
end

@test begin
    a = random(Exo2Real{BigInt})
    z = random(Exo3Real{BigInt})
    l = a + z
    r = z + a
    l == r
end

@test begin
    a = CayleyDickson.randomBigInt()
    z = random(Exo3Real{BigInt})
    l = a - z
    r = -(z - a)
    l == r
end

@test begin
    a = random(Exo2Real{BigInt})
    z = random(Exo3Real{BigInt})
    l = a - z
    r = -(z - a)
    l == r
end

@test begin
    x = random(Exo3Real{BigInt})
    y = random(Exo3Real{BigInt})
    l = x + y
    r = y + x
    l == r
end

@test begin
    x = random(Exo3Real{BigInt})
    y = random(Exo3Real{BigInt})
    l = x * y
    r = y * x
    l != r
end

@test begin
    x = random(Exo3Real{BigInt})
    y = random(Exo3Real{BigInt})
    l = x - y
    r = -(y - x)
    l == r
end

@test begin
    x = random(Exo3Real{BigInt})
    y = random(Exo3Real{BigInt})
    z = random(Exo3Real{BigInt})
    l = (x + y) + z
    r = x + (y + z)
    l == r
end

@test begin
    x = random(Exo3Real{BigInt})
    y = random(Exo3Real{BigInt})
    z = random(Exo3Real{BigInt})
    l = (x * y) * z
    r = x * (y * z)
    l != r
end

@test begin
    x = random(Exo3Real{BigInt})
    y = random(Exo3Real{BigInt})
    l = alternatorL(x, y)
    r = zero(Exo3Real{BigInt})
    l == r
end

@test begin
    x = random(Exo3Real{BigInt})
    y = random(Exo3Real{BigInt})
    l = alternatorR(x, y)
    r = zero(Exo3Real{BigInt})
    l == r
end

@test begin
    x = random(Exo3Real{BigInt})
    conj(conj(x)) == x
end

@test begin
    x = random(Exo3Real{Rational{BigInt}})
    inv(inv(x)) == x
end

@test begin
    x = random(Exo3Real{BigInt})
    y = random(Exo3Real{BigInt})
    l = conj(x * y)
    r = conj(y) * conj(x)
    l == r
end

@test begin
    x = random(Exo3Real{BigInt})
    y = random(Exo3Real{BigInt})
    l = conj(x + y)
    r = conj(x) + conj(y)
    l == r
end

@test begin
    x = random(Exo3Real{BigInt})
    y = random(Exo3Real{BigInt})
    l = conj(x - y)
    r = conj(x) - conj(y)
    l == r
end

@test begin
    x = random(Exo3Real{Rational{BigInt}})
    y = random(Exo3Real{Rational{BigInt}})
    l = inv(x * y)
    r = inv(y) * inv(x)
    l == r
end

@test begin
    x = random(Exo3Real{Rational{BigInt}})
    y = random(Exo3Real{Rational{BigInt}})
    l = inv(x / y)
    r = y / x
    l == r
end

@test begin
    x = random(Exo3Real{Rational{BigInt}})
    y = random(Exo3Real{Rational{BigInt}})
    l = inv(x \ y)
    r = y \ x
    l == r
end

@test begin
    x = random(Exo3Real{BigInt})
    y = random(Exo3Real{BigInt})
    z = random(Exo3Real{BigInt})
    l = x * (y + z)
    r = (x * y) + (x * z)
    l == r
end

@test begin
    a = CayleyDickson.randomBigInt()
    x = random(Exo3Real{BigInt})
    y = random(Exo3Real{BigInt})
    l = a * (x + y)
    r = (a * x) + (a * y)
    l == r
end

@test begin
    a = CayleyDickson.randomBigInt()
    x = random(Exo3Real{BigInt})
    l = a * x
    r = x * a
    l == r
end

@test begin
    a = random(Exo2Real{BigInt})
    x = random(Exo3Real{BigInt})
    l = commutator(a, x)
    r = zero(Exo3Real{BigInt})
    l != r
end

@test begin
    x = random(Exo3Real{BigInt})
    y = random(Exo3Real{BigInt})
    z = random(Exo3Real{BigInt})
    l = (x + y) * z
    r = (x * z) + (y * z)
    l == r
end

@test begin
    x = random(Exo3Real{BigInt})
    y = random(Exo3Real{BigInt})
    z = random(Exo3Real{BigInt})
    l = x * (y - z)
    r = (x * y) - (x * z)
    l == r
end

@test begin
    x = random(Exo3Real{BigInt})
    y = random(Exo3Real{BigInt})
    z = random(Exo3Real{BigInt})
    l = (x - y) * z
    r = (x * z) - (y * z)
    l == r
end

@test begin
    x = random(Exo3Real{BigInt})
    quadrance(x) > big(0) || quadrance(x) < big(0) || quadrance(x) == big(0)
end

@test begin
    x = random(Exo3Real{BigInt})
    y = random(Exo3Real{BigInt})
    l = quadrance(x * y)
    r = quadrance(x) * quadrance(y)
    l == r
end

@test_throws ErrorException begin
    inv(Exo3Real(0, 0.0))
end

@test_throws ErrorException begin
    random(Exo3Real{Int}) / 0
end

@test_throws ErrorException begin
    0 \ random(Exo3Real{Int})
end

@test begin
    a = CayleyDickson.randomBigRational()
    z = random(Exo3Real{Rational{BigInt}})
    l = z / a
    r = a \ z
    l == r
end

@test begin
    a = CayleyDickson.randomBigRational()
    z = random(Exo3Real{Rational{BigInt}})
    l = a / z
    r = z \ a
    l == r
end

@test begin
    z = random(Exo3Real{BigInt})
    z == z + zero(Exo3Real{BigInt})
end

@test begin
    z = random(Exo3Real{BigInt})
    z == z * one(Exo3Real{BigInt})
end

@test begin
    z = random(Exo3Real{BigInt})
    z == z + zero(z)
end

@test begin
    z = random(Exo3Real{BigInt})
    z == z * one(z)
end

@test begin
    x = random(Exo3Real{BigInt})
    cloak(cloak(x)) == x
end

@test begin
    x = random(Exo3Real{BigInt})
    dagger(dagger(x)) == x
end

@test begin
    x = random(Exo3Real{BigInt})
    star(star(x)) == x
end

@test begin
    x = random(Exo3Real{Rational{BigInt}})
    selfstar(selfstar(x)) == selfstar(x)
end

@test begin
    x = random(Exo3Real{Rational{BigInt}})
    anti_selfstar(anti_selfstar(x)) == anti_selfstar(x)
end

@test begin
    A = Exo3Real(0, 1)
    B = Exo3Real(0, 0, 1)
    C = Exo3Real(0, 0, 0, 0, 1)
    l = jacobiator(A, B, C)
    r = zero(Exo3Real{Int})
    l != r
end

@test begin
    A = Exo3Real(0, 1)
    B = Exo3Real(0, 0, 1)
    C = Exo3Real(0, 0, 0, 0, 1)
    l = star(star(A * B) * C) + star(star(B * C) * A) + star(star(C * A) * B)
    r = zero(Exo3Real{Int})
    l == r
end


text:
@testset "Pagerank" begin
    function dense_pagerank_solver(g::AbstractGraph, α=0.85::Real)
        p = fill(1 / nv(g), nv(g))
        danglingnodes = outdegree(g) .== 0
        M = Matrix{Float64}(adjacency_matrix(g))
        M = M'
        M[:, danglingnodes] .= sum(danglingnodes) ./ nv(g)
        M = M * Diagonal(1 ./ sum(M, dims=1)[:])
        @assert all(1.01 .>= sum(M, dims=1) .>= 0.999)
       # v = inv(I-β*M) * ((1-β)/nv(g) * ones(nv(g), 1))
        v = inv(I - α * M) * ((1 - α) / nv(g) * ones(nv(g), 1))
        return v
    end

    g5 = SimpleDiGraph(4)
    add_edge!(g5, 1, 2); add_edge!(g5, 2, 3); add_edge!(g5, 1, 3); add_edge!(g5, 3, 4)
    g6 = SimpleGraph(4)
    add_edge!(g6, 1, 2); add_edge!(g6, 2, 3); add_edge!(g6, 1, 3); add_edge!(g6, 3, 4)
    for α in [0.75, 0.85]

        for g in testdigraphs(g5)
            @test pagerank(g)[3] ≈ 0.318 atol = 0.001
            @test length(@inferred(pagerank(g))) == nv(g)
            @test_throws ErrorException pagerank(g, 2)
            @test_throws ErrorException pagerank(g, α, 2)
            @test isapprox(pagerank(g, α), dense_pagerank_solver(g, α), atol=0.001)
        end

        for g in testgraphs(g6)
            @test length(@inferred(pagerank(g))) == nv(g)
            @test_throws ErrorException pagerank(g, 2)
            @test_throws ErrorException pagerank(g, α, 2)
            @test isapprox(pagerank(g, α), dense_pagerank_solver(g, α), atol=0.001)
        end
    end
end

ext: jl


text:
@testset "neighbor features" begin
    scene=Frame(Entity{VehicleState, BicycleModel, Int},100)
    roadway=gen_straight_roadway(3, 200.0, lane_width=3.0)
    push!(scene,Entity(VehicleState(VecSE2(30.0,3.0,0.0), roadway, 0.0),
        BicycleModel(VehicleDef(AgentClass.CAR, 4.826, 1.81)),1))
    push!(scene,Entity(VehicleState(VecSE2(20.0,3.0,0.0), roadway, 0.0),
        BicycleModel(VehicleDef(AgentClass.CAR, 4.826, 1.81)),1))
    push!(scene,Entity(VehicleState(VecSE2(40.0,3.0,0.0), roadway, 0.0),
        BicycleModel(VehicleDef(AgentClass.CAR, 4.826, 1.81)),1))
    push!(scene,Entity(VehicleState(VecSE2(20.0,0.0,0.0), roadway, 0.0),
        BicycleModel(VehicleDef(AgentClass.CAR, 4.826, 1.81)),1))
    push!(scene,Entity(VehicleState(VecSE2(40.0,0.0,0.0), roadway, 0.0),
        BicycleModel(VehicleDef(AgentClass.CAR, 4.826, 1.81)),1))
    push!(scene,Entity(VehicleState(VecSE2(20.0,6.0,0.0), roadway, 0.0),
        BicycleModel(VehicleDef(AgentClass.CAR, 4.826, 1.81)),1))
    push!(scene,Entity(VehicleState(VecSE2(40.0,6.0,0.0), roadway, 0.0),
        BicycleModel(VehicleDef(AgentClass.CAR, 4.826, 1.81)),1))
    @test get_neighbor_fore_along_lane(scene,1,roadway) == NeighborLongitudinalResult(3,10.0)
    @test get_neighbor_rear_along_lane(scene,1,roadway) == NeighborLongitudinalResult(2,10.0)
    @test get_neighbor_fore_along_left_lane(scene,1,roadway) == NeighborLongitudinalResult(7,10.0)
    @test get_neighbor_rear_along_left_lane(scene,1,roadway) == NeighborLongitudinalResult(6,10.0)
    @test get_neighbor_fore_along_right_lane(scene,1,roadway) == NeighborLongitudinalResult(5,10.0)
    @test get_neighbor_rear_along_right_lane(scene,1,roadway) == NeighborLongitudinalResult(4,10.0)

    trajdata = get_test_trajdata(roadway)
    @test get_neighbor_fore_along_lane(get!(Scene(), trajdata, 1), 1, roadway) == NeighborLongitudinalResult(2, 3.0)
    @test get_neighbor_fore_along_lane(get!(Scene(), trajdata, 1), 2, roadway) == NeighborLongitudinalResult(nothing, 250.0)
    @test get_neighbor_fore_along_lane(get!(Scene(), trajdata, 2), 1, roadway) == NeighborLongitudinalResult(2, 4.0)
    @test get_neighbor_fore_along_lane(get!(Scene(), trajdata, 2), 2, roadway) == NeighborLongitudinalResult(nothing, 250.0)

    roadway = gen_stadium_roadway(1)
    scene = Scene(2)
    scene.n = 2
    def = VehicleDef(AgentClass.CAR, 2.0, 1.0)

    function place_at!(i, s)
        roadproj = proj(VecSE2(0.0,0.0,0.0), roadway)
        roadind = RoadIndex(roadproj.curveproj.ind, roadproj.tag)
        roadind = move_along(roadind, roadway, s)
        frenet = Frenet(roadind, roadway[roadind].s, 0.0, 0.0)
        state = VehicleState(frenet, roadway, 0.0)
        scene[i] = Vehicle(state, def, i)
    end

    place_at!(1, 0.0)
    place_at!(2, 0.0)
    foreinfo = get_neighbor_fore_along_lane(scene, 1, roadway, max_distance_fore=Inf)
    @test isapprox(foreinfo.Δs, 0.0)

    place_at!(2, 100.0)
    foreinfo = get_neighbor_fore_along_lane(scene, 1, roadway, max_distance_fore=Inf)
    @test foreinfo.ind == 2
    @test isapprox(foreinfo.Δs, 100.0)

    foreinfo = get_neighbor_fore_along_lane(scene, 2, roadway, max_distance_fore=Inf)
    @test foreinfo.ind == 1
    @test isapprox(foreinfo.Δs, 277.07, atol=1e-2)

    place_at!(2, 145.0)
    foreinfo = get_neighbor_fore_along_lane(scene, 1, roadway, max_distance_fore=Inf)
    @test foreinfo.ind == 2
    @test isapprox(foreinfo.Δs, 145.0, atol=1e-5)

    place_at!(2, 240.0)
    foreinfo = get_neighbor_fore_along_lane(scene, 1, roadway, max_distance_fore=Inf)
    @test foreinfo.ind == 2
    @test isapprox(foreinfo.Δs, 240.0, atol=1e-5)

    place_at!(1, 240.0)
    foreinfo = get_neighbor_fore_along_lane(scene, 1, roadway, max_distance_fore=Inf)
    @test foreinfo.ind == 2
    @test isapprox(foreinfo.Δs, 0.0, atol=1e-5)

    ################################################
    # get_frenet_relative_position

    place_at!(1, 0.0)
    place_at!(2, 10.0)
    frp = get_frenet_relative_position(scene[2].state.posG, scene[1].state.posF.roadind, roadway, max_distance_fore=Inf)
    @test frp.origin == scene[1].state.posF.roadind
    @test frp.target.ind.i == 1
    @test isapprox(frp.target.ind.t, 0.1, atol=0.01)
    @test frp.target.tag == scene[1].state.posF.roadind.tag
    @test isapprox(frp.Δs, 10.0, atol=1e-2)
    @test isapprox(frp.t, 0.0, atol=1e-5)
    @test isapprox(frp.ϕ, 0.0, atol=1e-7)

    place_at!(1, 10.0)
    place_at!(2, 20.0)
    frp = get_frenet_relative_position(scene[2].state.posG, scene[1].state.posF.roadind, roadway, max_distance_fore=Inf)
    @test frp.origin == scene[1].state.posF.roadind
    @test frp.target.ind.i == 1
    @test isapprox(frp.target.ind.t, 0.2, atol=0.01)
    @test frp.target.tag == scene[1].state.posF.roadind.tag
    @test isapprox(frp.Δs, 10.0, atol=1e-2)
    @test isapprox(frp.t, 0.0, atol=1e-5)
    @test isapprox(frp.ϕ, 0.0, atol=1e-7)

    place_at!(1, 0.0)
    place_at!(2, 120.0)
    frp = get_frenet_relative_position(scene[2].state.posG, scene[1].state.posF.roadind, roadway, max_distance_fore=Inf)
    @test frp.target.tag == LaneTag(1, 1)
    @test scene[1].state.posF.roadind.tag == LaneTag(6, 1)
    @test isapprox(frp.Δs, 120.0, atol=1e-2)
    @test isapprox(frp.t, 0.0, atol=1e-5)
    @test isapprox(frp.ϕ, 0.0, atol=1e-7)

    place_at!(1, 0.0)
    place_at!(2, 250.0)
    frp = get_frenet_relative_position(scene[2].state.posG, scene[1].state.posF.roadind, roadway, max_distance_fore=Inf)
    @test frp.target.tag != scene[1].state.posF.roadind.tag
    @test isapprox(frp.Δs, 250.0, atol=1e-2)
    @test isapprox(frp.t, 0.0, atol=1e-5)
    @test isapprox(frp.ϕ, 0.0, atol=1e-7)
end

@testset begin "features interface"
    roadway = gen_straight_roadway(3, 1000.0, lane_width=1.0)
    rec = SceneRecord(1, 0.1, 5)
    update!(rec, Scene([
            Vehicle(VehicleState(VecSE2( 0.0,0.0,0.0), roadway, 10.0), VehicleDef(AgentClass.CAR, 5.0, 2.0), 1),
            Vehicle(VehicleState(VecSE2(10.0,0.0,0.0), roadway, 10.0), VehicleDef(AgentClass.CAR, 5.0, 2.0), 2),
        ]))

    @test convert(Float64, get(IS_COLLIDING, rec, roadway, 1)) == 0.0
    @test convert(Float64, get(IS_COLLIDING, rec, roadway, 2)) == 0.0

    @test isapprox(convert(Float64, get(MARKERDIST_LEFT, rec, roadway, 1)), 0.5)
    @test isapprox(convert(Float64, get(MARKERDIST_LEFT, rec, roadway, 2)), 0.5)
    @test isapprox(convert(Float64, get(MARKERDIST_RIGHT, rec, roadway, 1)), 0.5)
    @test isapprox(convert(Float64, get(MARKERDIST_RIGHT, rec, roadway, 2)), 0.5)
    @test isapprox(convert(Float64, get(MARKERDIST_LEFT_LEFT, rec, roadway, 1)), 1.5)
    @test isapprox(convert(Float64, get(MARKERDIST_LEFT_LEFT, rec, roadway, 2)), 1.5)
    @test isnan(convert(Float64, get(MARKERDIST_RIGHT_RIGHT, rec, roadway, 1)))
    @test isnan(convert(Float64, get(MARKERDIST_RIGHT_RIGHT, rec, roadway, 2)))

    @test isapprox(convert(Float64, get(DIST_FRONT, rec, roadway, 1)), 10.0 - 5.0)
    @test convert(Float64, get(DIST_FRONT, rec, roadway, 2, censor_hi=100.0)) == 100.0

    update!(rec, Scene([
            Vehicle(VehicleState(VecSE2( 1.0,0.0,0.0), roadway, 10.0), VehicleDef(AgentClass.CAR, 5.0, 2.0), 1),
            Vehicle(VehicleState(VecSE2(10.0,0.0,0.0), roadway, 10.0), VehicleDef(AgentClass.CAR, 5.0, 2.0), 2),
            Vehicle(VehicleState(VecSE2(12.0,1.0,0.0), roadway, 10.0), VehicleDef(AgentClass.CAR, 5.0, 2.0), 3),
            Vehicle(VehicleState(VecSE2( 0.0,1.0,0.0), roadway, 10.0), VehicleDef(AgentClass.CAR, 5.0, 2.0), 4),
        ]))

    @test isapprox(convert(Float64, get(MARKERDIST_LEFT, rec, roadway, 3)), 0.5)
    @test isapprox(convert(Float64, get(MARKERDIST_RIGHT, rec, roadway, 3)), 0.5)
    @test isapprox(convert(Float64, get(MARKERDIST_LEFT_LEFT, rec, roadway, 3)), 1.5)
    @test isapprox(convert(Float64, get(MARKERDIST_RIGHT_RIGHT, rec, roadway, 3)), 1.5)

    @test isapprox(convert(Float64, get(DIST_FRONT, rec, roadway, 1)), 9.0 - 5.0)
    @test isapprox(convert(Float64, get(DIST_FRONT_LEFT, rec, roadway, 1)), 11.0 - 5.0)
    @test convert(Float64, get(DIST_FRONT_RIGHT, rec, roadway, 1)) == 100.0
    @test isapprox(convert(Float64, get(DIST_FRONT, rec, roadway, 4)), 12.0 - 5.0)
    @test convert(Float64, get(DIST_FRONT_LEFT, rec, roadway, 4)) == 100.0
    @test isapprox(convert(Float64, get(DIST_FRONT_RIGHT, rec, roadway, 4)), 1.0 - 5.0)

    @test convert(Float64, get(IS_COLLIDING, rec, roadway, 1)) == 1.0
    @test convert(Float64, get(IS_COLLIDING, rec, roadway, 2)) == 1.0
    @test convert(Float64, get(IS_COLLIDING, rec, roadway, 3)) == 1.0
    @test convert(Float64, get(IS_COLLIDING, rec, roadway, 4)) == 1.0

    update!(rec, Scene([
            Vehicle(VehicleState(VecSE2(  0.0,0.0,0.0), roadway, 10.0), VehicleDef(AgentClass.CAR, 5.0, 2.0), 1),
            Vehicle(VehicleState(VecSE2(-10.0,1.0,0.0), roadway, 10.0), VehicleDef(AgentClass.CAR, 5.0, 2.0), 2),
            Vehicle(VehicleState(VecSE2(  0.0,1.0,0.5), roadway, 10.0), VehicleDef(AgentClass.CAR, 5.0, 2.0), 3),
        ]))

    @test convert(Float64, get(IS_COLLIDING, rec, roadway, 1)) == 1.0
    @test convert(Float64, get(IS_COLLIDING, rec, roadway, 2)) == 0.0
    @test convert(Float64, get(IS_COLLIDING, rec, roadway, 3)) == 1.0

    @test get(LatLonAccel, rec, roadway, 1) == LatLonAccel(0.0, 0.0)
    @test get(AccelTurnrate, rec, roadway, 1) == AccelTurnrate(0.0, 0.0)
    # super-simple test
    for f in allfeatures()
        get(f, rec, roadway, 3)
    end
end

struct DumbFeatureExtractor <: AbstractFeatureExtractor end
@testset begin "features extractor"
    roadway = gen_straight_roadway(3, 1000.0, lane_width=1.0)
    rec = SceneRecord(1, 0.1, 5)
    update!(rec, Scene([
            Vehicle(VehicleState(VecSE2( 0.0,0.0,0.0), roadway, 10.0), VehicleDef(AgentClass.CAR, 5.0, 2.0), 1),
            Vehicle(VehicleState(VecSE2(10.0,0.0,0.0), roadway, 10.0), VehicleDef(AgentClass.CAR, 5.0, 2.0), 2),
        ]))

    ext = DumbFeatureExtractor()
    @test rec_length(ext) == 1
    @test_throws Exception length(ext)
    @test_throws Exception pull_features!(ext, Array{Float64}(undef, 0), rec, roadway, 1)

    ext = FeatureExtractor(AbstractFeature[SPEED, POSFT], 2)
    @test rec_length(ext) == 2
    @test length(ext) == 2
    @test pull_features!(ext, [NaN,NaN], rec, roadway, 1) == [10.0, 0.0]

    subset_ext = SubsetExtractor(ext, [1])
    @test rec_length(subset_ext) == 2
    @test length(subset_ext) == 1
    @test pull_features!(subset_ext, [NaN], rec, roadway, 1) == [10.0]

    subset_ext.subset[1] = 2
    @test pull_features!(subset_ext, [NaN], rec, roadway, 1) == [0.0]

    stan_ext = StandardizingExtractor(ext, [1.0, 2.0], [1.0, 4.0])
    @test rec_length(stan_ext) == 2
    @test length(stan_ext) == 2
    @test pull_features!(stan_ext, [NaN,NaN], rec, roadway, 1) == [(10.0-1.0)/1.0, (0.0-2.0)/4.0]
end

## Lidar sensor tests
#=
What needs to be checked:
1. that vehicles at various angles are detected
2. that the correct distance is returned
3. that when a ray encounters multiple vehicles, the correct vehicle is identified
4. that the range rate is correctly returned
5. that the max range of the ray is respected
=#

@testset "lidar sensor" begin
    #=
    Create a scene with 7 vehicles, an ego in the middle and 6 surrounding
    vehicles. Use this setup to check that lidar sensors work for various
    range and angle settings.
    =#
    num_veh = 7
    ego_index = 1
    roadway = gen_straight_roadway(4, 400.)
    scene = Scene(num_veh)
    # order: ego, fore, rear, left, right, fore_fore, fore_fore_fore
    speeds = [10., 15., 15., 0., -5, 20., 20.]
    positions = [200., 220., 150., 200., 200., 240., 350.]
    lanes = [2,2,2,4,1,2,2]
    for i in 1:num_veh
        lane = roadway.segments[1].lanes[lanes[i]]
        road_idx = RoadIndex(proj(VecSE2(0.0, 0.0, 0.0), lane, roadway))
        veh_state = VehicleState(Frenet(road_idx, roadway), roadway, speeds[i])
        veh_state = move_along(veh_state, roadway, positions[i])
        veh_def = VehicleDef(AgentClass.CAR, 2., 2.)
        push!(scene, Vehicle(veh_state, veh_def, i))
    end

    # basic lidar with sufficient range for all vehicles
    nbeams = 4
    lidar = LidarSensor(nbeams, max_range=200., angle_offset=0.)
    observe!(lidar, scene, roadway, ego_index)
    # order: right, fore, left, back
    @test isapprox(lidar.ranges[1], 2.0, atol=4)
    @test isapprox(lidar.ranges[2], 19.0, atol=4)
    @test isapprox(lidar.ranges[3], 5.0, atol=4)
    @test isapprox(lidar.ranges[4], 49.0, atol=4)
    @test isapprox(lidar.range_rates[1], 0.0, atol=4)
    @test isapprox(lidar.range_rates[2], 5.0, atol=4)
    @test isapprox(lidar.range_rates[3], 0.0, atol=4)
    @test isapprox(lidar.range_rates[4], -5.0, atol=4)

    # angles are such that all the vehicles are missed
    nbeams = 4
    lidar = LidarSensor(nbeams, max_range=200., angle_offset=π/4)
    observe!(lidar, scene, roadway, ego_index)
    # order: right, fore, left, back
    @test isapprox(lidar.ranges[1], 200.0, atol=4)
    @test isapprox(lidar.ranges[2], 200.0, atol=4)
    @test isapprox(lidar.ranges[3], 200.0, atol=4)
    @test isapprox(lidar.ranges[4], 200.0, atol=4)
    @test isapprox(lidar.range_rates[1], 0.0, atol=4)
    @test isapprox(lidar.range_rates[2], 0.0, atol=4)
    @test isapprox(lidar.range_rates[3], 0.0, atol=4)
    @test isapprox(lidar.range_rates[4], 0.0, atol=4)

    # range short enough that it misses the rear vehicle
    nbeams = 4
    lidar = LidarSensor(nbeams, max_range=20., angle_offset=0.)
    observe!(lidar, scene, roadway, ego_index)
    # order: right, fore, left, back
    @test isapprox(lidar.ranges[1], 2.0, atol=4)
    @test isapprox(lidar.ranges[2], 19.0, atol=4)
    @test isapprox(lidar.ranges[3], 5.0, atol=4)
    @test isapprox(lidar.ranges[4], 20.0, atol=4)
    @test isapprox(lidar.range_rates[1], 0.0, atol=4)
    @test isapprox(lidar.range_rates[2], 5.0, atol=4)
    @test isapprox(lidar.range_rates[3], 0.0, atol=4)
    @test isapprox(lidar.range_rates[4], 0.0, atol=4)

end


text:
# deactivate plot GUI, which is not available in Docker
ENV["GKSwstype"] = "100"

# instantiate project
import Pkg
Pkg.activate(@__DIR__)
Pkg.instantiate()

const TEST_LONG = true   # if true, the longer test suite is run; may take > 1 hour
                         # and requires at least 16gb RAM
const TARGET_FOLDER = "result"
const RESULTS_FILE = "results.csv"

function main()
    if !isdir(TARGET_FOLDER)
        mkdir(TARGET_FOLDER)
    end
    global io = open(joinpath(TARGET_FOLDER, RESULTS_FILE), "w")

    println("Running AFF benchmarks...")

    # Heat 3D benchmark
    println("###\nRunning Heat 3D benchmark\n###")
    include("models/Heat3D/heat3d_benchmark.jl")

    # ISS benchmark
    println("###\nRunning ISS benchmark\n###")
    include("models/ISS/ISS_benchmark.jl")

    # Gearbox benchmark
    println("###\nRunning Gearbox benchmark\n###")
    include("models/Gearbox/gearbox_benchmark.jl")

    # Spacecraft benchmark
    if TEST_LONG
        println("###\nRunning spacecraft benchmark\n###")
        include("models/Spacecraft/Spacecraft_benchmark.jl")
    end

    # Building benchmark
    println("###\nRunning building benchmark\n###")
    include("models/Building/building_benchmark.jl")

    # Platoon benchmark
    println("###\nRunning platoon benchmark\n###")
    include("models/Platoon/Platoon_benchmark.jl")

    # Electromechanic break benchmark
    if TEST_LONG
        println("###\nRunning electromechanic-break benchmark\n###")
        include("models/EMBrake/embrake_benchmark.jl")
    end

    # Powertrain benchmark
    println("###\nRunning Powertrain benchmark\n###")
    include("models/Powertrain/Powertrain_benchmark.jl")

    # Clamped benchmark
    println("###\nRunning Clamped benchmark\n###")
    include("models/Clamped/clamped_benchmark.jl")

    print(io, "\n")
    println("Finished running benchmarks.")
    close(io)
    nothing
end

main()


text:
@par function (F::QfForcing)(s::@par(AbstractSimulation))

    cPlaneNumber = 1
    dt = get_dt(s)

    _kf = getKf(F)
    maxdk = getmaxdk(F)
    avgWaveNumInShell2d = getavgk(F)
    nShells2d = length(avgWaveNumInShell2d)

    R = F.R

    s1 = fill!(F.forcex,0)
    s2 = fill!(F.forcey,0)

    u1 = s.u.c.x
    u2 = s.u.c.y

    Ef = F.Ef
    # // get current spectrum for k=1 plane
    calculate_u1u2_mag_spectrum!(Ef,s.u, cPlaneNumber)

  #  # -------------------- Horizontal Forcing --------------------
    ff = 0.0
    fv = 0.0

    powerFraction = 0.99
    @inbounds for j=YRANGE
        conjFactX=1.0
        for i=XRANGE
            k = fsqrt(muladd(KX[i],KX[i],KY[j]^2))
            n = round(Int,k/maxdk) + 1
            if (1 < n <= nShells2d)
                if avgWaveNumInShell2d[n] <= _kf
                    umag = fsqrt(abs2(u1[i,j,1]) + abs2(u2[i,j,1]))
                    Pp = R[n]
                    #fkn = (fsqrt(Ef[n]^2 + 2*dt*Pp*powerFraction) - Ef[n])
                    fkn = dt*Pp*powerFraction/Ef[n]
                    s1[i,j,1] = (u1[i,j,1]/umag)*fkn
                    s2[i,j,1] = (u2[i,j,1]/umag)*fkn
                    #ff += (abs2(s1[i,j,1]) + abs2(s2[i,j,1]))*conjFactX
                    fv += (proj(u1[i,j,1],s1[i,j,1]/dt) + proj(u2[i,j,1],s2[i,j,1]/dt))*conjFactX
                end
            end
            conjFactX=2.0
        end
    end
    ff = ff*0.5/dt
    fv = fv
    horizontalPower = ff+fv
    F.hp[] = horizontalPower

    # # -------------------- Vertical Forcing --------------------
    # /*
    #   Force the points with kx=ky=0 and kz=1,2,3,4,5 so that the power
    #   added is (1-powerFraction) of that added to the horizontal */
    ff = 0.0
    fv = 0.0
    umagsq = 0.0
    #     // get the sum of the existing velocites
    @inbounds for k = 4:6
        umagsq += abs2(u1[1,1,k]) + abs2(u2[1,1,k]) + abs2(u1[1,1,NZ-k+2]) + abs2(u2[1,1,NZ-k+2])
    end

    #f = dt*sqrt(2.0*horizontalPower*(1.0-powerFraction))/umag
    f = dt*horizontalPower*(1-powerFraction)/(umagsq*powerFraction)

    @inbounds for k=4:6
        s1[1,1,k] = u1[1,1,k]*f
        s2[1,1,k] = u2[1,1,k]*f
        s1[1,1,NZ-k+2] = u1[1,1,NZ-k+2]*f
        s2[1,1,NZ-k+2] = u2[1,1,NZ-k+2]*f
        #ff += abs2(s1[1,1,k]) + abs2(s2[1,1,k]) + abs2(s1[1,1,NZ-k+2]) + abs2(s2[1,1,NZ-k+2])
        fv += proj(u1[1,1,k],s1[1,1,k]/dt) + proj(u2[1,1,k],s2[1,1,k]/dt) + proj(u1[1,1,NZ-k+2],s1[1,1,NZ-k+2]/dt) + proj(u2[1,1,NZ-k+2],s2[1,1,NZ-k+2]/dt)
    end
    ff = ff*0.5/dt
    fv = fv
    verticalPower = ff+fv
    F.vp[] = verticalPower
    return 1/dt
end


text:
module DEDataArrays

import StaticArrays, RecursiveArrayTools
import ArrayInterface
import SciMLBase

using LinearAlgebra

using DocStringExtensions

"""
$(TYPEDEF)
"""
abstract type DEDataArray{T,N} <: AbstractArray{T,N} end
const DEDataVector{T} = DEDataArray{T,1}
const DEDataMatrix{T} = DEDataArray{T,2}

# iteration
Base.iterate(A::DEDataArray) = iterate(A.x)
Base.iterate(A::DEDataArray, state) = iterate(A.x, state)

# size
Base.length(A::DEDataArray) = length(A.x)
Base.size(A::DEDataArray) = size(A.x)

# indexing
@inline function Base.getindex(A::DEDataArray, I...)
    @boundscheck checkbounds(A.x, I...)
    @inbounds return A.x[I...]
end
@inline function Base.setindex!(A::DEDataArray, x, I...)
    @boundscheck checkbounds(A.x, I...)
    @inbounds A.x[I...] = x
end
Base.axes(A::DEDataArray) = axes(A.x)
Base.LinearIndices(A::DEDataArray) = LinearIndices(A.x)
Base.IndexStyle(::Type{<:DEDataArray}) = Base.IndexLinear()

Base.copy(A::DEDataArray) = deepcopy(A)

# zero data arrays
@generated function Base.zero(A::DEDataArray)
    assignments = [s == :x ? :(zero(A.x)) :
                   (sq = Meta.quot(s); :(deepcopy(getfield(A, $sq))))
                   for s in fieldnames(A)]
    :(SciMLBase.parameterless_type(A)($(assignments...)))
end

# similar data arrays
@generated function Base.similar(A::DEDataArray, ::Type{T}, dims::NTuple{N,Int}) where {T,N}
    assignments = [s == :x ? :(typeof(A.x) <: StaticArrays.StaticArray ? similar(A.x, T, StaticArrays.Size(A.x)) : similar(A.x, T, dims)) :
                   (sq = Meta.quot(s); :(deepcopy(getfield(A, $sq))))
                   for s in fieldnames(A)]
    :(SciMLBase.parameterless_type(A)($(assignments...)))
end

@generated function Base.similar(A::DEDataArray, ::Type{T}) where {T}
    assignments = [s == :x ? :(typeof(A.x) <: StaticArrays.StaticArray ? similar(A.x, T, StaticArrays.Size(A.x)) : similar(A.x, T)) :
                   (sq = Meta.quot(s); :(deepcopy(getfield(A, $sq))))
                   for s in fieldnames(A)]
    :(SciMLBase.parameterless_type(A)($(assignments...)))
end

@generated function Base.similar(A::DEDataArray) where {T}
    assignments = [s == :x ? :(typeof(A.x) <: StaticArrays.StaticArray ? similar(A.x) : similar(A.x)) :
                   (sq = Meta.quot(s); :(deepcopy(getfield(A, $sq))))
                   for s in fieldnames(A)]
    :(SciMLBase.parameterless_type(A)($(assignments...)))
end

"""
    recursivecopy!(dest::T, src::T) where {T<:DEDataArray}

Recursively copy fields of `src` to `dest`.
"""
@generated function RecursiveArrayTools.recursivecopy!(dest::T, src::T) where {T<:DEDataArray}
    fields = fieldnames(src)

    expressions = Vector{Expr}(undef, length(fields))

    @inbounds for i = 1:length(fields)
        f  = fields[i]
        Tf = src.types[i]
        qf = Meta.quot(f)

        if !ArrayInterface.ismutable(Tf)
            expressions[i] = :( dest.$f = getfield( src, $qf ) )
        elseif Tf <: AbstractArray
            expressions[i] = :( RecursiveArrayTools.recursivecopy!(dest.$f, getfield( src, $qf ) ) )
        else
            expressions[i] = :( dest.$f = deepcopy( getfield( src, $qf ) ) )
        end
    end

    :($(expressions...); dest)
end

"""
    copy_fields(arr:AbstractArray, template::DEDataArray)

Create `DEDataArray` that wraps `arr` with all other fields set to a deep copy of the
value in `template`.

copy_fields!(dest::T, src::T2) where {T<:DEDataArray,T2<:DEDataArray}

Replace all fields of `dest` except of its wrapped array with a copy of the
value in `src`. Arrays are recursively copied.
"""
@generated function copy_fields(arr::AbstractArray, template::DEDataArray)
    assignments = [s == :x ? :(arr) :
                   (sq = Meta.quot(s); :(deepcopy(getfield(template, $sq))))
                   for s in fieldnames(template)]
    :(SciMLBase.parameterless_type(template)($(assignments...)))
end

@generated function copy_fields!(dest::T, src::T2) where
    {T<:DEDataArray,T2<:DEDataArray}

    fields = fieldnames(src)

    expressions = Vector{Expr}(undef, length(fields))

    @inbounds for i = 1:length(fields)
        f  = fields[i]
        Tf = src.types[i]
        qf = Meta.quot(f)

        if f == :x
            expressions[i] = :( )
        elseif !ArrayInterface.ismutable(Tf)
            expressions[i] = :( dest.$f = getfield( src, $qf ) )
        elseif Tf <: AbstractArray
            expressions[i] = :( RecursiveArrayTools.recursivecopy!(dest.$f, getfield( src, $qf ) ) )
        else
            expressions[i] = :( dest.$f = deepcopy( getfield( src, $qf ) ) )
        end
    end

    :($(expressions...); dest)
end

################# Overloads for stiff solvers ##################################

LinearAlgebra.ldiv!(A::DEDataArray,F::Factorization, B::DEDataArray) = ldiv!(A.x,F,B.x)
LinearAlgebra.ldiv!(F::Factorization, B::DEDataArray) = ldiv!(F, B.x)
LinearAlgebra.ldiv!(F::Factorization,A::Base.ReshapedArray{T1,T2,T3,T4}) where {T1,T2,T3<:DEDataArray,T4} = ldiv!(F,vec(A.parent.x))
Base.:+(::LinearAlgebra.UniformScaling,x::DEDataArray) = copy_fields(I + x.x,x)

Base.unsafe_convert(::Type{Ptr{T}}, a::DEDataArray{T}) where {T} = Base.unsafe_convert(Ptr{T}, getfield(a,:x))
ArrayInterface.zeromatrix(x::DEDataArray) = ArrayInterface.zeromatrix(x.x)

################# Broadcast ####################################################

const DEDataArrayStyle = Broadcast.ArrayStyle{DEDataArray}
Base.BroadcastStyle(::Type{<:DEDataArray}) = Broadcast.ArrayStyle{DEDataArray}()
Base.BroadcastStyle(::Broadcast.ArrayStyle{DEDataArray},::Broadcast.ArrayStyle) = Broadcast.ArrayStyle{DEDataArray}()
Base.BroadcastStyle(::Broadcast.ArrayStyle,::Broadcast.ArrayStyle{DEDataArray}) = Broadcast.ArrayStyle{DEDataArray}()
Base.similar(bc::Broadcast.Broadcasted{Broadcast.ArrayStyle{DEDataArray}},::Type{ElType}) where ElType = similar(find_dedata(bc))

find_dedata(bc::Base.Broadcast.Broadcasted) = find_dedata(bc.args)
function find_dedata(args::Tuple)
  !isempty(args) && find_dedata(find_dedata(args[1]), Base.tail(args))
end
find_dedata(x) = x
find_dedata(a::DEDataArray, rest) = a
find_dedata(::Any, rest) = find_dedata(rest)

@inline function Base.copy(bc::Broadcast.Broadcasted{DEDataArrayStyle})
    out = find_dedata(bc)
    copy_fields(copy(unpack(bc)), out)
end

# drop DEData part
@inline unpack(bc::Broadcast.Broadcasted{Style}) where Style = Broadcast.Broadcasted{Style}(bc.f, unpack_args(bc.args))
@inline unpack(bc::Broadcast.Broadcasted{DEDataArrayStyle}) = Broadcast.Broadcasted(bc.f, unpack_args(bc.args))
unpack(x) = x
unpack(x::DEDataArray) = x.x

@inline unpack_args(args::Tuple) = (unpack(args[1]), unpack_args(Base.tail(args))...)
unpack_args(args::Tuple{Any}) = (unpack(args[1]),)
unpack_args(::Any, args::Tuple{}) = ()

# Broadcasting checks for aliasing with Base.dataids but the fallback
# for AbstractArrays is very slow. Instead, we just call dataids on the
# wrapped buffer
Base.dataids(A::DEDataArray) = Base.dataids(A.x)

export DEDataArray, DEDataVector, DEDataMatrix

end


text:
@everywhere Random.seed!(42)
weights = rand(400)
weights = weights/sum(weights)
test_sys_resample = SMC.resample(weights, method = :systematic)
test_multi_resample = SMC.resample(weights, method = :multinomial)
test_poly_resample = SMC.resample(weights, method = :polyalgo)

saved_sys_resample = load("$path/../../reference/resample.jld2", "sys")
saved_multi_resample = load("$path/../../reference/resample.jld2", "multi")
saved_poly_resample = load("$path/../../reference/resample.jld2", "poly")

####################################################################
@testset "Resampling methods" begin
    @test test_sys_resample == saved_sys_resample
    @test test_multi_resample == saved_multi_resample
    @test test_poly_resample == saved_poly_resample
end


text:
# Test file for LargVis
# Sample test:
@assert 1 == 1


text:
"""
    isdependent returns true if eq is dependent on x
"""
isdependent(eq, x) = !isequal(expand_derivatives(Differential(x)(eq)), 0)

"""
    is_number(x) returns true if x is a concrete numerical type
"""
is_number(x::T) where T<:Integer = true
is_number(x::T) where T<:Real = true
is_number(x::T) where T<:Complex = true
is_number(x::T) where T<:Rational = true
is_number(x) = false

is_proper(x) = is_number(x) && !isnan(x) && !isinf(x)

"""
    check_poly checks if a Symbolic expression is a polynomial
"""
# function check_poly(eq, x)
#     if deg(eq, x) == 0 return :not_poly end
#     coefs = values(collect_powers(value(eq), x))
#     if any(z -> isdependent(restore_x(z,x),x), coefs) return :not_poly end
#     if any(z -> z isa Complex, coefs) return :complex_poly end
#     return :real_poly
# end

# is_linear_poly(eq, x) = check_poly(eq, x) != :not_poly && deg(eq, x) == 1

function leading(eq, x)
    l = 0
    k₀ = -1

    for (k,coef) in collect_powers(value(eq), x)
        if !is_number(coef) return 0 end
        if k > k₀
            k₀ = k
            l = coef
        end
    end
    l
end

"""
    deg(p) returns the degree of p if p is a polynomial
"""
deg(p::SymbolicUtils.Add, x) = maximum(deg(t,x) for t in arguments(p))
deg(p::SymbolicUtils.Mul, x) = sum(deg(t,x) for t in arguments(p); init=0)

function deg(p::SymbolicUtils.Pow, x)
    if !isequal(arguments(p)[1], x) return 0 end
    return arguments(p)[2]
end

deg(p::SymbolicUtils.Sym, x) = isequal(p, x) ? 1 : 0
deg(p::SymbolicUtils.Term, x) = 0
deg(p, x) = 0

"""
    var(p) returns the unique variable of an expression (is exists)
"""
function var(p)
    vars = get_variables(p)
    if length(vars) == 1
        return vars[1]
    elseif length(vars) == 0
        return nothing
    else
        error("expects only one variable")
    end
end

###############################################################################

# pox (power-of-x) is a symbolic function to keep track of the powers of x
# pox(k,n) means k*x^n
@syms pox(k, n)

is_pox(x) = istree(x) && operation(x)==pox
is_not_pox(x) = !is_pox(x)

get_coef(p) = is_pox(p) ? arguments(p)[1] : p
get_power(p) = is_pox(p) ? arguments(p)[2] : 0

replace_x(eq, x) = substitute(eq, Dict(x => pox(1,1)))

function restore_x(eq, x)
    r = @rule pox(~k, ~n) => ~k*x^~n
    Prewalk(PassThrough(r))(eq)
end

iscomplex(x) = x isa Complex

count_rule1 = @rule ^(pox(~k, ~n1), ~n2) => isequal(~k,1) ? pox(1, ~n1 * ~n2) : pox(^(~k,~n2), ~n1 * ~n2)
count_rule2 = @rule pox(~k1, ~n1) * pox(~k2, ~n2) => pox(~k1 * ~k2, ~n1 + ~n2)
count_rule3 = @acrule pox(~k, ~n) * ~u::is_not_pox => pox(~k * ~u, ~n)

"""
    collect_powers separates the powers of x in eq (a polynomial) and returns
    a dictionary of power => term
"""
function collect_powers(eq, x)
    eq = expand(expand_derivatives(eq))
    eq = replace_x(eq, x)
    #eq = Prewalk(PassThrough(count_rule1))(eq)
    eq = Fixpoint(Prewalk(PassThrough(Chain([count_rule1, count_rule2, count_rule3]))))(eq)

    if !istree(eq)
        return Dict{Any, Any}(0 => eq)
    elseif is_pox(eq)
        return Dict{Any, Any}(get_power(eq) => get_coef(eq))
    else
        eqs = Dict{Any, Any}()
        for term in arguments(eq)
            n = get_power(term)
            if haskey(eqs, n)
                eqs[n] = eqs[n] + get_coef(term)
            else
                eqs[n] = get_coef(term)
            end
        end

        return eqs
    end
end

###############################################################################

function fibonacci_spiral(n, i)
    ϕ = (1 + sqrt(5)) / 2
    θ, r = mod((i-1)/ϕ, 1), (i-1)/n
    sqrt(r)*cis(2π*θ)
end


text:
# Defines different storages for zarr arrays. Currently only regular files (DirectoryStore)
# and Dictionaries are supported
abstract type AbstractStore end

#Define the interface
"""
  storagesize(d::AbstractStore)

This function shall return the size of all data files in a store.
"""
function storagesize end

"""
    zname(d::AbstractStore)

Returns the name of the current variable.
"""
function zname end

"""
    Base.getindex(d::AbstractStore,i::String)

Returns the data stored in the given key as a Vector{UInt8}
"""
Base.getindex(d::AbstractStore,i::AbstractString) = error("getindex not implemented for store $(typeof(d))")

"""
    Base.setindex!(d::AbstractStore,v,i::String)

Writes the values in v to the given store and key.
"""
Base.setindex!(d::AbstractStore,v,i::AbstractString) = error("setindex not implemented for store $(typeof(d))")

"""
    subdirs(d::AbstractStore)

Returns a list of keys for children stores in the given store.
"""
function subdirs end

"""
    Base.keys(d::AbstractStore)

Returns the keys of files in the given store.
"""
Base.keys(d::AbstractStore) = error("keys function not implemented for store $(typeof(d))")

"""
    newsub(d::AbstractStore, name::String)

Create a new Store as a child of the given store `d` with given `name`. Returns the new created substore.
"""
function newsub end

"""
    getsub(d::AbstractStore, name::String)

Returns the child store of name `name`.
"""
function getsub end

"""
    Base.delete!(d::AbstractStore, k::String)

Deletes the given key from the store.
"""

citostring(i::CartesianIndex) = join(reverse((i - one(i)).I), '.')

Base.getindex(s::AbstractStore, i::CartesianIndex) = s[citostring(i)]
Base.delete!(s::AbstractStore, i::CartesianIndex) = delete!(s, citostring(i))

maybecopy(x) = copy(x)
maybecopy(x::String) = x

function getattrs(s::AbstractStore)
  atts = s[".zattrs"]
  if atts === nothing
    Dict()
  else
    JSON.parse(replace(String(maybecopy(atts)),": NaN,"=>": \"NaN\","))
  end
end
function writeattrs(s::AbstractStore, att::Dict)
  b = IOBuffer()
  JSON.print(b,att)
  s[".zattrs"] = take!(b)
  att
end

is_zgroup(s::AbstractStore) = isinitialized(s,".zgroup")
is_zarray(s::AbstractStore) = isinitialized(s,".zarray")

isinitialized(s::AbstractStore, i::CartesianIndex)=isinitialized(s,citostring(i))

isinitialized(s::AbstractStore, i) = s[i] !== nothing

getmetadata(s::AbstractStore) = Metadata(String(s[".zarray"]))
function writemetadata(s::AbstractStore, m::Metadata)
  met = IOBuffer()
  JSON.print(met,m)
  s[".zarray"] = take!(met)
  m
end

Base.setindex!(s::AbstractStore,v,i::CartesianIndex) = s[citostring(i)]=v

Base.isempty(s::AbstractStore) = isempty(keys(s)) && isempty(subdirs(s))

include("directorystore.jl")
include("dictstore.jl")
include("s3store.jl")
include("consolidated.jl")
include("http.jl")


text:
module Dieter

# %%
# using Reexport
using DataFrames
using DataFramesMeta
using CSV
using SQLite
# using Tables
# @reexport using DataFrames, CSV

# using Feather
# using Arrow  # review whether Arrow needs to be a dependency
# # using Missings

using JuMP
# import MathOptInterface
# const MOI = MathOptInterface
# @reexport using JuMP
using ProgressMeter
# using Plots
# @reexport using Plots
# using ColorSchemes
# @reexport using ColorSchemes

# %% Fundamental structures
include("util.jl")
include("struct.jl")
include("base.jl")

# include("spacetime.jl")

# %% Model components
include("ev.jl")
include("heat.jl")
include("p2g.jl")

# %%
const hoursInYear = 8760     # Units: hours/yr
const H2_energy_density = 39405.6 # Units : MWh / kg
const infeas_cost = 1000     # Currency (cost)

include("model.jl")

# %%
export DieterModel,
InitialiseDieterModel,
create_relation,
initialise_data_file_dict!,
check_files_exist,
parse_data_to_model!,
build_model!,
solve_model!,
generate_results!

# The following code is taken from PowerModels.jl:
# the follow items are also exported for user-friendliness when calling
# `using Dieter`
# so that users do not need to import JuMP to use a solver via `with_optimizer`
import JuMP: with_optimizer
export with_optimizer

# import MathOptInterface: TerminationStatusCode
import JuMP.MOI: TerminationStatusCode
export TerminationStatusCode

# import MathOptInterface: ResultStatusCode
import JuMP.MOI: ResultStatusCode
export ResultStatusCode

for status_code_enum in [TerminationStatusCode, ResultStatusCode]
    for status_code in instances(status_code_enum)
        @eval import MathOptInterface: $(Symbol(status_code))
        @eval export $(Symbol(status_code))
    end
end

end # module Dieter


text:
function four04()
    html("""

    <image src = "images/cover.png"></image>

 It looks like this directory could not be found.

    Thank you for visiting here at Empea. Thank you for sharing the world. We are
    sorry you are lost, or we made a mistake. Would you like to report this issue?

    """)
    Button()
end

end # module


text:
##########
## join ##
##########

function joinSortedIdx_inner(td1::AbstractTimedata,
                             td2::AbstractTimedata)
    ## only works for sorted and increasing index!

    if typeof(idx(td1)) != typeof(idx(td2))
        error("both indices must be of same type to join")
    end

    ## find common dates
    equalIds = findCommonDates(idx(td1), idx(td2))

    df = [td1.vals[equalIds[:, 1], :] td2.vals[equalIds[:, 2], :]]
    td = Timedata(df, idx(td1)[equalIds[:, 1]])
end

function joinSortedIdx_outer(td1::AbstractTimedata,
                             td2::AbstractTimedata)
    ## only works for sorted and increasing index!

    if typeof(idx(td1)) != typeof(idx(td2))
        error("both indices must be of same type to join")
    end

    equalIds = findAllOccurringDates(idx(td1), idx(td2))
    nObs = size(equalIds, 1)
    nVars1 = size(td1, 2)
    nVars2 = size(td2, 2)

    ## init output DataFrame with NAs
    df = DataFrame([eltypes(td1.vals), eltypes(td2.vals)],
                   [names(td1.vals), names(td2.vals)],
                   nObs)
    idxs = Array(eltype(idx(td1)), nObs)

    for ii=1:nObs
        if isna(equalIds[ii, 1]) & !isna(equalIds[ii, 2])
            rowInd = equalIds[ii, 2]
            idxs[ii] = idx(td2)[rowInd]
            for jj=1:nVars1
                df[ii, jj] = NA
            end
            for jj=(nVars1+1):(nVars1+nVars2)
                df[ii, jj] = td2.vals[rowInd, (jj-nVars1)]
            end
        elseif !isna(equalIds[ii, 1]) & isna(equalIds[ii, 2])
            rowInd = equalIds[ii, 1]
            idxs[ii] = idx(td1)[rowInd]
            for jj=1:nVars1
                df[ii, jj] = td1.vals[rowInd, jj]
            end
            for jj=(nVars1+1):(nVars1+nVars2)
                df[ii, jj] = NA
            end
        elseif !isna(equalIds[ii, 1]) & !isna(equalIds[ii, 2])
            rowInd1 = equalIds[ii, 1]
            rowInd2 = equalIds[ii, 2]
            idxs[ii] = idx(td1)[rowInd1]
            for jj=1:nVars1
                df[ii, jj] = td1.vals[rowInd1, jj]
            end
            for jj=(nVars1+1):(nVars1+nVars2)
                df[ii, jj] = td2.vals[rowInd2, (jj-nVars1)]
            end
        end
    end

    td = Timedata(df, idxs)
end

function joinSortedIdx_right(td1::AbstractTimedata,
                             td2::AbstractTimedata)

    ## only works for sorted and increasing index!
    td = joinSortedIdx_left(td2, td1)
    nVars2 = size(td2, 2)
    td2 = [td[:, (nVars2+1):end] td[:, 1:nVars2]]
end

function joinSortedIdx_left(td1::AbstractTimedata,
                            td2::AbstractTimedata)
    ## only works for sorted and increasing index!

    if typeof(idx(td1)) != typeof(idx(td2))
        error("both indices must be of same type to join")
    end

    nIds1 = size(td1, 1)

    ## find dates of 1 in 2
    equalIds = findDatesOfLeftInRight(idx(td1), idx(td2))

    ## get number of variables
    (nVars1, nVars2) = (size(td1, 2), size(td2, 2))

    ## init output DataFrame with NAs
    df = DataFrame([eltypes(td1.vals), eltypes(td2.vals)],
                   [names(td1.vals), names(td2.vals)],
                   nIds1)

    ## successively fill values: run through rows and columns
    for ii=1:nIds1
        for jj=1:nVars1
            df[ii, jj] = td1.vals[ii, jj] # values of 1 must be filled
                                        # always
        end

        if !isna(equalIds[ii])
            ## append values of 2
            for jj=1:nVars2
                df[ii, (nVars1 + jj)] = td2.vals[equalIds[ii], jj]
            end
        end
    end
    tdJoint = Timedata(df, idx(td1))        # contains NAs
end

##################################
## sorted join helper functions ##
##################################

function findAllOccurringDates(idx1, idx2)
    (nIds1, nIds2) = (length(idx1), length(idx2))

    ## some preallocation
    equalIds = DataArray(Int, max(nIds1, nIds2), 2)
    ii = 1; jj = 1; nDatesFound = 1

    ## find none of both indices is completely used
    while (ii <= nIds1) & (jj <= nIds2)
        ## test for equality
        if idx1[ii] == idx2[jj]
            ## note down respective index
            equalIds = addValueToDataArray!(equalIds, @data([ii jj]),
                                 nDatesFound)
            ii += 1; jj += 1; nDatesFound += 1

            ## if no equality: increase smaller index
        elseif idx1[ii] < idx2[jj]
            equalIds = addValueToDataArray!(equalIds, @data([ii NA]),
                                 nDatesFound)
            ii += 1; nDatesFound += 1
        elseif idx1[ii] > idx2[jj]
            equalIds = addValueToDataArray!(equalIds, @data([NA jj]),
                                nDatesFound)
            jj += 1; nDatesFound += 1

        end
    end

    ## append rest
    if jj > nIds2
        while ii <= nIds1
            equalIds = addValueToDataArray!(equalIds, @data([ii NA]),
                                nDatesFound)
            ii += 1; nDatesFound += 1
        end
    elseif ii > nIds1
        while jj <= nIds2
            equalIds = addValueToDataArray!(equalIds, @data([NA jj]),
                                nDatesFound)
            jj += 1; nDatesFound += 1
        end
    end
    return equalIds
end

function addValueToDataArray!(equalIds, DArow, nDatesFound)
    if nDatesFound > size(equalIds, 1)
        equalIds = vcat(equalIds, DArow)
    else
        equalIds[nDatesFound, :] = DArow
    end
    equalIds
end

function getIndices(equalIds, idx1, idx2)
    nObs = size(equalIds, 1)
    ids = Array(eltype(idx1), nObs)

    for ii=1:nObs
        if !isna(equalIds[ii, 1])
            ids[ii] = idx1[equalIds[ii, 1]]
        else
            ids[ii] = idx2[equalIds[ii, 2]]
        end
    end
    ids
end

function getIndices(equalIds, idx)
    nObs = size(equalIds, 1)
    ids = DataArray(eltype(idx), nObs)

    for ii=1:nObs
        if !isna(equalIds[ii])
            ids[ii] = idx[equalIds[ii]]
        end
    end
    ids
end

function findDatesOfLeftInRight(idx1, idx2)
    ## for each date in idx1, find its entry in idx2. If date does not
    ## exist, mark it with NA.

    nIds1 = length(idx1)
    nIds2 = length(idx2)

    ## running variables
    ii = 1
    jj = 1

    ## preallocation with NAs
    equalIds = DataArray(typeof(ii), nIds1)

    ## find common dates
    while ii <= nIds1
        ## test for equality
        if idx1[ii] == idx2[jj]
            ## note down respective index
            equalIds[ii] = jj
            ii += 1
            jj += 1
            if jj > nIds2
                break
            end
            ## if no equality: increase smaller index
        elseif idx1[ii] < idx2[jj]
            equalIds[ii] = NA
            ii += 1
        elseif idx1[ii] > idx2[jj]
            jj += 1
            if jj > nIds2
                break
            end
        end
    end
    return equalIds
end

function findCommonDates(idx1, idx2)

    nIds1 = length(idx1)
    nIds2 = length(idx2)

    ## running variables
    ii = 1
    jj = 1

    ## preallocation with NAs
    equalIds = [[] []]

    ## find common dates
    while ii <= nIds1
        ## test for equality
        if idx1[ii] == idx2[jj]
            ## note down respective index
            equalIds = vcat(equalIds, [ii jj])
            ii += 1
            jj += 1
            if jj > nIds2
                break
            end
            ## if no equality: increase smaller index
        elseif idx1[ii] < idx2[jj]
            ii += 1
        elseif idx1[ii] > idx2[jj]
            jj += 1
            if jj > nIds2
                break
            end
        end
    end
    return equalIds
end


text:
abstract type SharpRD end

abstract type VarianceEstimator end

struct EickerHuberWhite <: VarianceEstimator end
"""
    Homoskedastic()

Variance estimator for the RD estimator that assumes homoskedasticity.
"""
struct Homoskedastic <: VarianceEstimator end

abstract type BandwidthSelector end

bandwidth(h::Number, args...) = h
_string(::Number) = ""

abstract type BinwidthSelector end

binwidth(b::Number, args...) = b

abstract type WorstCaseCurvatureSelector end

curvature(b::Number, args...) = b


text:
module ExoticDataStructures

    import Base: length, isempty, start, next, done
                 show, dump, empty!, getindex, setindex!, get, get!
                 in, haskey, keys, merge, copy, cat,
                 push!, pop!, shift!, unshift!,
                 union!, delete!, similar, sizehint,
                 isequal, hash,
                 map, reverse

    export RedBlackTree

    include("redblacktree.jl")

end # module


text:
#= module xim_basics.jl
   started at 2020-05-20
   aim to contain support functions to xim package

=#
#module xim_basics

#export xim_basics_testsum
#export xim_basics_load_packages

#*************************************************
function xim_basics_testsum(x,y)
  c=x+y
  c
end
#xim_basics_testsum(a,b)=a+b
#*************************************************
function xim_basics_load_packages()
  #using Pkg;
  #Pkg.add("Images");
  #Pkg.add("ImageView");
  #Pkg.add("TestImages");
  #Pkg.add("FileIO");
  #Pkg.add("Colors");
  #Pkg.add("QuartzImageIO");

  #using Images, ImageView, TestImages, FileIO, Colors, QuartzImageIO;
  println("impossible?")
end
#*************************************************

max_stars_repo_name: NelsonMCLopes/xim

max_issues_repo_name: NelsonMCLopes/xim

max_forks_repo_name: NelsonMCLopes/xim


text:
ansi(x) = string("\x1b[", x, "m")
clear = ansi(0)
bold = ansi(1)
black = ansi(30)
red = ansi(31)
green = ansi(32)
yellow = ansi(33)
blue = ansi(34)
magenta = ansi(35)
cyan = ansi(36)
white = ansi(37)

"""
    log(str...)

    Simply prints the arguments provided. This may be reworked later to provide
    more robust logging capabilities (e.g. to files)
"""
log(str...) = println(str...)

"""
    warn(str...)

    Prints the arguments provided, but with a yellow color and a bold "WARNING"
    in front. This may be reworked later to provide more robust logging
    capabilities (e.g. to files)
"""
warn(str...) = println(bold, yellow, "WARNING: ", clear, yellow, str..., clear)

"""
    error(str...)

    Prints the arguments provided, but with a red color and a bold "ERROR" in
    front. This may be reworked later to provide more robust logging
    capabilities (e.g. to files)
"""
error(str...) = println(bold, red, "ERROR: ", clear, red, str..., clear)


text:
using Dates: Date, UTInstant, Day

struct DateStrategy <: SearchStrategy
    min::Date
    max::Date

    DateStrategy(min::Nothing, max::Nothing) = new(Date(1,1,1), Date(9999,12,31))
    DateStrategy(min::Nothing, max::Date) = new(Date(1,1,1), max)
    DateStrategy(min::Date, max::Nothing) = new(min, Date(9999,12,31))
    function DateStrategy(min::Date, max::Date)
        min <= max || error("DateStrategy has unsatisfiable bounds")
        new(min, max)
    end
end

function draw(strategy::DateStrategy)
    days = rand(strategy.min.instant.periods.value:strategy.max.instant.periods.value)
    Date(UTInstant(Day(days)))
end

dates(; min=nothing, max=nothing) = DateStrategy(min, max)
Base.eltype(::DateStrategy) = Date


text:
#contour qubit

using CairoMakie
let
    function tα_qubit(β, ψ1, ψ2, fα, f)
        2 + 2 * β - cos(ψ1) - cos(ψ2) - 2 * β * cos(π * fα) * cos(2 * π * f + π * fα - ψ1 - ψ2)
    end
    ψ1 = ψ2 = LinRange(0, 4 * π, 100)
    z = [tα_qubit(0.61, x, y, 0.2, 0.1) for x in ψ1, y in ψ2]

    fig = Figure(resolution = (600, 400), backgroundcolor = :honeydew)
    ax = Axis(fig[1, 1], aspect = 1, xlabel = "ψ1", ylabel = "ψ2")
    cls = CairoMakie.contour!(ψ1, ψ2, z; colormap = :tableau_summer, levels = 20, linewidth = 1.5)
    Colorbar(fig[1, 2], cls, label = "α-q")
    limits!(ax, 0, 4π, 0, 4π)
    colsize!(fig.layout, 1, Aspect(1, 1.0))
    display(fig)
end

max_stars_repo_name: Ellon-M/visualizations

max_issues_repo_name: Ellon-M/visualizations

max_forks_repo_name: Ellon-M/visualizations


text:
"""
    metaoperator

Top level operator.
**jula> help("Metaoperator")
"""
function metaoperator(command, logicset::LogicalCombo; verbose = true)
    counter!()

    logicsetcopy = deepcopy(logicset)

    metaset = "XOR|IFF|IF|THEN|AND|OR"

    (sum(logicset[:]) == 0) && return logicset
    !occursin(Regex("(?:\b| )([><=|!+\\-^&]{4}|$metaset)(?:\b| )"), command) &&
      return superoperator(command, logicset, verbose=verbose)

    m = match(Regex("(^.*)(?:\b| )([><=|!+\\-^&]{4}|$metaset)(?:\b| )(.*?\$)"), command)
    left, operator, right = m.captures

    ℧left  = metaoperator(left , logicset, verbose=verbose)[:]
    ℧right = metaoperator(right, logicset, verbose=verbose)[:]
    ℧η = logicsetcopy[:]

    (operator ∈ ["&&&&", "AND"]) && (℧η = logicset[:] .& (℧left .& ℧right))
    (operator ∈ ["====", "IFF"]) && (℧η = logicset[:] .& (℧left .== ℧right))
    (operator ∈ ["===>", "THEN"]) && (℧η[℧left] = logicset[:][℧left]  .& ℧right[℧left])
    (operator ∈ ["!==>"]) && (℧η[.!℧left] = logicset[:][.!℧left]  .& ℧right[.!℧left])
    (operator ∈ ["<===", "IF"]) && (℧η[℧right] = logicset[:][℧right] .& ℧left[℧right])
    (operator ∈ ["||||", "OR"]) && (℧η = logicset[:] .& (℧left .| ℧right))

    (operator ∈ ["^^^^", "XOR", "!==="]) &&
      (℧η = logicset[:] .& ((℧left .& .!℧right) .| (.!℧left .& ℧right)))

    logicsetcopy[:] = ℧η
    logicsetcopy
end


text:
using HDF5

"""Read an IMS measurement file (in all available formats)"""
function imsread(file_name::AbstractString; normalize=true, asfloat=true)
    isfile(file_name) || error(@sprintf("`%s` is not a readable file", file_name))
    ims = ishdf5(file_name) ? imsread_hdf5(file_name) : imsread_csv(file_name)
    if normalize
        if length(ims["R"]) != size(ims["S"], 1)
            ims["S"] = ims["S"]'
        end
        if mean(ims["S"]) < 0
            ims["S"] = -ims["S"]
        end
    end
    if asfloat
        ims["S"] = map(Float64, ims["S"])
    end
    return ims
end

"""CSV parser"""
function imsread_csv(file_name::AbstractString)
    data = readcsv(file_name, Float64, skipstart=132, use_mmap=true)
    S = data[:,3:end]'
    D = data[:,2]
    T = data[:,1]

    tR = split(readline(file_name, 131), ",")
    @assert startswith(tR[1], "\\")
    @asserteq strip(tR[2]) "tR"
    R = float(tR[3:end])
    @asserteq size(S,1) length(R)

    return Dict(@qvs S D T R)
end

"""HDF5 parser"""
function imsread_hdf5(h5_fname::AbstractString)
    h5open(h5_fname, "r") do file
        if HDF5.exists(file, "IMS")
            file = g_open(file, "IMS")
        end
        return Dict(k => read(file, k) for k in names(file))
    end
end


text: println("Hello world.")


text:
using CUDA
const pretrained_folder = "./pretrained_discriminators"

# Dictionary for pretrained discriminator paths
pretrained = Dict(
            "clickbait" => "clickbait_classifier_head_1024.bson",
            "toxicity"  => "toxicity_classifier_head_768.bson",
            "sentiment" => "sentiment_classifier_head_1024.bson"
        )

"""
struct ClassifierHead
    linear_layer::Dense
    embed_size::Int
    class_size::Int
end

Struct for ClassifiedHead, defined with a single linear layer and two paramters: embed_size-> Size of Embedding, class_size->Number of classes.

"""
struct ClassifierHead
    linear_layer::Dense
    embed_size::Int
    class_size::Int
end

"""
    ClassifierHead(class_size::Int=1, embed_size::Int=768; load_from_pretrained=false, discrim=nothing, file_name=nothing, path=nothing)

Function to initiate ClassifierHead layer, with defined class_size and embed_size. If `load_from_pretrained` is set to be true, load model from pretrained models (based on `discrim` name) or from specified `path`.
"""
function ClassifierHead(class_size::Int=1, embed_size::Int=768; load_from_pretrained=false, discrim=nothing, file_name=nothing, path=nothing)
    global pretrained_folder
    if load_from_pretrained==true
        if discrim==nothing
            error("Discriminator not provided")
        else
            if discrim in keys(pretrained)
                path_f = joinpath(@__DIR__, pretrained_folder, pretrained[discrim])
            else
                print("Discriminator $discrim not in pretrained discriminators, checking if `path` provided")
                if isnothing(path)
                    error("`path` not provided for custom discriminator bson file. Try `ClassifierHead(...; path=...)` or `getDiscriminator(...; ... , path=... )`")
                else
                    path_f = path
                end
            end
            BSON.@load path_f config_metadata weights
            layer = Dense(weights["weight"], weights["bias"])
            class_size = length(weights["bias"])
        end
    else
        layer = Dense(embed_size, class_size)
        config_metadata=nothing
    end
    ClassifierHead(layer, embed_size, class_size), config_metadata
end

Flux.@functor ClassifierHead (linear_layer,)

"""
    (m::ClassifierHead)(hidden_state)

Method for passing hidden state through Linear layer of ClassifierHead `m`, return the final logits.

"""
function (m::ClassifierHead)(hidden_state)
    lm_logits = m.linear_layer(hidden_state)
    lm_logits
end

"""
struct DiscriminatorV1
    cl_head::ClassifierHead
    model
    embed_size::Int
end

Struct to contain the model and ClassifierHead of PPLM model, to be used in masked training.

"""
struct DiscriminatorV1
    cl_head::ClassifierHead
    model
    embed_size::Int
end

Flux.@functor DiscriminatorV1 (cl_head,)

"""
    avg_representation(m::DiscriminatorV1, input_ids; args=nothing)

Function to create average representation (without mask) of Input tokens `input_ids`, using the DiscriminatorV1 `m`.

"""
function avg_representation(m::DiscriminatorV1, input_ids; args=nothing)
    if args!=nothing && args.cached == true
        return input_ids
    else
        outputs = m.model(input_ids; output_attentions=false,
                    output_hidden_states=true, use_cache=false)
        hidden = output.hidden_states
        hidden = sum(hidden, dims = 2)
        hidden_ = dropdims(hidden, dims=2)
        return hidden_
    end
end

# DiscriminatorV1 -> without mask on hidden_state
function (m::DiscriminatorV1)(x; args=nothing)
    hidden_state = avg_representation(m, x; args)
    logits = m.cl_head(hidden_state)
    logits, dims=1
end

"""
struct DiscriminatorV2
    cl_head::ClassifierHead
    model
    embed_size::Int
end

Struct to contain the model and ClassifierHead of PPLM model, to be used in masked training.

"""
struct DiscriminatorV2
    cl_head::ClassifierHead
    model
    embed_size::Int
end

Flux.@functor DiscriminatorV2 (cl_head,)

"""
    avg_representation(m::DiscriminatorV2, input_ids, mask; args=nothing)

Function to create masked average representation of Input tokens `input_ids`, using the DiscriminatorV2 `m`.

"""
function avg_representation(m::DiscriminatorV2, input_ids, mask; args=nothing)
    if args!=nothing && args.cached == true
        return input_ids
    else
        outputs = m.model(input_ids; output_attentions=false,
                    output_hidden_states=true, use_cache=false)
        hidden = outputs.hidden_states[end]
        hidden.*= mask
        hidden = sum(hidden, dims = 2) ./ sum(mask, dims=2)
        hidden_ = dropdims(hidden, dims=2)
        return hidden_
    end
end

function (m::DiscriminatorV2)(x, mask; args=nothing)
    hidden_state = avg_representation(m, x, mask; args)
    logits = m.cl_head(hidden_state)
    logits
end

"""
    get_discriminator(model; load_from_pretrained=false, discrim=nothing, file_name=nothing, version=2, class_size::Int=1, embed_size::Int=768, path=nothing)

Function to create discriminator based on provided model. Incase, `load_from_pretrained` is set to be true, loads ClassifierHead layer from pretrained models or `path` provided.

"""
function get_discriminator(model; load_from_pretrained=false, discrim=nothing, file_name=nothing, version=2, class_size::Int=1, embed_size::Int=768, path=nothing)
    global device
    model = model |> device
    cl_head, _ = ClassifierHead(class_size, embed_size; load_from_pretrained=load_from_pretrained, discrim=discrim, file_name=file_name, path=path) |> gpu
    if version == 1
        return DiscriminatorV1(cl_head, model, embed_size)
    else
        return DiscriminatorV2(cl_head, model, embed_size)
    end
end

"""
    save_classifier_head(cl_head; file_name=nothing, path=nothing, args=nothing, register_discrim=true, discrim_name="")

Function to save the ClassifiedHead as a BSON once the training is complete, based on the path provided. In case path is set as nothing, it saves the discriminators in `./pretrained_discriminators` folder relative to the package directory.
"""
function save_classifier_head(cl_head; file_name=nothing, path=nothing, args=nothing, register_discrim=true, discrim_name="")
    if path == nothing
        joinpath(@__DIR__,"./pretrained_discriminators")
    end
    if file_name == nothing
        file_name = string("custom_classifier_head_", cl_head.embed_size, "_.jl")
    end

    cl_head = cl_head |> cpu
    save_path = joinpath(path, file_name)
    config_metadata = args
    weights = Dict("weight"=> cl_head.linear_layer.W, "bias"=>cl_head.linear_layer.bias)
    if args==nothing
        @warn "Saving classifier without Hyperparameter information"
        config_metadata = Dict()
        config_metadata["embed_size"]=size(weights["weight"])[2]
        config_metadata["class_size"]=size(weights["weight"])[1]
        config_metadata["default_class"]=1
    else
        config_metadata = Dict()
        config_metadata=args
    end

    BSON.@save save_path config_metadata weights
end

"""
    save_discriminator(discrim, discrim_name="Custom"; file_name=nothing, path=nothing, args=nothing)

Function to save ClassifiedHead part of discriminator (by calling `save_classifier_head` function), which is the only trainable part of discriminator
"""
function save_discriminator(discrim, discrim_name="Custom"; file_name=nothing, path=nothing, args=nothing)
    if path == nothing
        path = joinpath(@__DIR__, "./pretrained_discriminators")
    end
    if file_name == nothing
        file_name = string("custom_classifier_head_", discrim.embed_size, "_.jl")
    end
    save_path = joinpath(path, file_name)
    println("Saving classifier head weights for the discriminator to $save_path")
    save_classifier_head(discrim.cl_head; file_name=file_name, path=path, args=args, discrim_name=discrim_name)
    print("Discriminator saved successfully")
end


text:
# Spacecraft Models
abstract type Spacecraft end

# Simple Spacecraft Model
struct SimpleSpacecraft <: Spacecraft

    # Initial mass
    initMass::Float64

    # Initial propellant mass
    initProp::Float64

    # Max thrust
    tMax::Float64

    # Specific impulse
    isp::Float64

    # Exaust velocity
    c::Float64

    function SimpleSpacecraft(initMass, initProp, tMax, isp)
        c = 9.81*isp
        new(initMass, initProp, tMax, isp, c)
    end
end


text:
using Plots
import Statistics: median, std
import MyterialColors: black, blue_light
using Term
using Term.progress
import InfiniteOpt: objective_value
install_term_logger()
# install_stacktrace()

using jcontrol
using jcontrol.visuals
import jcontrol: closest_point_idx, Track, State
import jcontrol.comparisons: ComparisonPoints
import jcontrol.bicycle: Bicycle

function compare(;  problemtype=:dynamics)
    # ---------------------------------- run MTM --------------------------------- #
    track = Track(; start_waypoint=4, keep_n_waypoints=-1)

    coptions = ControlOptions(;
    u_bounds=Bounds(10, 80),
    δ_bounds=Bounds(-50, 50, :angle),
    δ̇_bounds=Bounds(-6, 6),
    ω_bounds=Bounds(-400, 400, :angle),
    v_bounds=Bounds(-15, 15),
    Fu_bounds=Bounds(-2000, 3000),
    )

    bike = Bicycle(;
        l_f=3, l_r=2, c=4000, width=2.0
    )

    track, bike, _, solution = run_mtm(
        problemtype,  # model type
        2;  # supports density
        showtrials=nothing,
        control_options=coptions,
        bike=bike,
        track=track,
        n_iter=5000,
        fcond=State(; u=30, ψ=0, n=-.5),
        timed=false,
        showplots=false,
    )
    @info "Duration" solution.t[end]

    # plot model trajectory
    plt = draw(:arena)
    draw!(track)

    # -------------------------- do comparison with data ------------------------- #
    # load data
    trials = load_cached_trials(; keep_n = 100,)
    cpoints = ComparisonPoints(track; δs=5, trials=trials)
    fasttrials = filter(t -> t.duration <= solution.t[end], trials)

    # print percentage of fast trials (roudned)
    @info "Fast trials: $(round(100 * length(fasttrials) / length(trials)))% | ($(length(fasttrials))/$(length(trials)))"

    # show data
    draw!.(trials; lw=2, alpha=.15)
    draw!.(cpoints.points)
    plot_bike_trajectory!(solution, bike; showbike=false)

    scatter!(
        solution.x, solution.y, marker_z=solution.ω, c=:bwr, label=nothing
    )
    # do comparison
    speedplot = plot(; title="speed", xlabel="s (cm)", ylabel="speed cm/s", legend=false)
    uplot = plot(; title="u", xlabel="s (cm)", ylabel="u cm/s", legend=false)
    vplot = plot(; title="v", xlabel="s (cm)", ylabel="v cm/s", legend=false)
    ωplot = plot(; title="ang.vel.", xlabel="s (cm)", ylabel="avel rad/s", legend=false)

    @track for trial in trials[1:100]
        # plot trial kinematics
        plot!(speedplot, trial.s, trial.speed; color=black, alpha=0.7)
        plot!(uplot, trial.s, trial.u; color=black, alpha=0.7)
        plot!(vplot, trial.s, trial.v; color=black, alpha=0.7)
        plot!(ωplot, trial.s, trial.ω; color=black, alpha=0.7)
    end

    # plot model kinematics at CP
    for cp in cpoints.points
        # get the closest model point
        idxs = closest_point_idx(solution.x, cp.x, solution.y, cp.y)

        if problemtype==:kinematics
            speed = solution.u[idxs]
            u = solution.u[idxs] * cos(solution.β[idxs])
            v = solution.u[idxs] * sin(solution.β[idxs])
        else
            u, v = solution.u[idxs], solution.v[idxs]
            speed = sqrt(u^2 + v^2)
        end

        scatter!(speedplot, [cp.s], [speed]; color="red", ms=10)
        scatter!(uplot, [cp.s], [u]; color="red", ms=10)
        scatter!(vplot, [cp.s], [v]; color="red", ms=10)
        scatter!(ωplot, [cp.s], [solution.ω[idxs]]; color="red", ms=10)
    end

    # ------------------------------------ fin ----------------------------------- #
    h = histogram(
        map(t->t.duration, trials), color="black", label=nothing, xlim=[0, 15]
    )
    duration = solution.t[end]

    plot!(h, [duration, duration], [0, 200], lw=5, alpha=.8, color="red", label=nothing)

    l = @layout [
        a{0.5w} grid(5, 1)
    ]
    return display(
        plot(
            plt,
            speedplot,
            uplot,
            vplot,
            ωplot,
            h;
            layout=l,
            size=(1600, 1200),
        ),
    )
end

print("\n\n" * hLine("start"; style="bold green"))

compare()

print("\n\n" * hLine("finish"; style="bold blue"))

max_stars_repo_name: FedeClaudi/LocomotionControl

max_issues_repo_name: FedeClaudi/LocomotionControl

max_forks_repo_name: FedeClaudi/LocomotionControl


text:
# Everying Euler trails/tours

export euler

# This function finds an Euler trail/tour given initial
# vertices. Returns the trail if it exists of [] if it does not.

"""
`euler(G,u,v)` finds an Euler trail in `G` starting at `u` and ending
at `v` returned as a list of vertices (in the order they are visited
on the trail).

`euler(G,u)` finds an Euler tour beginning and ending at
`u`. Alternatively, call `euler(G)` and the initial/final vertex will
be selected for you.

Note: The algorithm will succeed even if there are isolated vertices
in the graph (just don't choose an isolated vertex as the first/last).

If no Euler trail/tour exists, an empty list is returned.
"""
function euler{T}(G::SimpleGraph{T}, u::T, v::T)
    notrail = T[]

    # perform basic checks:
    if ! (has(G,u) && has(G,v))
        error("One or both of these vertices is not in this graph")
    end

    # special case: if all vertices have degree zero
    if all( [deg(G,x)==0 for x in G.V] )
        if u==v
            return [u]
        else
            return notrail
        end
    end

    # if either vertex has degree zero, we're doomed
    if deg(G,u)==0 || deg(G,v)==0
        return notrail
    end

    # vertex degree checks
    if u==v
        for x in G.V
            if deg(G,x)%2 == 1
                return notrail
            end
        end
    else        # u != v
        for x in G.V
            if x==u || x==v
                if deg(G,x)%2==0
                    return notrail
                end
            else
                if deg(G,x)%2==1
                    return notrail
                end
            end
        end
    end

    # Remove isolates and check for connectivity
    GG = trim(G)
    if !is_connected(GG)
        return notrail
    end

    # all tests have been satisfied. Now find the trail in GG using a
    # helper function.
    return euler_work!(GG,u)
end

# special case: find an Euler tour from a specified vertex
function euler{T}(G::SimpleGraph{T},u::T)
    return euler(G,u,u)
end

# special case: find any Euler tour. If the graph is connected, any
# vertex will do but if there are isolated vertices, we don't want to
# pick one of those!
function euler{T}(G::SimpleGraph{T})
    if NV(G)==0
        return T[]
    end

    verts = vlist(G)

    # search for a vertex that isn't isolated
    for u in verts
        if deg(G,u) > 0
            return euler(G,u,u)
        end
    end

    # if we reach here, the graph only has isolated vertics. Let it
    # work from any isolated vertex (and likely fail).
    u = verts[1]
    return euler(G,u,u)

end

# private helper function for euler()
function euler_work!{T}(G::SimpleGraph{T}, u::T)
    trail = T[]

    while true
        # if last vertex
        if NV(G)==1
            append!(trail,[u])
            return trail
        end

        # get the neighbors of u
        Nu = G[u]

        # if only one neighbor delete and move on
        if length(Nu)==1
            w = Nu[1]
            delete!(G,u)
            append!(trail,[u])
            u = w
        else
            for w in Nu
                if ! is_cut_edge(G,u,w)
                    delete!(G,u,w)
                    append!(trail,[u])
                    u = w
                    break
                end
            end
        end
    end
    error("This can't happen")
end


text:
abstract type SGStkeBC end

"""
    NoSGStkeFlux() :: SGStkeBC

No energy flux across boundary
"""
struct NoSGStkeFlux <: SGStkeBC end
function atmos_sgstke_boundary_state!(
    nf,
    bc_sgstke::NoSGStkeFlux,
    atmos,
    _...,
)
    nothing
end
function atmos_sgstke_normal_boundary_flux_second_order!(
    nf,
    bc_sgstke::NoSGStkeFlux,
    atmos,
    _...,
)
    nothing
end

"""
    PrescribedSGStke(sgs) :: SGStkeBC

Dirichlet boundary condition
"""
struct PrescribedSGStke{SGS} <: SGStkeBC
    e_SGS::SGS
end
function atmos_sgstke_boundary_state!(
    nf,
    bc_sgstke::PrescribedSGStke,
    atmos,
    state⁺,
    args,
)
    @unpack aux⁻, state⁻, t = args
    state⁺.sgstke.ρe_SGS = state⁺.ρ * bs_sgstke.e_SGS(state⁻, aux⁻, t)
end
function atmos_sgstke_normal_boundary_flux_second_order!(
    nf,
    bc_sgstke::PrescribedSGStke,
    atmos,
    fluxᵀn,
    args,
)
    @unpack state⁻, aux⁻, diffusive⁻, hyperdiff⁻, t, n⁻ = args
    tend_type = Flux{SecondOrder}()
    _args⁻ = (;
        state = state⁻,
        aux = aux⁻,
        t,
        diffusive = diffusive⁻,
        hyperdiffusive = hyperdiff⁻,
    )
    pargs = merge(_args⁻, (precomputed = precompute(atmos, _args⁻, tend_type),))
    total_flux =
        Σfluxes(SGSTKE(), eq_tends(SGSTKE(), atmos, tend_type), atmos, pargs)
    nd_ρe_SGS = dot(n⁻, total_flux)
    # both sides involve projections of normals, so signs are consistent
    fluxᵀn.sgstke.ρe_SGS += nd_ρe_SGS
end

"""
    PrescribedSGStkeFlux(f_sgs) :: SGStkeBC
Prescribe the net inward SGS tke flux across the boundary by `f_sgs`, a function
with signature `f_sgs(state, aux, t)`, returning the flux (in W/m^2).
"""
struct PrescribedSGStkeFlux{FSGS} <: SGStkeBC
    f_sgs::FSGS
end
function atmos_sgstke_boundary_state!(
    nf,
    bc_sgstke::PrescribedSGStkeFlux,
    atmos,
    _...,
) end
function atmos_sgstke_normal_boundary_flux_second_order!(
    nf,
    bc_sgstke::PrescribedSGStkeFlux,
    atmos,
    fluxᵀn,
    args,
)
    @unpack state⁻, aux⁻, t = args

    # DG normal is defined in the outward direction
    # we want to prescribe the inward flux
    fluxᵀn.sgstke.ρe_SGS -= bc_sgstke.f_sgs(state⁻, aux⁻, t)
end


text:
#PROCESSING FILE
#All the mathematic and processing fonctions
#---------------------------------------------------------------------
# Analysis Parameters :
#	statParam -> Horizontal_sum_of_pixels(Array)|Vertical_sum_of_pixels(Array)
#	IdealParam -> mean of all statparams
#	FactorWeight -> weighting(std(statParam)) -> array 1-D of floats , closer to 1 == highly representative.
#	Similarity -> sum(abs(IdealParam-statParam).*FactorWeight) , closer to 0 is closer to the ideal pattern.
#	CategorizingLimit -> [mean(Similarity)+X*std(Similarity),mean(Similarity)-X*std(Similarity)]
# 		Using the "average" allow to avoid the offset due to the difference induced by the irrealistic shape of the ideal pattern.
#		the X is the amout of learning data included into the accepting range. Working like 1/2/3 sigma in gaussian repartition.

#Validated Footprint Coordinate -> Coordiante that locate a pattern that passed the partial comparaison, the comparaison , the filtering , and the Similarity check.
#Every Validated Footprint Coordinate is compared to other footprint to avoid redundancy
#Consequently, every Validated Footprint Coordinate is linked to the coressponding Similarity value
#When comparing to array of VFC to an another (of an another footprint), if double are noticed, the one with the smaller Similary value will be choosed.


text:
module TestConfoundersPreparation

using Test
using SnpArrays
using TMLEEpistasis
using DataFrames
using CSV

@testset "Various functions" begin
    # Test issnp
    @test TMLEEpistasis.issnp("A") == true
    @test TMLEEpistasis.issnp("AA") == false

    #Test bounds
    @test TMLEEpistasis.bounds(0, -10, 10) == (-10^4, 10^4)
    @test TMLEEpistasis.bounds(15000, 1000, 30000) == (1000, 30000)

    # Test notin_ldblocks
    ldblocks = DataFrame(chr        = [1, 3],
                        lower_bound = [10, 19],
                        upper_bound = [100, 109])
    ldblocks = groupby(ldblocks, :chr)

    snp_info = DataFrame(chromosome=[1, 1, 2, 3, 3],
                        position=[21, 5, 9, 102, 4])

    expected_snps = DataFrame(chromosome = [1, 2, 3],
                              position   = [5, 9, 4])
    @test filter(x -> TMLEEpistasis.notin_ldblocks(x, ldblocks), snp_info) == expected_snps
end

@testset "Test filter_chromosome" begin
    parsed_args = Dict(
        "input"  => SnpArrays.datadir("mouse"),
        "output" => joinpath("data", "filtered-mouse"),
        "qcfile" => joinpath("data", "ukbb", "qcfile.txt"),
        "ld-blocks" => joinpath("data", "VDR_LD_blocks.txt"),
        "maf-threshold" => 0.31
    )
    filter_chromosome(parsed_args)

    filtered = SnpData(parsed_args["output"])

    @test filtered.snp_info.snpid == ["rs13476318"]
    @test size(filtered.snparray) == (1940, 1)

    # Clean
    for ext in [".bed", ".bim", ".fam"]
        rm(parsed_args["output"]*ext)
    end
end

@testset "Test merge_beds" begin
    genotypes_dir = joinpath("data", "ukbb", "genotypes")
    chr_prefix = joinpath(genotypes_dir, "ukbb")

    parsed_args = Dict(
            "input"  => chr_prefix,
            "output" => joinpath(genotypes_dir, "ukbb_merged")
        )
    merge_beds(parsed_args)
    merged = SnpData(parsed_args["output"])

    @test length(unique(merged.snp_info.chromosome)) == 3
    # Clean

    for ext in [".bed", ".bim", ".fam"]
        rm(parsed_args["output"]*ext)
    end

end

@testset "Test adapt_flashpca" begin
    parsed_args = Dict(
        "input" => joinpath("data", "flashpca_output.txt"),
        "output" => "flashpca_after_adapt.csv",
    )
    adapt_flashpca(parsed_args)

    output = CSV.File(parsed_args["output"]) |> DataFrame

    @test names(output) == ["SAMPLE_ID",
                            "PC1",
                            "PC2",
                            "PC3",
                            "PC4",
                            "PC5",
                            "PC6"]
    @test size(output) == (2, 7)

    rm(parsed_args["output"])
end

end;

true


text:
module TestExamples2DAdvection

using Test
using Trixi

include("test_trixi.jl")

# pathof(Trixi) returns /path/to/Trixi/src/Trixi.jl, dirname gives the parent directory
EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "2d")

@testset "Linear scalar advection" begin
  @testset "elixir_advection_basic.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_basic.jl"),
      l2   = [9.14468177884088e-6],
      linf = [6.437440532947036e-5])
  end

  @testset "elixir_advection_extended.jl with polydeg=1" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_extended.jl"),
      l2   = [0.052738240907090894],
      linf = [0.08754218386076529],
      polydeg=1)
  end

  @testset "elixir_advection_restart.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_restart.jl"),
      l2   = [1.2148032454832534e-5],
      linf = [6.495644795245781e-5])
  end

  @testset "elixir_advection_mortar.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_mortar.jl"),
      l2   = [0.0016133676309497648],
      linf = [0.013960195840607481])
  end

  @testset "elixir_advection_amr.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_amr.jl"),
      l2   = [8.553422725808076e-5],
      linf = [0.0007746920671942853])
  end

  @testset "elixir_advection_amr_nonperiodic.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_amr_nonperiodic.jl"),
      l2   = [3.5536338520440995e-5],
      linf = [0.0005514506116067221])
  end

  @testset "elixir_advection_amr_solution_independent.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_amr_solution_independent.jl"),
      l2   = [7.944846759690297e-5],
      linf = [0.0007454287896710293])
  end

  @testset "elixir_advection_amr_visualization.jl" begin
    # To make CI tests work, disable showing a plot window with the GR backend of the Plots package
    # Xref: https://github.com/jheinen/GR.jl/issues/278
    # Xref: https://github.com/JuliaPlots/Plots.jl/blob/8cc6d9d48755ba452a2835f9b89d3880e9945377/test/runtests.jl#L103
    if !isinteractive()
      restore = get(ENV, "GKSwstype", nothing)
      ENV["GKSwstype"] = "100"
    end

    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_amr_visualization.jl"),
      l2   = [0.0010300631535183275],
      linf = [0.009109608720471729])

    # Restore GKSwstype to previous value (if it was set)
    if !isinteractive()
      if isnothing(restore)
        delete!(ENV, "GKSwstype")
      else
        ENV["GKSwstype"] = restore
      end
    end
  end

  @testset "elixir_advection_timeintegration.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_timeintegration.jl"),
      l2   = [2.8777424975442084e-5],
      linf = [0.00047677279220256774])
  end

  @testset "elixir_advection_timeintegration.jl with carpenter_kennedy_erk43" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_timeintegration.jl"),
      l2   = [2.9389967058539404e-5],
      linf = [0.0005734205598502823],
      ode_algorithm=Trixi.CarpenterKennedy2N43(),
      cfl = 1.0)
  end

  @testset "elixir_advection_timeintegration.jl with carpenter_kennedy_erk43 with maxiters=1" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_timeintegration.jl"),
      l2   = [1.2046995151616365e-5],
      linf = [9.814863096402338e-5],
      ode_algorithm=Trixi.CarpenterKennedy2N43(),
      cfl = 1.0,
      maxiters = 1)
  end

  @testset "elixir_advection_timeintegration.jl with parsani_ketcheson_deconinck_erk94" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_timeintegration.jl"),
      l2   = [2.8778180933565268e-5],
      linf = [0.0004768629208851405],
      ode_algorithm=Trixi.ParsaniKetchesonDeconinck3Sstar94())
  end

  @testset "elixir_advection_timeintegration.jl with parsani_ketcheson_deconinck_erk32" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_timeintegration.jl"),
      l2   = [4.9135265795604245e-5],
      linf = [0.0005480422381496791],
      ode_algorithm=Trixi.ParsaniKetchesonDeconinck3Sstar32(),
      cfl = 1.0)
  end

  @testset "elixir_advection_timeintegration.jl with parsani_ketcheson_deconinck_erk32 with maxiters=1" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_timeintegration.jl"),
      l2   = [1.2074110061802306e-5],
      linf = [9.826959330024553e-5],
      ode_algorithm=Trixi.ParsaniKetchesonDeconinck3Sstar32(),
      cfl = 1.0,
      maxiters = 1)
  end

  @testset "elixir_advection_callbacks.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_callbacks.jl"),
      l2   = [9.14468177884088e-6],
      linf = [6.437440532947036e-5])
  end
end

end # module


text:
module H5WriteBuffer
using HDF5
import Base.push!,
       Base.close,
       HDF5.flush

export FileBackedBuffer,
       push!,
       dset,
       flush,
       close

### This type is designed to streamline appending, in the
### single most-major dimension, to an array backed by an
### HDF5 file. It is not designed to work for appending to
### an array of dimension 1, because that's not useful for
### its original purpose, but it wouldn't be difficult to
### extend that.
type FileBackedBuffer{T,N}
    handle               # Keep HDF5 file handle.
    take_ownership       # Take ownership of file handle.
    A                    # Reference to array inside.
    selectall            # A tuple of colons to select inner dims.
    bsize::Int           # Present size of in-memory buffer.
    bcapacity::Int       # Present capacity of in-memory buffer.
    b::Array{T,N}        # Buffer.
end
function FileBackedBuffer(hdf5_handle, dset_name::String,
                          T::Type, inner_size::Tuple, buffer_size::Int;
                          take_ownership=false)
    if exists(hdf5_handle, dset_name)
        A = d_open(hdf5_handle, dset_name)
    else
        initial_size = (inner_size..., buffer_size)
        max_size = (inner_size..., -1)
        A = d_create(hdf5_handle, dset_name, T,
                     (initial_size, max_size),
                     "chunk", (inner_size...,buffer_size))
        set_dims!(A, (inner_size...,0))
    end
    FileBackedBuffer(hdf5_handle, take_ownership, A,
                     (Colon() for i=1:length(size(A))-1),
                     0, buffer_size,
                     Array(T, (inner_size...,buffer_size)))
end
function FileBackedBuffer(fname::String, dset_name::String,
                          T::Type, inner_size::Tuple, buffer_size::Int)
    handle = h5open(fname, isfile(fname) && (filesize(fname)>0) ? "r+" : "w")
    FileBackedBuffer(handle, dset_name::String,
                     T::Type, inner_size::Tuple, buffer_size::Int;
                     take_ownership=true)
end
function close{T}(this::FileBackedBuffer{T})
    if this.take_ownership
        close(this.handle)
    end
end
function push!{T}(this::FileBackedBuffer{T}, x)
    this.bsize += 1
    this.b[this.selectall...,this.bsize] = x

    if this.bsize == this.bcapacity
        flush(this)
    end
end

# Return a reference to the dataset backing this buffer.
function dset{T}(this::FileBackedBuffer{T})
    flush(this)
    return this.A
end

# Write in-memory buffer to file and clear it.
function flush{T}(this::FileBackedBuffer{T})
    if this.bsize > 0
        set_dims!(this.A, (size(this.A)[1:end-1]..., size(this.A)[end]+this.bsize))
        this.A[this.selectall...,end-this.bsize+1:end] = this.b[this.selectall...,1:this.bsize]
    end
    flush(this.handle)
    _clear!(this)
end

# Clear in-memory buffer.
function _clear!{T}(this::FileBackedBuffer{T})
    this.bsize = 0
end

end # module


text:
# Register numbering, bit mask is 1 << x for x the enum value
perf_regs_numbering = Dict{UInt32, Symbol}(
   0 => :rax,  1 => :rbx,  2 => :rcx,  3 => :rdx,  4 => :rsi,
   5 => :rdi,  6 => :rbp,  7 => :rsp,  8 => :rip,  9 => :rflags,
  10 =>  :cs, 11 =>  :ss, 12 =>  :ds, 13 =>  :es, 14 =>  :fs,
  15 =>  :gs, 16 =>  :r8, 17 =>  :r9, 18 => :r10, 19 => :r11,
  20 => :r12, 21 => :r13, 22 => :r14, 23 => :r15
)


text:
# plots kernel density summary of data, with 90% acceptance region, mean and median.
using KernelDensity, Plots

function npdensity(z)
    #pyplot()
    n = size(z,2)
    p = zeros(n)
    for i = 1:n
        x = z[:,i]
        y = kde(x)
        q05 = quantile(x,0.05)
        Plots.plot(range(minimum(x), stop=q05, length=100),z->pdf(y,z), color=:orange, fill=(0,0.5,:orange), label="q05")
        q95 = quantile(x,0.95)
        Plots.plot!(range(q05, stop=q95, length=100),z->pdf(y,z), color=:green, fill=(0,0.5,:green),label="90% non-rejection region")
        Plots.plot!(range(q95, stop=maximum(x), length=100),z->pdf(y,z), color = :orange, fill=(0,0.5,:orange),label="q95")
        m = mean(x)
        Plots.plot!([m,m],[0,pdf(y,m)],color=:blue, label="mean")
        m = median(x)
        if i == 1
            p = Plots.plot!([m,m],[0,pdf(y,m)],color=:black, label="median")
        else
            p = [p, Plots.plot!([m,m],[0,pdf(y,m)],color=:black, label="median")]
        end
    end
    return p
end

max_stars_repo_name: nilshg/Econometrics

max_issues_repo_name: nilshg/Econometrics

max_forks_repo_name: nilshg/Econometrics


text:
import REPL.Terminals: CSI
const inner_width = 45
const inner_liner = "─"^inner_width

function Base.show(io::IO, frame::PredictionFrame)
    print_frame_table((args...)->nothing, frame)
end

function Base.show(io::IO, frames::Vector{PredictionFrame})
    print_frame_table((args...)->nothing, frames)
end

print_frame_table(image_callback, frame::PredictionFrame) =
    print_frame_table(image_callback, [frame])
function print_frame_table(image_callback, frames::Vector{PredictionFrame})
    lines, cols = displaysize(stdin)
    per_row = div(cols-1, 47)
    remaining = length(frames)
    offset = 0
    buf = IOBuffer()
    io = IOContext(buf, :color => true)
    first_row = min(per_row, remaining)
    println(io, "┌",join(("─"^inner_width for i = 1:first_row),"┬"),"┐")
    while remaining > 0
        this_row = min(per_row, remaining)
        row_frames = frames[(1:this_row) .+ offset]
        if any(p->p.filename !== nothing, row_frames)
            print(io, "│")
            print(io, join(map(row_frames) do p
                sprint() do io
                    fname = p.filename === nothing ? "" : p.filename
                    if length(fname) > inner_width
                        fname = basename(fname)
                    end
                    if length(fname) > inner_width
                        fname = string("…", fname[(end-(inner_width-1)):end])
                    end
                    padding = inner_width - length(fname)
                    print(io, " "^floor(Int, padding/2), fname, " "^ceil(Int, padding/2))
                end
            end, "│"))
            println(io, "│")
        end
        image_callback(io, row_frames)
        println(io, "├",join(("─"^inner_width for i = 1:this_row),"┼"),"┤")
        max_predictions = maximum(map(p->length(p.prediction.sorted_predictions), row_frames))
        any_correct = Symbol[:no for i = 1:this_row]
        for row in 1:max_predictions
            print(io, "│")
            for (col, p) in enumerate(row_frames)
                this_correct = false
                pred = p.prediction.sorted_predictions[row]
                if p.ground_truth === nothing
                    any_correct[col] = :yes
                else
                    ground_truth = p.ground_truth
                    if typeof(pred[1]) != typeof(ground_truth)
                        if hasmethod(convert, Tuple{Type{typeof(pred[1])}, typeof(ground_truth)})
                            ground_truth = convert(typeof(pred[1]), ground_truth)
                        else
                            ground_truth = nothing
                        end
                    end
                    if ground_truth !== nothing
                        this_correct = pred[1] == ground_truth
                        this_correct && (any_correct[col] = :yes)
                    else
                        any_correct[col] = :mismatch
                    end
                end
                print_prediction_row(io, pred,
                    annotation = this_correct ? :correct : :none)
                print(io," │")
            end
            println(io)
        end
        if !all(x->x===:yes, any_correct)
            print(io, "│")
            for (col, (had, p)) in enumerate(zip(any_correct, row_frames))
                if had === :yes
                    print(io," "^inner_width, "│")
                else
                    print_prediction_row(io, p.ground_truth=>p.ground_truth_confidence,
                        annotation = had == :mismatch ? :mismatch : :incorrect)
                    print(io," │")
                end
            end
            println(io)
        end
        remaining -= this_row
        offset += this_row
        if remaining != 0
            seps = ['├']
            for i = 1:per_row
                push!(seps, i == per_row ?
                    remaining >= i ? '┤' : '┘' :
                    remaining >= i ? '┼' : '┴')
            end
            println(io, join(seps, inner_liner))
        else
            println(io, "└",join((inner_liner for i = 1:this_row),"┴"),"┘")
        end
        # Write out at the end of each row
        write(stdout, take!(buf))
    end
end

const block_eights = [ Char(0x2590-i) for i = 1:8 ]
function confidence_bar(width, confidence)
    nfull_blocks = trunc(Int, confidence*width)
    neights = trunc(Int, ((confidence*width) - nfull_blocks)*8)
    sprint() do io
        for i = 1:nfull_blocks
            print(io, block_eights[8])
        end
        if neights != 0
            print(io, block_eights[neights])
            nfull_blocks += 1
        end
        print(io, " "^(max(0, width - nfull_blocks)))
    end
end

function print_prediction_row(io::IO, r::Pair; annotation = :none)
    class, confidence = r
    label_text = sprint(print, class)
    label_text = split(label_text, ',')[1]
    # TODO: Maybe use charwidth at some point if we ever have unicode labels
    if length(label_text) > 20
        label_text = string(label_text[1:19],"…")
    end
    color = annotation == :none ? :default :
            annotation == :correct ? :green :
            annotation == :incorrect ? :red :
            annotation == :mismatch ? :yellow :
            error("Unknown annotation")
    printstyled(io, rpad(label_text, 22); color = color)
    print(io, '│', " ")
    if annotation !== :mismatch
        r = round(confidence*100, digits=1)
        printstyled(io, string(
            confidence_bar(15, confidence)),
            lpad(r ≈ 100 ? 100 : r, 4), '%';
            color = color)
    else
        print(io, " "^20)
    end
end

function Base.show(io::IO, p::Prediction)
    for pred in p.sorted_predictions
        print_prediction_row(io, pred)
        println(io)
    end
end