text:
using Test
using Expronicon

@testset "is_fn" begin
    @test is_fn(:(foo(x) = x))
    @test is_fn(:(x -> 2x))
end

@testset "uninferrable_typevars" begin
    def = @expr JLKwStruct struct Inferable1{T}
        x::Constaint{T, <(2)}
    end

    @test uninferrable_typevars(def) == []

    def = @expr JLKwStruct struct Inferable2{T}
        x::Constaint{Float64, <(2)}
    end

    @test uninferrable_typevars(def) == [:T]

    def = @expr JLKwStruct struct Inferable3{T, N}
        x::Int
        y::N
    end
    @test uninferrable_typevars(def) == [:T]

    def = @expr JLKwStruct struct Inferable4{T, N}
        x::T
        y::N
    end
    @test uninferrable_typevars(def) == []

    def = @expr JLKwStruct struct Inferable5{T, N}
        x::T
        y::Float64
    end

    @test uninferrable_typevars(def) == [:T, :N]
    @test uninferrable_typevars(def; leading_inferable=false) == [:N]
end

@testset "has_plain_constructor" begin
    def = @expr JLKwStruct struct Foo1{T, N}
        x::Int
        y::N
        Foo1{T, N}(x, y) where {T, N} = new{T, N}(x, y)
    end

    @test has_plain_constructor(def) == true

    def = @expr JLKwStruct struct Foo2{T, N}
        x::T
        y::N
        Foo2(x, y) = new{typeof(x), typeof(y)}(x, y)
    end

    @test has_plain_constructor(def) == false

    def = @expr JLKwStruct struct Foo3{T, N}
        x::Int
        y::N
        Foo3{T}(x, y) where T = new{T, typeof(y)}(x, y)
    end

    @test has_plain_constructor(def) == false

    def = @expr JLKwStruct struct Foo4{T, N}
        x::T
        y::N
        Foo4{T, N}(x::T, y::N) where {T, N} = new{T, N}(x, y)
    end

    @test has_plain_constructor(def) == false
end

@testset "is_kw_fn" begin
    @test is_kw_fn(:(
        function foo(x::Int; kw=1)
        end
    ))

    ex = :(function (x::Int; kw=1) end)
    @test is_kw_fn(ex)
    @test !is_kw_fn(true)

    @test !is_kw_fn(:(
        function foo(x::Int)
        end
    ))

    @test !is_kw_fn(:(
        function (x::Int)
        end
    ))
end

@testset "JLFunction(ex)" begin
    jlfn = JLFunction()
    @test jlfn.name === nothing

    @test_expr JLFunction function foo(x::Int, y::Type{T}) where {T <: Real}
        return x
    end

    def = @test_expr JLFunction function (x, y)
        return 2
    end
    @test is_kw_fn(def) == false

    def = @test_expr JLFunction function (x, y; kw=2)
        return "aaa"
    end
    @test is_kw_fn(def) == true

    @test_expr JLFunction (x, y)->sin(x)

    # canonicalize head when it's a block
    @test_expr JLFunction function (x::Int; kw=1) end

    ex = :(struct Foo end)
    @test_throws AnalysisError JLFunction(ex)
    ex = :(@foo(2, 3))
    @test_throws AnalysisError split_function_head(ex)

    ex = :(Foo(; a = 1) = new(a))
    @test JLFunction(ex).kwargs[1] == Expr(:kw, :a, 1)
end

@testset "JLStruct(ex)" begin
    @test JLField(;name=:x).name === :x
    @test JLField(;name=:x).type === Any
    @test JLStruct(;name=:Foo).name === :Foo

    ex = :(struct Foo
        x::Int
    end)

    jlstruct = JLStruct(ex)
    println(jlstruct)
    @test jlstruct.name === :Foo
    @test jlstruct.ismutable === false
    @test length(jlstruct.fields) == 1
    @test jlstruct.fields[1].name === :x
    @test jlstruct.fields[1].type === :Int
    @test jlstruct.fields[1].line isa LineNumberNode
    @test codegen_ast(jlstruct) == ex

    ex = :(mutable struct Foo{T, S <: Real} <: AbstractArray
        a::Float64

        function foo(x, y, z)
            new(1)
        end
    end)

    jlstruct = JLStruct(ex)
    println(jlstruct)
    @test jlstruct.ismutable == true
    @test jlstruct.name === :Foo
    @test jlstruct.typevars == Any[:T, :(S <: Real)]
    @test jlstruct.supertype == :AbstractArray
    @test jlstruct.misc[1] == ex.args[3].args[end]
    @test rm_lineinfo(codegen_ast(jlstruct)) == rm_lineinfo(ex)

    ex = quote
        """
        Foo
        """
        struct Foo
            "xyz"
            x::Int
            y

            Foo(x) = new(x)
            1 + 1
        end
    end
    ex = ex.args[2]
    jlstruct = JLStruct(ex)
    @test jlstruct.doc == "Foo\n"
    @test jlstruct.fields[1].doc == "xyz"
    @test jlstruct.fields[2].type === Any
    @test jlstruct.constructors[1].name === :Foo
    @test jlstruct.constructors[1].args[1] === :x
    @test jlstruct.misc[1] == :(1 + 1)
    ast = codegen_ast(jlstruct)
    @test ast.args[1] == GlobalRef(Core, Symbol("@doc"))
    @test ast.args[3] == "Foo\n"
    @test ast.args[4].head === :struct
    @test is_fn(ast.args[4].args[end].args[end-1])
    println(jlstruct)

    @test_throws AnalysisError split_struct_name(:(function Foo end))
end

@testset "JLKwStruct" begin
    @test JLKwField(;name=:x).name === :x
    @test JLKwField(;name=:x).type === Any
    @test JLKwStruct(;name=:Foo).name === :Foo

    def = @expr JLKwStruct struct ConvertOption
        include_defaults::Bool=false
        exclude_nothing::Bool=false
    end

    @test_expr codegen_ast_kwfn(def, :create) == quote
        function create(::Type{S}; include_defaults = false, exclude_nothing = false) where S <: ConvertOption
            ConvertOption(include_defaults, exclude_nothing)
        end
    end

    def = @expr JLKwStruct struct Foo1{N, T}
        x::T = 1
    end
    println(def)

    @test_expr codegen_ast_kwfn(def, :create) == quote
        function create(::Type{S}; x = 1) where {N, T, S <: Foo}
            Foo1{N, T}(x)
        end
        function create(::Type{S}; x = 1) where {N, S <: Foo{N}}
            Foo1{N}(x)
        end
    end

    @test_expr codegen_ast(def) == quote
        struct Foo1{N, T}
            x::T
        end
        function Foo1{N, T}(; x = 1) where {N, T}
            Foo1{N, T}(x)
        end
        function Foo1{N}(; x = 1) where N
            Foo1{N}(x)
        end
    end

    def = @expr JLKwStruct struct Foo2 <: AbstractFoo
        x = 1
        y::Int
    end

    @test_expr codegen_ast(def) == quote
        struct Foo2 <: AbstractFoo
            x
            y::Int
        end
        function Foo2(; x = 1, y)
            Foo2(x, y)
        end
    end

    ex = quote
        """
        Foo
        """
        mutable struct Foo
            "abc"
            a::Int = 1
            b

            Foo(x) = new(x)
            1 + 1
        end
    end
    ex = ex.args[2]
    jlstruct = JLKwStruct(ex)
    @test jlstruct.doc == "Foo\n"
    @test jlstruct.fields[1].doc == "abc"
    @test jlstruct.fields[2].name === :b
    @test jlstruct.constructors[1].name === :Foo
    @test jlstruct.misc[1] == :(1 + 1)
    println(jlstruct)

    def = @expr JLKwStruct struct Foo3
        a::Int = 1
        Foo3(;a = 1) = new(a)
    end

    @test_expr codegen_ast(def) == quote
        struct Foo3
            a::Int
            Foo3(; a = 1) = new(a)
        end
    end
end

@testset "codegen_match" begin
    ex = codegen_match(:x) do
        quote
            1 => true
            2 => false
            _ => nothing
        end
    end

    eval(codegen_ast(JLFunction(;name=:test_match, args=[:x], body=ex)))

    @test test_match(1) == true
    @test test_match(2) == false
    @test test_match(3) === nothing
end

@test sprint(print, AnalysisError("a", "b")) == "expect a expression, got b."

@testset "JLIfElse" begin
    jl = JLIfElse()
    jl.map[:(foo(x))] = :(x = 1 + 1)
    jl.map[:(goo(x))] = :(y = 1 + 2)
    jl.otherwise = :(error("abc"))
    println(jl)
    ex = codegen_ast(jl)
    @test ex.head === :if
    @test ex.args[1] == :(foo(x))
    @test ex.args[2].args[1] == :(x = 1 + 1)
    @test ex.args[3].head === :elseif
end

@testset "JLMatch" begin
    jl = JLMatch(:x)
    jl.map[1] = true
    jl.map[2] = :(sin(x))
    println(jl)
    ex = codegen_ast(jl)
    jl = JLFunction(;name=:test_match, args=[:x], body=ex)
    println(jl)
    eval(codegen_ast(jl))

    @test test_match(1) == true
    @test test_match(2) == sin(2)
    @test test_match(3) === nothing
end

@testset "JLFor" begin
    ex = :(for i in 1:10, j in 1:20,
            k in 1:10
        1 + 1
    end)
    jl = JLFor(ex)
    println(jl)
    @test codegen_ast(jl) == ex

    jl = JLFor(;vars=[:x], iterators=[:itr], kernel=:(x + 1))
    ex = codegen_ast(jl)
    @test ex.head === :for
    @test ex.args[1].args[1] == :(x = itr)
    @test ex.args[2] == :(x+1)

    ex = :(for i in 1:10
        1 + 1
    end)
    jl = JLFor(ex)
    println(jl)
    @test jl.vars == [:i]
    @test jl.iterators == [:(1:10)]
end


text:
print_all(Lexers.JuliaLexer, "julia")

tokentest(
    Lexers.JuliaLexer,
    "#= ... =#",
    COMMENT_MULTILINE => "#= ... =#",
)

tokentest(
    Lexers.JuliaLexer,
    """
    #= ...
       #= ... =#
    =#
    """,
    COMMENT_MULTILINE => "#= ...\n   #= ... =#\n=#",
    TEXT => "\n",
)

tokentest(
    Lexers.JuliaLexer,
    "# ...",
    COMMENT_SINGLE => "# ...",
)

tokentest(
    Lexers.JuliaLexer,
    "   # ...",
    TEXT => "   ",
    COMMENT_SINGLE => "# ...",
)

tokentest(
    Lexers.JuliaLexer,
    "[1]",
    PUNCTUATION => "[",
    NUMBER_INTEGER => "1",
    PUNCTUATION => "]",
)

tokentest(
    Lexers.JuliaLexer,
    "a[begin]",
    NAME => "a",
    PUNCTUATION => "[",
    KEYWORD => "begin",
    PUNCTUATION => "]",
)

tokentest(
    Lexers.JuliaLexer,
    "a[begin:end]",
    NAME => "a",
    PUNCTUATION => "[",
    KEYWORD => "begin",
    OPERATOR => ":",
    KEYWORD => "end",
    PUNCTUATION => "]",
)

tokentest(
    Lexers.JuliaLexer,
    "a[begin-1]",
    NAME => "a",
    PUNCTUATION => "[",
    KEYWORD => "begin",
    OPERATOR => "-",
    NUMBER_INTEGER => "1",
    PUNCTUATION => "]",
)

tokentest(
    Lexers.JuliaLexer,
    "(1,2.0)",
    PUNCTUATION => "(",
    NUMBER_INTEGER => "1",
    PUNCTUATION => ",",
    NUMBER_FLOAT => "2.0",
    PUNCTUATION => ")",
)

tokentest(
    Lexers.JuliaLexer,
    "1 .+ 2 // 3 × !1",
    NUMBER_INTEGER => "1",
    TEXT => " ",
    OPERATOR => ".+",
    TEXT => " ",
    NUMBER_INTEGER => "2",
    TEXT => " ",
    OPERATOR => "//",
    TEXT => " ",
    NUMBER_INTEGER => "3",
    TEXT => " ",
    OPERATOR => "×",
    TEXT => " ",
    OPERATOR => "!",
    NUMBER_INTEGER => "1",
)

tokentest(
    Lexers.JuliaLexer,
    ":symbol",
    STRING_CHAR => ":symbol",
)

tokentest(
    Lexers.JuliaLexer,
    " :symbol",
    TEXT => " ",
    STRING_CHAR => ":symbol",
)

tokentest(
    Lexers.JuliaLexer,
    "(:symbol,",
    PUNCTUATION => "(",
    STRING_CHAR => ":symbol",
    PUNCTUATION => ",",
)

tokentest(
    Lexers.JuliaLexer,
    ":type",
    STRING_CHAR => ":type",
)

tokentest(
    Lexers.JuliaLexer,
    "[:_1, :2]",
    PUNCTUATION => "[",
    STRING_CHAR => ":_1",
    PUNCTUATION => ",",
    TEXT => " ",
    OPERATOR => ":",
    NUMBER_INTEGER => "2",
    PUNCTUATION => "]",
)

tokentest(
    Lexers.JuliaLexer,
    "if x in y",
    KEYWORD => "if",
    TEXT => " ",
    NAME => "x",
    TEXT => " ",
    KEYWORD_PSEUDO => "in",
    TEXT => " ",
    NAME => "y",
)

tokentest(
    Lexers.JuliaLexer,
    "x in y ? false : true",
    NAME => "x",
    TEXT => " ",
    KEYWORD_PSEUDO => "in",
    TEXT => " ",
    NAME => "y",
    TEXT => " ",
    OPERATOR => "?",
    TEXT => " ",
    KEYWORD_CONSTANT => "false",
    TEXT => " ",
    OPERATOR => ":",
    TEXT => " ",
    KEYWORD_CONSTANT => "true",
)

tokentest(
    Lexers.JuliaLexer,
    "let x = 1; 2x end",
    KEYWORD => "let",
    TEXT => " ",
    NAME => "x",
    TEXT => " ",
    OPERATOR => "=",
    TEXT => " ",
    NUMBER_INTEGER => "1",
    PUNCTUATION => ";",
    TEXT => " ",
    NUMBER_INTEGER => "2",
    NAME => "x",
    TEXT => " ",
    KEYWORD => "end",
)

tokentest(
    Lexers.JuliaLexer,
    "local t; global s; const c",
    KEYWORD_DECLARATION => "local",
    TEXT => " ",
    NAME => "t",
    PUNCTUATION => ";",
    TEXT => " ",
    KEYWORD_DECLARATION => "global",
    TEXT => " ",
    NAME => "s",
    PUNCTUATION => ";",
    TEXT => " ",
    KEYWORD_DECLARATION => "const",
    TEXT => " ",
    NAME => "c",
)

tokentest(
    Lexers.JuliaLexer,
    "module M end",
    KEYWORD => "module",
    TEXT => " ",
    NAME => "M",
    TEXT => " ",
    KEYWORD => "end",
)

tokentest(
    Lexers.JuliaLexer,
    "baremodule M end",
    KEYWORD => "baremodule",
    TEXT => " ",
    NAME => "M",
    TEXT => " ",
    KEYWORD => "end",
)

tokentest(
    Lexers.JuliaLexer,
    "abstract type A <: T end",
    KEYWORD => "abstract type",
    TEXT => " ",
    NAME => "A",
    TEXT => " ",
    OPERATOR => "<:",
    TEXT => " ",
    NAME => "T",
    TEXT => " ",
    KEYWORD => "end"
)

tokentest(
    Lexers.JuliaLexer,
    "mutable struct Point{T}\n    x::T\nend",
    KEYWORD => "mutable struct",
    TEXT => " ",
    NAME_FUNCTION => "Point",
    PUNCTUATION => "{",
    NAME => "T",
    PUNCTUATION => "}",
    TEXT => "\n    ",
    NAME => "x",
    OPERATOR => "::",
    NAME => "T",
    TEXT => "\n",
    KEYWORD => "end",
)

tokentest(
    Lexers.JuliaLexer,
    "Base.@inline f(x) = 2x",
    NAME => "Base",
    OPERATOR => ".",
    NAME_DECORATOR => "@inline",
    TEXT => " ",
    NAME_FUNCTION => "f",
    PUNCTUATION => "(",
    NAME => "x",
    PUNCTUATION => ")",
    TEXT => " ",
    OPERATOR => "=",
    TEXT => " ",
    NUMBER_INTEGER => "2",
    NAME => "x",
)

tokentest(
    Lexers.JuliaLexer,
    "struct Empty{T} end",
    KEYWORD => "struct",
    TEXT => " ",
    NAME_FUNCTION => "Empty",
    PUNCTUATION => "{",
    NAME => "T",
    PUNCTUATION => "}",
    TEXT => " ",
    KEYWORD => "end",
)

tokentest(
    Lexers.JuliaLexer,
    "struct Empty{T} where T <: S end",
    KEYWORD => "struct",
    TEXT => " ",
    NAME_FUNCTION => "Empty",
    PUNCTUATION => "{",
    NAME => "T",
    PUNCTUATION => "}",
    TEXT => " ",
    KEYWORD_PSEUDO => "where",
    TEXT => " ",
    NAME => "T",
    TEXT => " ",
    OPERATOR => "<:",
    TEXT => " ",
    NAME => "S",
    TEXT => " ",
    KEYWORD => "end",
)

tokentest(
    Lexers.JuliaLexer,
    "primitive type Float16 <: AbstractFloat 16 end",
    KEYWORD => "primitive type",
    TEXT => " ",
    NAME => "Float16",
    TEXT => " ",
    OPERATOR => "<:",
    TEXT => " ",
    NAME => "AbstractFloat",
    TEXT => " ",
    NUMBER_INTEGER => "16",
    TEXT => " ",
    KEYWORD => "end",
)

tokentest(
    Lexers.JuliaLexer,
    """
    macro something(x...)
        # ...
    end
    """,
    KEYWORD => "macro",
    TEXT => " ",
    NAME_FUNCTION => "something",
    PUNCTUATION => "(",
    NAME => "x",
    OPERATOR => "...",
    PUNCTUATION => ")",
    TEXT => "\n    ",
    COMMENT_SINGLE => "# ...",
    TEXT => "\n",
    KEYWORD => "end",
    TEXT => "\n",
)

tokentest(
    Lexers.JuliaLexer,
    "' ', '\\n', '\\'', '\"', '\\u1234', '⻆'",
    STRING_CHAR => "' '",
    PUNCTUATION => ",",
    TEXT => " ",
    STRING_CHAR => "'\\n'",
    PUNCTUATION => ",",
    TEXT => " ",
    STRING_CHAR => "'\\''",
    PUNCTUATION => ",",
    TEXT => " ",
    STRING_CHAR => "'\"'",
    PUNCTUATION => ",",
    TEXT => " ",
    STRING_CHAR => "'\\u1234'",
    PUNCTUATION => ",",
    TEXT => " ",
    STRING_CHAR => "'⻆'",
)

tokentest(
    Lexers.JuliaLexer,
    """
    "\\\"\\\" ..."
    """,
    STRING => "\"",
    STRING_ESCAPE => "\\\"\\\"",
    STRING => " ...\"",
    TEXT => "\n",
)

tokentest(
    Lexers.JuliaLexer,
    """
    " \$x \$(let x = y + 1; x^2; end)..."
    """,
    STRING => "\" ",
    STRING_INTERPOL => "\$x",
    STRING => " ",
    STRING_INTERPOL => "\$",
    PUNCTUATION => "(",
    KEYWORD => "let",
    TEXT => " ",
    NAME => "x",
    TEXT => " ",
    OPERATOR => "=",
    TEXT => " ",
    NAME => "y",
    TEXT => " ",
    OPERATOR => "+",
    TEXT => " ",
    NUMBER_INTEGER => "1",
    PUNCTUATION => ";",
    TEXT => " ",
    NAME => "x",
    OPERATOR => "^",
    NUMBER_INTEGER => "2",
    PUNCTUATION => ";",
    TEXT => " ",
    KEYWORD => "end",
    PUNCTUATION => ")",
    STRING => "...\"",
    TEXT => "\n",
)

tokentest(
    Lexers.JuliaLexer,
    """r"[a-z]+\$xyz"m""",
    STRING_OTHER => "r\"[a-z]+\$xyz\"",
    NAME => "m",
)

tokentest(
    Lexers.JuliaLexer,
    "run(`foo \$(bar)`)",
    NAME_FUNCTION => "run",
    PUNCTUATION => "(",
    STRING_BACKTICK => "`foo ",
    STRING_INTERPOL => "\$",
    PUNCTUATION => "(",
    NAME => "bar",
    PUNCTUATION => ")",
    STRING_BACKTICK => "`",
    PUNCTUATION => ")",
)

tokentest(
    Lexers.JuliaLexer,
    """
    "\\\\"
    """,
    STRING => "\"",
    STRING_ESCAPE => "\\\\",
    STRING => "\"",
    TEXT => "\n",
)

tokentest(
    Lexers.JuliaLexer,
    """
    1_000_000_000 + 1.0e-9 * .121 / 1121.
    1f0 - 1E-12
    0b100_101_111
    0o12123535252
    0x4312afAfabC
    1234
    1_2_3_4
    """,
    NUMBER_INTEGER => "1_000_000_000",
    TEXT => " ",
    OPERATOR => "+",
    TEXT => " ",
    NUMBER_FLOAT => "1.0e-9",
    TEXT => " ",
    OPERATOR => "*",
    TEXT => " ",
    NUMBER_FLOAT => ".121",
    TEXT => " ",
    OPERATOR => "/",
    TEXT => " ",
    NUMBER_FLOAT => "1121.",
    TEXT => "\n",
    NUMBER_FLOAT => "1f0",
    TEXT => " ",
    OPERATOR => "-",
    TEXT => " ",
    NUMBER_FLOAT => "1E-12",
    TEXT => "\n",
    NUMBER_BIN => "0b100_101_111",
    TEXT => "\n",
    NUMBER_OCT => "0o12123535252",
    TEXT => "\n",
    NUMBER_HEX => "0x4312afAfabC",
    TEXT => "\n",
    NUMBER_INTEGER => "1234",
    TEXT => "\n",
    NUMBER_INTEGER => "1_2_3_4",
    TEXT => "\n",
)

tokentest(
    Lexers.JuliaLexer,
    "√",
    OPERATOR => "√",
)

tokentest(
    Lexers.JuliaLexer,
    "√α",
    OPERATOR => "√",
    NAME => "α",
)

ext: jl


text:
# include("../src/LarSurf.jl")
# using LinearAlgebraicRepresentation
# Lar = LinearAlgebraicRepresentation
using LarSurf
using Pio3d
using Plasm
# using LarSurf

threshold = 4000
pth = Pio3d.datasets_join_path("medical/orig/sample-data/nrn4.pklz")
datap = Pio3d.read3d(pth)

data3d = datap["data3d"]
# segmentation = convert(Array{Int8, 2}, data3d .> threshold)
segmentation = data3d .> threshold

segmentation_linear_int = [Int8(segmentation[k]) for k=1:length(segmentation)]

data_size = LarSurf.size_as_array(size(data3d))
b3, larmodel = LarSurf.get_boundary3(data_size , true)

V = larmodel[1]
VV, FV, CV = larmodel[2]

FVk = LarSurf.Lar.characteristicMatrix(FV)
segmentation_arr2 = reshape(segmentation_linear_int, length(segmentation), 1)

faces = b3' * segmentation_arr2

selected_faces = FVk .* faces
FVill = [selected_faces[k, :] for k=1:size(selected_faces,1) if any(selected_faces[k, :] .!= 0)]
# V, topology_unused, LarSurf.cuboidGrid(dat_size)
Plasm.view(V, FVill)(2, 2, 2)

# block_size = [5, 5, 5]
# margin_size = 0

## Artifical data
# data3d = ones(Int16, 3, 3, 3)
# data3d[2, 2, 2] = 10
# data3d

# voxelsize_mm = [0.5, 1.0, 2.]
# threshold=0

# include("../src/LarSurf.jl")

# import SparseArrays.spzeros
# using Plasm
# using LinearAlgebraicRepresentation
# Lar = LinearAlgebraicRepresentation

# block_size = [2, 2, 2]
# block_size = size(data3d)
# b3 = LarSurf.get_boundary3(block_size)

# blocks_number, blocks_number_axis = LarSurf.number_of_blocks_per_axis(
#     size(data3d), block_size)
#
#
# block1 = LarSurf.get_block(
#     data3d, block_size, margin_size, blocks_number_axis, false, 1
# )
#
# verts, trifaces = LarSurf.import_data3d(data3d, voxelsize_mm, 4000)
#
# println(verts)
#
#
# println(trifaces)
#
#
# println("======= To LAR ========")
#
# V, EV, FE = LarSurf.to_lar(verts, trifaces)
#
# println(V)
# println(EV)
# println(FE)
#
# println("V: ", typeof(V), size(V))
# println("EV: ", typeof(EV))
# print("FV: ", typeof(FE))
#
# @show EV;
#
#
# # convert to list of lists
# # trifaces_list = [trifaces[k, :] for k=1:size(trifaces,1)]
# trifaces_list = [trifaces[k, 3:-1:1] for k=1:size(trifaces,1)]
# Vt = permutedims(V, [2,1])
#
# # Plasm.view(Lar.cuboid([1,1,1]))
#
# # Plasm.view(Vt, trifaces_list)
# Plasm.viewexploded(V, selected_faces)(2,2,2)


text:
using Documenter, FDM

makedocs(
    modules=[FDM],
    format=:html,
    pages=[
        "Home" => "index.md",
        "API" => "pages/api.md"
    ],
    sitename="FDM.jl",
    authors="Invenia Labs",
    assets=[
        "assets/invenia.css",
    ],
)

deploydocs(
    repo = "github.com/invenia/FDM.jl.git",
    julia = "1.0",
    target = "build",
    deps = nothing,
    make = nothing,
)


text:
#=
Testing of frequencyfilt.jl

Hard to test these functions other than visually.  This script simply
tries to run then all to make sure that they at least run

Set disp = true to display images for visual verification
=#

disp = false

using ImagePhaseCongruency, Test, TestImages, AbstractFFTs

if disp
    using PyPlot
end

println("Testing frequencyfilt...")

function mypause()
    println("Hit return to continue ")
    a = readline()
    return nothing
end

if disp; PyPlot.set_cmap(PyPlot.ColorMap("gray")); end

lena = Float64.(testimage("lena_gray"))

rows = 101
cols = 200
(f, fx, fy) = filtergrids((rows, cols))
if disp; imshow(f); title("filtergrids() f"); mypause(); end

if disp; imshow(fx); title("filtergrids() fx"); mypause(); end
if disp; imshow(fy); title("filtergrids() fy"); mypause(); end

f = filtergrid((rows, cols))
if disp; imshow(f); title("filtergrid() f"); mypause(); end

(H1, H2, f) = monogenicfilters((rows, cols))
if disp; imshow(real(-im*H1)); title("H1"); mypause(); end
if disp; imshow(real(-im*H2)); title("H2"); mypause(); end
if disp; imshow(f); title("monogenicfilters() f"); mypause(); end

(H, f) = packedmonogenicfilters((rows, cols))
if disp; imshow(real(H)); title("packedmonogenicfilters() H1"); mypause(); end
if disp; imshow(imag(H)); title("packedmonogenicfilters() H2"); mypause(); end

sze = (rows,cols)
cutin = 0.1
cutoff = 0.2
n = 2
boost = 2

f = lowpassfilter(sze, cutoff, n)
if disp; imshow(f); title("lowpassfilter() cutoff 0.2"); mypause(); end

f = bandpassfilter(sze, cutin, cutoff, n)
if disp; imshow(f); title("bandpassfilter() 0.1 0.2"); mypause(); end

f = highboostfilter(sze, cutoff, n, boost)
if disp; imshow(f); title("highboostfilter() 0.2 "); mypause(); end

f = highpassfilter(sze, cutoff, n)
if disp; imshow(f); title("highpassfilter() 0.2"); mypause(); end

f = 0.2
fo = 0.3
sigmaOnf = 0.55
@test abs(loggabor(0, fo, sigmaOnf) - 0) < eps()
@test abs(loggabor(fo, fo, sigmaOnf) - 1) < eps()

(f, fx, fy) = filtergrids((rows, cols))
(sintheta, costheta) = gridangles(f, fx, fy)
if disp; imshow(sintheta); title("filtergrids() sintheta"); mypause(); end
if disp; imshow(costheta); title("filtergrids() costheta"); mypause(); end

angl = pi/4
wavelen = pi/8
flter = cosineangularfilter(angl, wavelen, sintheta, costheta)
if disp; imshow(flter); title("cosineangularfilter"); mypause(); end

thetaSigma = .4
flter = gaussianangularfilter(angl, thetaSigma, sintheta, costheta)
if disp; imshow(flter); title("gaussianangularfilter"); mypause(); end

(P, S, p, s) = perfft2(lena)
if disp; imshow(lena); title("lena"); mypause(); end
if disp; imshow(p); title("periodic lena"); mypause(); end
if disp; imshow(s); title("lena - periodic lena"); mypause(); end

s1 = geoseries(0.5, 2, 4)
s2 = geoseries((0.5, 4), 4)
@test sum(abs.(s1 .- [0.5000,    1.0000,    2.0000,    4.0000])) < eps()
@test sum(abs.(s2 .- [0.5000,    1.0000,    2.0000,    4.0000])) < eps()

nothing


text:
using Genie.Router, Genie.Renderer.Html, Genie.Renderer.Json

route("/") do
  html("Welcome to Genie! This site is served by Julia $(Base.VERSION) with $(Threads.nthreads()) threads.")
end

route("/random") do
  json(rand(Int, parse(Int64, @params(:num, "1"))))
end

max_stars_repo_name: essenciary/performance-tests

max_issues_repo_name: essenciary/performance-tests

max_forks_repo_name: essenciary/performance-tests


text:
"""
    p_func()

Lorem ipsum dolor sit amet.
"""
function p_func()
  "X"
end

"""
    p_func(cur_obj::GType)

Lorem ipsum dolor sit amet.
"""
function p_func(cur_obj::GType)
  "Y"
end

"""
    p_func(cur_obj::SType)

Lorem ipsum dolor sit amet.
"""
function p_func(cur_obj::SType)
  "Z"
end


text:
export ft_seteventnotification
export FT_EVENT_RXCHAR, FT_EVENT_MODEM_STATUS
export FT_EVENT_LINE_STATUS

#Notification Events (see ft_seteventnotification)
const FT_EVENT_RXCHAR = 1
const FT_EVENT_MODEM_STATUS = 2
const FT_EVENT_LINE_STATUS = 4

function ft_seteventnotification(ft_handle::Culong,
                                 event_mask::Integer,
                                 event_handle::Culong)
  eh = Ref{Culong}(event_handle)
  ft_status = ccall((:FT_SetEventNotification, d2xx),
                   Cuint,
                   (Culong, Cuint,Ref{Culong}),
                   ft_handle, event_mask, eh)
  checkstatus(ft_status)
  return nothing
end

#=
const _k32 = "Kernel32.dll"
hevent = ccall((:CreateEventW,_k32),
              Culonglong,
              (Cchar, Cchar, Cchar, Ptr{UInt8}),
              C_NULL,0x00,0x00,"")

=#


text:
# test to see if update partition works
# Claudia May 2015

if !isdefined(:individualtest) individualtest = false; end

if(individualtest)
    include("../src/types.jl")
    include("../src/functions.jl")
end

tree = "(((((((1,2),3),4),5),(6,7)),(8,9)),10);"

#seed = 2738
seed = 56326

currT0 = readTopologyLevel1(tree);
srand(seed);
besttree = deepcopy(currT0);
success,hybrid,flag,nocycle,flag2,flag3 = addHybridizationUpdate!(besttree);
success
#printEdges(besttree)
#printNodes(besttree)
writeTopologyLevel1(besttree,true)
net = deepcopy(besttree);
length(net.partition)
[n.number for n in net.partition[1].edges] == [15] || error("wrong partition")
[n.number for n in net.partition[2].edges] == [11] || error("wrong partition")
[n.number for n in net.partition[3].edges] == [10] || error("wrong partition")
[n.number for n in net.partition[4].edges] == [9,7,5,3,1,2,4,6,8] || error("wrong partition")
[n.number for n in net.partition[5].edges] == [17] || error("wrong partition")
[n.number for n in net.partition[6].edges] == [14] || error("wrong partition")

success = false
success,hybrid,flag,nocycle,flag2,flag3 = addHybridizationUpdate!(besttree);
success
#printEdges(besttree)
#printNodes(besttree)
@test_nowarn writeTopologyLevel1(besttree,true)
net = deepcopy(besttree);
length(net.partition)
[n.number for n in net.partition[1].edges] == [15] || error("wrong partition")
[n.number for n in net.partition[2].edges] == [11] || error("wrong partition")
[n.number for n in net.partition[3].edges] == [10] || error("wrong partition")
[n.number for n in net.partition[4].edges] == [17] || error("wrong partition")
[n.number for n in net.partition[5].edges] == [14] || error("wrong partition")
[n.number for n in net.partition[6].edges] == [3,1,2] || error("wrong partition")
[n.number for n in net.partition[7].edges] == [21] || error("wrong partition")
[n.number for n in net.partition[8].edges] == [8] || error("wrong partition")
[n.number for n in net.partition[9].edges] == [6] || error("wrong partition")
[n.number for n in net.partition[10].edges] == [4] || error("wrong partition")

deleteHybridizationUpdate!(net,net.node[21]);
length(net.partition)
length(net.partition) == 6 || error("wrong partition")
# 15,11,10,[9,7,5,3,1,2,4,6,8],17,14
[n.number for n in net.partition[1].edges] == [15] || error("wrong partition")
[n.number for n in net.partition[2].edges] == [11] || error("wrong partition")
[n.number for n in net.partition[3].edges] == [10] || error("wrong partition")
[n.number for n in net.partition[4].edges] == [17] || error("wrong partition")
[n.number for n in net.partition[5].edges] == [14] || error("wrong partition")
[n.number for n in net.partition[6].edges] == [3,1,2,21,8,6,4,9,7] || error("wrong partition")
#printNodes(net)

deleteHybridizationUpdate!(net,net.node[18]);
length(net.partition) == 0 || error("wrong partition")


text:
function logistic(r,x)
    r*x*(1.0-x)
end

x0 = 0.2
r = 3.7
results = [(x0,logistic(r,x0))];
# We've already seen arrays: [1,2,3,4]
# Here I'm using tuples: (1,2,3,4)
# Results is an array filled with tuples!
N = 100; # or any big number
for i in 1:(N-1)
  xold = results[end][2]
  xnew = logistic(r,xold)
  append!(results,[(xold,xnew)])
end

print(results)

function logisticmap(x0=0.2,r=4,n=10)
    results = [(x0,logistic(r,x0))];
    N = 100;
    for i in 1:(N-1)
      xold = results[end][2]
      xnew = logistic(r,xold)
      append!(results,[(xold,xnew)])
      print(results)
  end
  results
end

using PlotlyJS
using Distributions

labs=["1","2","3","4","5"]
N=100;
x=1:N;
y = for i in 1:(N-1)
  xold = results[end][2]
  xnew = logistic(r,xold)
  append!(results,[(xold,xnew)])
end

layout = layout(
    title = "Logistical Map"
    xaxis_title="Xn"
    yaxis_title="Xn+1"
    font_family = "Arial",
    font_size = 12
)

plot()

labs=["1"]
function logisticmap(x0,r=4,n=10)
    results = [(x0,logistic(r,x0))];
    N = 100;
    for i in 1:(N-1)
      xold = results[end][2]
      xnew = logistic(r,xold)
      append!(results,[(xold,xnew)])
      print(results)
  end
  results
end
x = xold
y = xnew

layout = layout(
    title = "Logistical Map"
    xaxis_title="Xn"
    yaxis_title="Xn+1"
    font_family ="Arial",
    font_size = 12
)

plot()

function logisticmaping()
    trace1 = scatter(;x=xold, y=xnew, mode="lines+markers")
    plot([trace1])
end

function logisticmaping()
    function logisticmap(x0=0.2,r=4,n=10)
        results = [(x0,logistic(r,x0))];
        N = 100;
        for i in 1:(N-1)
          xold = results[end][2]
          xnew = logistic(r,xold)
          append!(results,[(xold,xnew)])
          print(results)
      end
      results
    end
    trace1 = scatter(;x=xold, y=xnew, mode="lines+markers")
    plot([trace1])
end


text:
# Autogenerated wrapper script for CSFML_jll for x86_64-apple-darwin
export libcsfml_audio, libcsfml_graphics, libcsfml_network, libcsfml_system, libcsfml_window

using SFML_jll
JLLWrappers.@generate_wrapper_header("CSFML")
JLLWrappers.@declare_library_product(libcsfml_audio, "@rpath/libcsfml-audio.2.5.dylib")
JLLWrappers.@declare_library_product(libcsfml_graphics, "@rpath/libcsfml-graphics.2.5.dylib")
JLLWrappers.@declare_library_product(libcsfml_network, "@rpath/libcsfml-network.2.5.dylib")
JLLWrappers.@declare_library_product(libcsfml_system, "@rpath/libcsfml-system.2.5.dylib")
JLLWrappers.@declare_library_product(libcsfml_window, "@rpath/libcsfml-window.2.5.dylib")
function __init__()
    JLLWrappers.@generate_init_header(SFML_jll)
    JLLWrappers.@init_library_product(
        libcsfml_audio,
        "lib/libcsfml-audio.2.5.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libcsfml_graphics,
        "lib/libcsfml-graphics.2.5.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libcsfml_network,
        "lib/libcsfml-network.2.5.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libcsfml_system,
        "lib/libcsfml-system.2.5.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libcsfml_window,
        "lib/libcsfml-window.2.5.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This code is licensed under MIT license (see LICENSE.md for full details)
#############################################################################
# EAGO
# A development environment for robust and global optimization
# See https://github.com/PSORLab/EAGO.jl
#############################################################################
# Defines functions used to parse the input optimization problem into
# a solvable form including routines used to classify input problems as a
# LP, SOCP, MILP, MISOCP, and convex problem types.
#############################################################################

#include(joinpath(@__DIR__, "parse_bridge_nl", "parse_bridge_nl.jl"))
include(joinpath(@__DIR__, "label_branch.jl"))

const INTEGER_FUNC_SET{T <: Real} = Union{Tuple{MOI.SingleVariable, MOI.ZeroOne},
                                          Tuple{MOI.SingleVariable, MOI.Integer},
                                          Tuple{MOI.SingleVariable, MOI.Semicontinuous},
                                          Tuple{MOI.SingleVariable, MOI.Semiinteger},
                                          Tuple{MOI.VectorOfVariables, MOI.SOS1{T}},
                                          Tuple{MOI.VectorOfVariables, MOI.SOS2{T}}}

const LINEAR_FUNC_SET{T <: Real} = Union{Tuple{MOI.SingleVariable, MOI.LessThan{T}},
                                         Tuple{MOI.SingleVariable, MOI.GreaterThan{T}},
                                         Tuple{MOI.SingleVariable, MOI.EqualTo{T}},
                                         Tuple{MOI.SingleVariable, MOI.Interval{T}},
                                         Tuple{SAF{T}, MOI.LessThan{T}},
                                         Tuple{SAF{T}, MOI.GreaterThan{T}},
                                         Tuple{SAF{T}, MOI.EqualTo{T}},
                                         Tuple{SAF{T}, MOI.Interval{T}},
                                         Tuple{MOI.VectorOfVariables,    MOI.Nonnegatives},
                                         Tuple{MOI.VectorOfVariables,    MOI.Nonpositives},
                                         Tuple{MOI.VectorOfVariables,    MOI.Zeros},
                                         Tuple{MOI.VectorAffineFunction{T}, MOI.Nonnegatives},
                                         Tuple{MOI.VectorAffineFunction{T}, MOI.Nonpositives},
                                         Tuple{MOI.VectorAffineFunction{T}, MOI.Zeros}}

const CONE_FUNC_SET{T <: Real} =  Union{Tuple{MOI.VectorOfVariables, MOI.NormInfinityCone},
                                        Tuple{MOI.VectorOfVariables, MOI.NormOneCone},
                                        Tuple{MOI.VectorOfVariables, MOI.SecondOrderCone},
                                        Tuple{MOI.VectorOfVariables, MOI.RotatedSecondOrderCone},
                                        Tuple{MOI.VectorOfVariables, MOI.GeometricMeanCone},
                                        Tuple{MOI.VectorOfVariables, MOI.ExponentialCone},
                                        Tuple{MOI.VectorOfVariables, MOI.DualExponentialCone},
                                        Tuple{MOI.VectorOfVariables, MOI.RelativeEntropyCone},
                                        Tuple{MOI.VectorOfVariables, MOI.NormSpectralCone},
                                        Tuple{MOI.VectorOfVariables, MOI.NormNuclearCone},
                                        Tuple{MOI.VectorOfVariables, MOI.PowerCone{T}},
                                        Tuple{MOI.VectorOfVariables, MOI.DualPowerCone{T}},
                                        }

const SDP_FUNC_SET = Union{Tuple{MOI.VectorOfVariables, MOI.PositiveSemidefiniteConeTriangle},
                           Tuple{MOI.VectorOfVariables, MOI.PositiveSemidefiniteConeSquare},
                           Tuple{MOI.VectorOfVariables, MOI.RootDetConeTriangle},
                           Tuple{MOI.VectorOfVariables, MOI.RootDetConeSquare},
                           Tuple{MOI.VectorOfVariables, MOI.LogDetConeTriangle},
                           Tuple{MOI.VectorOfVariables, MOI.LogDetConeSquare}}

_in_prob(::Val{LP}, S, T)   = Tuple{S, T} <: LINEAR_FUNC_SET{Float64}
_in_prob(::Val{MILP}, S, T) = _in_prob(Val{LP}(), S, T) | (Tuple{S, T} <: INTEGER_FUNC_SET{Float64})
_in_prob(::Val{SOCP}, S, T) = _in_prob(Val{LP}(), S, T) | (Tuple{S, T} <: CONE_FUNC_SET)
_in_prob(::Val{SDP}, S, T)  = _in_prob(Val{SOCP}(), S, T) | (Tuple{S, T} <: SDP_FUNC_SET)

for F in (LP, MILP, SOCP, SDP)
    @eval function _parse_classify_problem(::Val{$F}, m::Optimizer)
        ip = _input_model(m)
        lc = MOI.get(ip, MOI.ListOfConstraints())
        lcf = filter(x -> !_in_prob(Val{$F}(), x[1], x[2]), lc)
        flag = isempty(lcf)
        @show flag
        @show lc
        @show lcf
        flag &= _input_nlp_data(m) === nothing
        @show flag
        if flag
            m._problem_type = $F
        end
        return flag
    end
end

#=
Labels branching and nonbranching variables, remakes the optimizer for correct
dimension of NodeBB.
=#
function _parse_classify_problem(::Val{MINCVX}, m::Optimizer{T}) where T<:AbstractFloat
    m._problem_type = MINCVX

    d = Dict{Int,Bool}()

    sqf_lt_list = MOI.get(m, MOI.ListOfConstraintIndices{SQF{T},LT{T}}())
    sqf_gt_list = MOI.get(m, MOI.ListOfConstraintIndices{SQF{T},GT{T}}())
    sqf_et_list = MOI.get(m, MOI.ListOfConstraintIndices{SQF{T},ET{T}}())
    foreach(ci -> _label_branch!(d, m, ci), sqf_lt_list)
    foreach(ci -> _label_branch!(d, m, ci), sqf_gt_list)
    foreach(ci -> _label_branch!(d, m, ci), sqf_et_list)

    soc_list = MOI.get(m, MOI.ListOfConstraintIndices{VECVAR,SOC_CONE}())
    exp_list = MOI.get(m, MOI.ListOfConstraintIndices{VECVAR,EXP_CONE}())
    pow_list = MOI.get(m, MOI.ListOfConstraintIndices{VECVAR,POW_CONE{T}}())
    psd_list = MOI.get(m, MOI.ListOfConstraintIndices{VECVAR,PSD_CONE}())
    foreach(ci -> _label_branch!(d, m, ci, start = 2), soc_list)
    foreach(ci -> _label_branch!(d, m, ci), exp_list)
    foreach(ci -> _label_branch!(d, m, ci), pow_list)
    foreach(ci -> _label_branch!(d, m, ci), psd_list)

    #Get branch variables from nonlinear expressions
    has_nl_obj = _label_nl(d, _input_nlp_data(m))

    branch_num = has_nl_obj ? 1 : 0
    branch_num += length(keys(d))

    m._solver = GlobalOptimizer{branch_num, T, typeof(_ext_type(m))}()
    m._solver._ext_type = _ext_type(m)

    return true
end

"""
Classifies the problem type
"""
function _parse_classify_problem!(m::Optimizer)
    #_parse_bridge_nl!(m)
    _parse_classify_problem(Val{LP}(), m)       && return
    _parse_classify_problem(Val{MILP}(), m)     && return
    _parse_classify_problem(Val{SOCP}(), m)     && return
    _parse_classify_problem(Val{SDP}(), m)      && return
    _parse_classify_problem(Val{MINCVX}(), m)   && return
    return
end


text:
function predict(prior::Distribution, fitted_obj, prior_μs, ss::NormalSamples)
    predict.(prior, fitted_obj, prior_μs, ss)
end

function predict(prior::Normal, A, prior_μ, s::NormalSample)
        z =  response(s)
        noise_var = var(s)
        λ_shrink = A/(noise_var + A)
        λ_shrink*z + (1-λ_shrink)*prior_μ
end


text:
using StructuredOptimization
using AbstractOperators
using AbstractOperators.BlockArrays
using ProximalOperators
using LinearAlgebra, Random
using DSP, FFTW
using Test

Random.seed!(0)

@testset "StructuredOptimization" begin

@testset "Calculus" begin
  include("test_proxstuff.jl")
end

@testset "Syntax" begin
  include("test_variables.jl")
  include("test_expressions.jl")
  include("test_AbstractOp_binding.jl")
  include("test_terms.jl")
end

@testset "Problem construction" begin
  include("test_problem.jl")
  include("test_build_minimize.jl")
end

@testset "Integration tests" begin
  include("test_usage_small.jl")
  include("test_usage.jl")
end

end


text:
module FEFMM_Tests
    using Test
    using FEFMM
    using UnicodePlots
    using LinearAlgebra: norm

    include("testmediums.jl")

    function convergence_test(hlist, test_fun, test_name)
        @testset "$test_name" begin
            linorm = zeros(length(hlist))
            l2norm = zeros(length(hlist))
            for i = 1:length(hlist)
                (x0, k2, t_exact) = test_fun(hlist[i])
                t_pred = FEFMM.fefmm(k2, repeat([hlist[i]], ndims(k2)), x0)
                err = t_exact.-t_pred
                linorm[i] = norm(err, Inf)
                l2norm[i] = norm(err, 2)/sqrt(length(err))
            end

            plt = lineplot(log10.(hlist),
                        log10.(linorm),
                        title = test_name,
                        name = "L_inf norm",
                        xlabel = "log_10(h)",
                        ylabel = "log_10(Error)",
                        xlim = [floor(minimum(log10.(hlist))), ceil(maximum(log10.(hlist)))],
                        ylim = [floor(minimum(log10.([linorm; l2norm]))), ceil(maximum(log10.([linorm; l2norm])))],
                        color = :blue)

            lineplot!(plt, log10.(hlist), log10.(l2norm), name = "mean L2 norm", color = :red)
            println(plt)
            @test reverse(linorm) == sort(linorm) #test for decreasing error with decreasing h
            @test reverse(l2norm) == sort(l2norm)
        end
    end

    @testset "FEFMM Tests" begin
        @testset "Solve Piecewise Quadratic" begin
            @test isapprox(FEFMM.solve_piecewise_quadratic(1.0,1.0,1.0), 2.0)
            @test isapprox(FEFMM.solve_piecewise_quadratic(1.0,1.0,1.0,1.0,1.0), 1.0+1/sqrt(2))
            @test isapprox(FEFMM.solve_piecewise_quadratic(1.0,1.0,100.0,1.0,1.0), 2.0)
            @test isapprox(FEFMM.solve_piecewise_quadratic(2.0,1.0,1.0,100.0,1.0), 1.0)
            @test isapprox(FEFMM.solve_piecewise_quadratic(1.0,1.0,1.0,1.0,1.0,1.0,1.0), 1.0+1/sqrt(3))
            @test isapprox(FEFMM.solve_piecewise_quadratic(1.0,2.0,3.0,100.0,10.0,1.0,1.0), 2.0/3.0)
            @test isapprox(FEFMM.solve_piecewise_quadratic(1.0,2.0,3.0,10.0,100.0,1.0,1.0), 2.0/3.0)
            @test isapprox(FEFMM.solve_piecewise_quadratic(1.0,2.0,3.0,10.0,1.0,100.0,1.0), 1.0)
            @test isapprox(FEFMM.solve_piecewise_quadratic(1.0,2.0,3.0,1.0,10.0,100.0,1.0), 2.0)
        end

        @testset "Basic FEFMM Checks" begin
            κ21d = ones(101)
            dx1d = [0.1]
            xs1d = CartesianIndex(1)
            τa1d = 0:0.1:10
            @test all(isapprox.(τa1d, FEFMM.fefmm(κ21d, dx1d, xs1d)))
            κ22d = ones(101,101)
            dx2d = [0.1,0.1]
            xs2d = CartesianIndex(1,1)
            τaxis = 0:0.1:10
            τa2d = sqrt.(τaxis'.^2 .+ τaxis.^2)
            @test all(isapprox.(τa2d, FEFMM.fefmm(κ22d, dx2d, xs2d)))
            @test all(isapprox.(Float32.(τa2d), Float32.(FEFMM.fefmm(κ22d, Float32.(dx2d), xs2d))))
            κ23d = ones(101,101,101)
            dx3d = [0.1,0.1,0.1]
            xs3d = CartesianIndex(1,1,1)
            τa3d = zeros(101,101,101)
            for k in 1:101
                for j in 1:101
                    for i in 1:101
                        τa3d[i,j,k] = sqrt(τaxis[i]^2+τaxis[j]^2+τaxis[k]^2)
                    end
                end
            end
            @test all(isapprox.(τa3d, FEFMM.fefmm(κ23d, dx3d, xs3d)))
            @test all(isapprox.(Float32.(τa3d), Float32.(FEFMM.fefmm(κ23d, Float32.(dx3d), xs3d))))
        end

        @testset "Convergence Tests" begin
            convergence_test([1/20, 1/40, 1/80, 1/160, 1/320], const_k2_2D, "Constant Gradient of Squared Slowness 2D")
            convergence_test([1/20, 1/40, 1/80, 1/160, 1/320], const_v2_2D, "Constant Gradient of Squared Velocity 2D")
            convergence_test([1/20, 1/40, 1/80, 1/160, 1/320], gaussian_factor_2D, "Gaussian Factor 2D")
            convergence_test([1/20, 1/40, 1/80], const_k2_3D, "Constant Gradient of Squared Slowness 3D")
            convergence_test([1/20, 1/40, 1/80], const_v2_3D, "Constant Gradient of Squared Velocity 3D")
            convergence_test([1/20, 1/40, 1/80], gaussian_factor_3D, "Gaussian Factor 3D")
        end
    end
end


text:
function soscascade(signal::AbstractVector, coeffs::AbstractMatrix, state::AbstractMatrix)

    nSamples = size(signal, 1)
    nChannels = size(coeffs, 1)

    output = similar(signal, nSamples, nChannels)

    for n in 1:nSamples

        inp = signal[n]

        for i in 1:nChannels

            a0 = coeffs[i, 1]
            a1 = coeffs[i, 2]
            a2 = coeffs[i, 3]
            b1 = coeffs[i, 4]
            b2 = coeffs[i, 5]

            output[n, i] = a0 * inp + state[i, 1]

            state[i, 1] = a1 * inp - b1 * output[n, i] + state[i, 2]
            state[i, 2] = a2 * inp - b2 * output[n, i]

            inp = output[n, i]
        end

    end

    output, state

end


text:
function test_perf(G::Group)
    s = 0
    for g in G
        s += 1^g
    end
    return s
end

@testset "test_perf/benchmark iteration" begin

    G = SymmetricGroup(8)
    K = PermGroup([perm"(1,5,6,2,4,8)", perm"(1,3,6)(2,5,7,4)(8)"])
    @test test_perf(G) == test_perf(K) == 181440
    @info "Native iteration over S8 group:"
    @btime test_perf($G)
    # 1.429 ms (80644 allocations: 7.38 MiB) → 181440
    @info "Iteration over K ≅ S8 PermGroup:"
    @btime test_perf($K)
    # 17.457 ms (161282 allocations: 11.69 MiB) → 181440
end

@testset "GAP Docs examples/benchmarks" begin

    cube4 = Perm.([
       [1,9,3,11,5,13,7,15,2,10,4,12,6,14,8,16],
       [1,2,3,4,9,10,11,12,5,6,7,8,13,14,15,16],
       [1,2,5,6,3,4,7,8,9,10,13,14,11,12,15,16],
       [16,8,14,6,12,4,10,2,15,7,13,5,11,3,9,1],
       [3,11,1,9,7,15,5,13,4,12,2,10,8,16,6,14]]
       )
    G = PermGroup(cube4)
    @info "Schreier-Sims for Rubik cube-4 group:"
    @btime schreier_sims($(gens(G)));
    # 93.937 μs (2100 allocations: 225.73 KiB)
    @test order(G) == 384

    rubik_gens = [
    perm"( 1, 3, 8, 6)( 2, 5, 7, 4)( 9,33,25,17)(10,34,26,18)(11,35,27,19)",
    perm"( 9,11,16,14)(10,13,15,12)( 1,17,41,40)( 4,20,44,37)( 6,22,46,35)",
    perm"(17,19,24,22)(18,21,23,20)( 6,25,43,16)( 7,28,42,13)( 8,30,41,11)",
    perm"(25,27,32,30)(26,29,31,28)( 3,38,43,19)( 5,36,45,21)( 8,33,48,24)",
    perm"(33,35,40,38)(34,37,39,36)( 3, 9,46,32)( 2,12,47,29)( 1,14,48,27)",
    perm"(41,43,48,46)(42,45,47,44)(14,22,30,38)(15,23,31,39)(16,24,32,40)"
    ]

    rubik = PermGroup(rubik_gens)
    @info "Schreier-Sims for Rubik cube-9 group:"
    @btime schreier_sims($(gens(rubik)));
    # 43.665 ms (611137 allocations: 129.55 MiB)
    @test order(rubik) == 43252003274489856000
    @test order(rubik) == order(Int128, rubik) # fits Int128

    @testset "SL(4,7)" begin
        a = perm"""
    (  2,  3,  5,  9, 17, 32)(  4,  7, 13, 25, 48, 89)(  6, 11, 21, 40, 76,137)
    (  8, 15, 28, 54, 99,170)( 10, 19, 36, 69,127,113)( 12, 23, 44, 83,148,195)
    ( 14, 27, 52, 96,166,250)( 16, 30, 58,106,179,266)( 18, 34, 65,120,200,290)
    ( 20, 38, 72,131, 39, 74)( 22, 42, 79,142,225,125)( 24, 46, 86,153,220,115)
    ( 26, 50, 92,160,242,232)( 29, 56, 37, 71,130,211)( 31, 60,110,172,260,336)
    ( 33, 63,116,196,287,358)( 35, 67,124,119, 43, 81)( 41, 78,140,223,308, 84)
    ( 45, 80,144,159,117,177)( 47, 53, 93,162,244,194)( 49, 91,158,241,321,264)
    ( 51, 94, 64,118, 68,123)( 55,101,152,121,201,265)( 57,104,178,217,304,171)
    ( 59,108,185,245,324,374)( 61,112,192,285,357,379)( 62,114,191,283,279,351)
    ( 66,122,202, 73, 88,141)( 70,129,209,297,126, 97)( 75,135,216,291,342,247)
    ( 77,139,221,306,334,154)( 82,138,173,262,330,269)( 85,151,235,317,284,205)
    ( 87, 98,168,254,197,289)( 90,157,239,320,259,253)( 95,155,207,294,335,318)
    (100,163,169,256,331,181)(102,174,136,218,204,226)(103,176,175,263,325,156)
    (105,180,268,219,150,228)(107,183,273,257,165,248)(109,187,278,189,261,338)
    (111,190,282,300,366,362)(128,208,296,364,212,145)(132,143,161,243,315,372)
    (133,214,301,288,314,199)(134,206,292,361,376,370)(146,229,313,237,303,369)
    (147,231,164,246,193,286)(149,234,316,373,377,356)(167,252,327,328,305,323)
    (182,271,344,332,238,319)(184,274,186,276,341,382)(198,203,210,222,240,230)
    (213,299,295,215,267,298)(224,309,360,394,375,258)(227,270,343,293,340,307)
    (233,310,255,329,359,339)(236,251,322,311,371,352)(249,326,345,365,386,333)
    (272,346,383,378,275,348)(280,353,389,355,392,350)(302,367,312,337,380,363)
    (347,385,349,387,396,384)(354,391,398)(388,393,390)"""

        b = perm"""(  1,  2,  4,  8, 16, 31, 61,113,193,118,198,183,274, 69,128, 94,164,247,317,309,370,319,367,200,291, 67,125,206,293,362, 63,117, 74,134,215,302,368,392,391,  3,  6, 12, 24, 47, 88,155,238,276,287,359, 99,171,259,335,260,337,381,353,390,  9, 18, 35, 68,126,207,295,363, 40, 77, 54,100,172,261,127,104,179,267,341,283,356,153,236,244,323,218,160,214,262,209,168,255,330,377,361,374,346,384,380, 48, 90, 83,149, 96, 91,159,229,308,216,303,180,269,342,130, 50, 93,163,245,192,116,197, 58,107,184, 76,138,220,305,289,144,157,240, 60,111,191,284, 23, 45, 85,152,176,265, 92,161,110,189,281,355,393, 17, 33, 64,119,199,254,313,364,301,252,328,250,246,325,375,372,344,366,290,360,148,233,315,273,187,279,352, 72,132,213,300, 36, 70, 86,154,237,318,324,326, 65,121, 81,146,230,108,186,277,350,388,397,399,400)
        (  5, 10, 20, 39, 75,136,162,135,217,223,158,129,210,298,282, 11, 22, 43, 82,147,232,292,336,379, 21, 41, 28, 55,102,175,264,286,202,101,173,122,203,271,345,358,327,124,205, 52, 97,167,253,263,339,142,226,311,371,394,211,234,243,185,275,349,190, 89,156,170,258,334,294,343,383,396,357,  7, 14, 15, 29, 57,105,181,270,338,120,174,106,182,272,347,386, 13, 26, 51, 95,165,249, 19, 37, 46, 87,151,139,222,307,348,385,278, 25, 49, 44, 84,150,225,310,321,329,376,248,285,351,369,131,212,178,166,251,201, 42, 80,145,228,297,316,373,320,306,241,221,235,239,242,322,296,140, 78,141,224,304,331,299,365,114,194, 30, 59,109,188,280,354, 32, 62,115,195,268,231,314, 79,143,227,312,196,288, 38, 73,133,208, 71, 56,103,177, 27, 53, 98,169,257,333,137,219,266,340,112, 34, 66,123,204,256,332,378,387,382,395,389,398)"""

        SL_4_7 = PermGroup([a,b])
        @test order(SL_4_7) == 2317591180800
        @info "Schreier-Sims for SL(4,7):"
        @btime schreier_sims($(gens(SL_4_7)));
        # gap> G := Group([a,b]);; StabChain(G);; time;
        # GAP: ~15 ms vs
        # 220.512 ms (383827 allocations: 444.93 MiB)
    end

    @testset "DirectProduct example" begin
        a = perm"""
        (  1,401,400,399,398,397,396,395,394,393,392,391,390,389,388,387,386,
         385,384,383,382,381,380,379,378,377,376,375,374,373,372,371,370,369,
         368,367,366,365,364,363,362,361,360,359,358,357,356,355,354,353,352,
         351,350,349,348,347,346,345,344,343,342,341,340,339,338,337,336,335,
         334,333,332,331,330,329,328,327,326,325,324,323,322,321,320,319,318,
         317,316,315,314,313,312,311,310,309,308,307,306,305,304,303,302,301,
         300,299,298,297,296,295,294,293,292,291,290,289,288,287,286,285,284,
         283,282,281,280,279,278,277,276,275,274,273,272,271,270,269,268,267,
         266,265,264,263,262,261,260,259,258,257,256,255,254,253,252,251,250,
         249,248,247,246,245,244,243,242,241,240,239,238,237,236,235,234,233,
         232,231,230,229,228,227,226,225,224,223,222,221,220,219,218,217,216,
         215,214,213,212,211,210,209,208,207,206,205,204,203,202,201,200,199,
         198,197,196,195,194,193,192,191,190,189,188,187,186,185,184,183,182,
         181,180,179,178,177,176,175,174,173,172,171,170,169,168,167,166,165,
         164,163,162,161,160,159,158,157,156,155,154,153,152,151,150,149,148,
         147,146,145,144,143,142,141,140,139,138,137,136,135,134,133,132,131,
         130,129,128,127,126,125,124,123,122,121,120,119,118,117,116,115,114,
         113,112,111,110,109,108,107,106,105,104,103,102,101,100, 99, 98, 97,
          96, 95, 94, 93, 92, 91, 90, 89, 88, 87, 86, 85, 84, 83, 82, 81, 80,
          79, 78, 77, 76, 75, 74, 73, 72, 71, 70, 69, 68, 67, 66, 65, 64, 63,
          62, 61, 60, 59, 58, 57, 56, 55, 54, 53, 52, 51, 50, 49, 48, 47, 46,
          45, 44, 43, 42, 41, 40, 39, 38, 37, 36, 35, 34, 33, 32, 31, 30, 29,
          28, 27, 26, 25, 24, 23, 22, 21, 20, 19, 18, 17, 16, 15, 14, 13, 12,
          11, 10,  9,  8,  7,  6,  5,  4,  3,  2)"""
        b = perm"""
        (  1, 20,400,381)(  2, 40,399,361)(  3, 60,398,341)(  4, 80,397,321)
        (  5,100,396,301)(  6,120,395,281)(  7,140,394,261)(  8,160,393,241)
        (  9,180,392,221)( 10,200,391,201)( 11,220,390,181)( 12,240,389,161)
        ( 13,260,388,141)( 14,280,387,121)( 15,300,386,101)( 16,320,385, 81)
        ( 17,340,384, 61)( 18,360,383, 41)( 19,380,382, 21)( 22, 39,379,362)
        ( 23, 59,378,342)( 24, 79,377,322)( 25, 99,376,302)( 26,119,375,282)
        ( 27,139,374,262)( 28,159,373,242)( 29,179,372,222)( 30,199,371,202)
        ( 31,219,370,182)( 32,239,369,162)( 33,259,368,142)( 34,279,367,122)
        ( 35,299,366,102)( 36,319,365, 82)( 37,339,364, 62)( 38,359,363, 42)
        ( 43, 58,358,343)( 44, 78,357,323)( 45, 98,356,303)( 46,118,355,283)
        ( 47,138,354,263)( 48,158,353,243)( 49,178,352,223)( 50,198,351,203)
        ( 51,218,350,183)( 52,238,349,163)( 53,258,348,143)( 54,278,347,123)
        ( 55,298,346,103)( 56,318,345, 83)( 57,338,344, 63)( 64, 77,337,324)
        ( 65, 97,336,304)( 66,117,335,284)( 67,137,334,264)( 68,157,333,244)
        ( 69,177,332,224)( 70,197,331,204)( 71,217,330,184)( 72,237,329,164)
        ( 73,257,328,144)( 74,277,327,124)( 75,297,326,104)( 76,317,325, 84)
        ( 85, 96,316,305)( 86,116,315,285)( 87,136,314,265)( 88,156,313,245)
        ( 89,176,312,225)( 90,196,311,205)( 91,216,310,185)( 92,236,309,165)
        ( 93,256,308,145)( 94,276,307,125)( 95,296,306,105)(106,115,295,286)
        (107,135,294,266)(108,155,293,246)(109,175,292,226)(110,195,291,206)
        (111,215,290,186)(112,235,289,166)(113,255,288,146)(114,275,287,126)
        (127,134,274,267)(128,154,273,247)(129,174,272,227)(130,194,271,207)
        (131,214,270,187)(132,234,269,167)(133,254,268,147)(148,153,253,248)
        (149,173,252,228)(150,193,251,208)(151,213,250,188)(152,233,249,168)
        (169,172,232,229)(170,192,231,209)(171,212,230,189)(190,191,211,210)"""
        c = perm"""(402,403,404,405,406)"""
        d = perm"""(402,403)"""

        G = PermGroup([a,b,c,d])
        @test order(G) == 192480
        @info "Schreier-Sims for a direct-product group:"
        @btime schreier_sims($(gens(G)))
        # GAP: 22 ms vs
        # 34.596 ms (53786 allocations: 62.93 MiB)
        @time test_perf(G)
    end
end


text:
#   Copyright (c) 2020-2022 Velexi Corporation
#
#   Licensed under the Apache License, Version 2.0 (the "License");
#   you may not use this file except in compliance with the License.
#   You may obtain a copy of the License at
#
#       http://www.apache.org/licenses/LICENSE-2.0
#
#   Unless required by applicable law or agreed to in writing, software
#   distributed under the License is distributed on an "AS IS" BASIS,
#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
#   See the License for the specific language governing permissions and
#   limitations under the License.

"""
Unit tests for the +(x, y) function
"""

# --- Imports

# Standard library
import InteractiveUtils.subtypes
using Test

# GeometricAlgebra.jl
using GeometricAlgebra

# --- File inclusions

# Test utilities
include("test_utils.jl")

# --- Tests

# ------ M::Multivector

@testset "+(M::Multivector, N::Multivector)" begin
    @test_skip 1
end

@testset "+(M::Multivector, N::Blade)" begin
    @test_skip 1
end

@testset "+(M::Multivector, N::Pseudoscalar)" begin
    @test_skip 1
end

@testset "+(M::Multivector, N::Scalar)" begin
    @test_skip 1
end

@testset "+(M::Multivector, N::One)" begin
    @test_skip 1
end

@testset "+(M::Multivector, N::Zero)" begin
    @test_skip 1
end

@testset "+(M::Multivector, N::Real)" begin
    @test_skip 1
end

@testset "+(M::Multivector, N::Vector)" begin
    @test_skip 1
end

# ------ B::Blade

@testset "+(B::Blade, M::Multivector)" begin
    @test_skip 1
end

@testset "+(B::Blade, C::Blade)" begin
    @test_skip 1
end

@testset "+(B::Blade, C::Pseudoscalar)" begin
    @test_skip 1
end

@testset "+(B::Blade, C::Scalar)" begin
    @test_skip 1
end

@testset "+(B::Blade, C::One)" begin
    @test_skip 1
end

@testset "+(B::Blade, C::Zero)" begin
    test_dim = 10
    B = Blade(rand(test_dim, 3))
    C = Zero()
    @test B + C === B
end

@testset "+(B::Blade, C::Real)" begin
    @test_skip 1
end

@testset "+(B::Blade, C::Vector)" begin
    @test_skip 1
end

# ------ B::Pseudoscalar

@testset "+(B::Pseudoscalar, M::Multivector)" begin
    @test_skip 1
end

@testset "+(B::Pseudoscalar, C::Blade)" begin
    @test_skip 1
end

@testset "+(B::Pseudoscalar, C::Pseudoscalar)" begin
    # --- Preparations

    test_dim = 10

    test_value_1 = get_random_value(1)  # add 1 to keep value away from 0
    B = Pseudoscalar(test_dim, test_value_1)

    # --- Tests

    # B.value != -C.value
    test_value_2 = abs(test_value_1) + rand() + 1
    test_value_2 = rand() > 0.5 ? test_value_2 : -test_value_2
    C = Pseudoscalar(test_dim, test_value_2)

    B_plus_C = B + C
    @test B_plus_C isa Pseudoscalar
    @test B_plus_C.dim == test_dim
    @test B_plus_C.value == (test_value_1 + test_value_2)

    # B.value == -C.value
    C = Pseudoscalar(test_dim, -test_value_1)

    B_plus_C = B + C
    @test B_plus_C isa Zero
end

@testset "+(B::Pseudoscalar, C::Scalar)" begin
    @test_skip 1
end

@testset "+(B::Pseudoscalar, C::One)" begin
    @test_skip 1
end

@testset "+(B::Pseudoscalar, C::Zero)" begin
    test_dim = 10
    test_value = get_random_value(1)  # add 1 to keep value away from 0

    B = Pseudoscalar(test_dim, test_value)
    C = Zero()
    @test B + C === B
end

@testset "+(B::Pseudoscalar, C::Real)" begin
    @test_skip 1
end

@testset "+(B::Pseudoscalar, C::Vector)" begin
    @test_skip 1
end

# ------ B::Scalar

@testset "+(B::Scalar, M::Multivector)" begin
    @test_skip 1
end

@testset "+(B::Scalar, C::Blade)" begin
    @test_skip 1
end

@testset "+(B::Scalar, C::Pseudoscalar)" begin
    @test_skip 1
end

@testset "+(B::Scalar, C::Scalar)" begin
    # --- Preparations

    test_value_1 = get_random_value(2)  # add 2 to keep value away from 0 and 1
    B = Scalar(test_value_1)

    # --- Tests

    # B != -C, B != -C + 1
    test_value_2 = abs(test_value_1) + rand() + 2
    test_value_2 = rand() > 0.5 ? test_value_2 : -test_value_2
    C = Scalar(test_value_2)

    B_plus_C = B + C
    @test B_plus_C isa Scalar
    @test B_plus_C.value == test_value_1 + test_value_2

    # B == -C
    C = Scalar(-test_value_1)

    B_plus_C = B + C
    @test B_plus_C isa Zero

    # B == -C + 1
    C = Scalar(-test_value_1 + 1)

    B_plus_C = B + C
    @test B_plus_C isa One
end

@testset "+(B::Scalar, C::One)" begin
    # --- Preparations
    C = One()

    # --- Tests

    # B != -1
    test_value = get_random_value(2)  # add 2 to keep value away from 0 and 1
    B = Scalar(test_value)

    B_plus_C = B + C
    @test B_plus_C isa Scalar
    @test B_plus_C.value == 1 + test_value

    # B == -1
    B = Scalar(-1)

    B_plus_C = B + C
    @test B_plus_C isa Zero
end

@testset "+(B::Scalar, C::Zero)" begin
    test_value = get_random_value(2)  # add 2 to keep value away from 0 and 1
    B = Scalar(test_value)

    C = Zero()

    @test B + C === B
end

@testset "+(B::Scalar, C::Real)" begin
    # --- Preparations

    test_value_1 = get_random_value(2)  # add 2 to keep value away from 0 and 1
    B = Scalar(test_value_1)

    # --- Tests

    # C != -B, C != -B + 1
    test_value_2 = abs(test_value_1) + rand() + 2
    test_value_2 = rand() > 0.5 ? test_value_2 : -test_value_2
    C = test_value_2

    B_plus_C = B + C
    @test B_plus_C isa Scalar
    @test B_plus_C.value == test_value_1 + test_value_2

    # C == -B
    C = -test_value_1

    B_plus_C = B + C
    @test B_plus_C isa Zero

    # C == -B + 1
    C = -test_value_1 + 1

    B_plus_C = B + C
    @test B_plus_C isa One
end

@testset "+(B::Scalar, C::Vector)" begin
    @test_skip 1
end

# ------ B::One

@testset "+(B::One, M::Multivector)" begin
    @test_skip 1
end

@testset "+(B::One, C::Blade)" begin
    @test_skip 1
end

@testset "+(B::One, C::Pseudoscalar)" begin
    @test_skip 1
end

@testset "+(B::One, C::Scalar)" begin
    # --- Preparation

    B = One()

    # --- Tests

    # C != -1
    test_value = get_random_value(2)  # add 2 to keep value away from 0, 1, and -1
    C = Scalar(test_value)

    B_plus_C = B + C
    @test B_plus_C isa Scalar
    @test B_plus_C.value == 1 + test_value

    # C == -1
    C = Scalar(-1)

    B_plus_C = B + C
    @test B_plus_C isa Zero
end

@testset "+(B::One, C::One)" begin
    B = One()
    C = One()

    B_plus_C = B + C
    @test B_plus_C isa Scalar
    @test B_plus_C.value == 2
end

@testset "+(B::One, C::Zero)" begin
    B = One()
    C = Zero()
    @test B + C === B
end

@testset "+(B::One, C::Real)" begin
    # --- Preparation

    B = One()

    # --- Tests

    # C != -1, C != 0
    test_value = get_random_value(2)  # add 2 to keep value away from 0 and -1
    C = test_value

    B_plus_C = B + C
    @test B_plus_C isa Scalar
    @test B_plus_C.value == 1 + test_value

    # C == -1
    C = -1

    B_plus_C = B + C
    @test B_plus_C isa Zero

    # C == 0
    C = 0

    B_plus_C = B + C
    @test B_plus_C === B
end

@testset "+(B::One, C::Vector)" begin
    @test_skip 1
end

# ------ B::Zero

@testset "+(B::Zero, M::Multivector)" begin
    @test_skip 1
end

@testset "+(B::Zero, C::Blade)" begin
    test_dim = 10
    B = Zero()
    C = Blade(rand(test_dim, 3))
    @test B + C === C
end

@testset "+(B::Zero, C::Pseudoscalar)" begin
    test_dim = 10
    test_value = get_random_value(1)  # add 1 to keep value away from 0

    B = Zero()
    C = Pseudoscalar(test_dim, test_value)
    @test B + C === C
end

@testset "+(B::Zero, C::Scalar)" begin
    B = Zero()

    test_value = get_random_value(2)  # add 2 to keep value away from 0 and 1
    C = Scalar(test_value)

    @test B + C === C
end

@testset "+(B::Zero, C::One)" begin
    B = Zero()
    C = One()
    @test B + C === C
end

@testset "+(B::Zero, C::Zero)" begin
    B = Zero()
    C = Zero()
    @test B + C === B
end

@testset "+(B::Zero, C::Real)" begin
    # --- Preparation

    B = Zero()

    # --- Tests

    # C != 0, C != 1
    test_value = get_random_value(2)  # add 2 to keep value away from 0 and 1
    C = test_value

    B_plus_C = B + C
    @test B_plus_C isa Scalar
    @test B_plus_C.value == C

    # C == 0
    C = 0

    B_plus_C = B + C
    @test B_plus_C isa Zero

    # C == 1
    C = 1

    B_plus_C = B + C
    @test B_plus_C isa One
end

@testset "+(B::Zero, C::Vector)" begin
    B = Zero()
    C = rand(5)
    @test B + C === C
end

# ------ B::Real

@testset "+(B::Real, M::Multivector)" begin
    @test_skip 1
end

@testset "+(B::Real, C::Blade)" begin
    @test_skip 1
end

@testset "+(B::Real, C::Pseudoscalar)" begin
    @test_skip 1
end

@testset "+(B::Real, C::Scalar)" begin
    # --- Preparation

    test_value_1 = get_random_value(2)  # add 2 to keep value away from 0 and 1
    C = Scalar(test_value_1)

    # --- Tests

    # B != -C
    test_value_2 = abs(test_value_1) + rand() + 2
    test_value_2 = rand() > 0.5 ? test_value_2 : -test_value_2
    B = test_value_2

    B_plus_C = B + C
    @test B_plus_C isa Scalar
    @test B_plus_C.value == test_value_1 + test_value_2

    # B == -C
    B = -test_value_1

    B_plus_C = B + C
    @test B_plus_C isa Zero
end

@testset "+(B::Real, C::One)" begin
    # --- Preparation

    C = One()

    # --- Tests

    # B != 0, B != -1
    test_value = get_random_value(2)  # add 2 to keep value away from 0 and -1
    B = test_value

    B_plus_C = B + C
    @test B_plus_C isa Scalar
    @test B_plus_C.value == 1 + test_value

    # B == 0
    B = 0

    B_plus_C = B + C
    @test B_plus_C isa One

    # B == -1
    B = -1

    B_plus_C = B + C
    @test B_plus_C isa Zero
end

@testset "+(B::Real, C::Zero)" begin
    # --- Preparation

    C = Zero()

    # --- Tests

    # B != 0, B != 1
    test_value = get_random_value(2)  # add 2 to keep value away from 0 and 1
    B = test_value

    B_plus_C = B + C
    @test B_plus_C isa Scalar
    @test B_plus_C.value == B

    # B == 0
    B = 0

    B_plus_C = B + C
    @test B_plus_C isa Zero

    # B == 1
    B = 1

    B_plus_C = B + C
    @test B_plus_C isa One
end

# ------ B::Vector

@testset "+(B::Vector, M::Multivector)" begin
    @test_skip 1
end

@testset "+(B::Vector, C::Blade)" begin
    @test_skip 1
end

@testset "+(B::Vector, C::Pseudoscalar)" begin
    @test_skip 1
end

@testset "+(B::Vector, C::Scalar)" begin
    @test_skip 1
end

@testset "+(B::Vector, C::One)" begin
    @test_skip 1
end

@testset "+(B::Vector, C::Zero)" begin
    B = rand(5)
    C = Zero()
    @test B + C === B
end


text:
@inbounds begin
function _vcmp(x1::Vector{U}, x2::Vector{U}) where {U<:BitUnsigned}
    l1, l2 = length(x1), length(x2)
    if l2 == l1
        for k = l1:-1:1
            if x1[k] != x2[k]
                return cmp(x1[k], x2[k])
            end
        end
        return 0
    else
        return cmp(l1, l2)
    end
end

function _vleftshift!(x::Vector{U}, n::UInt) where {U<:BitUnsigned}
    n==0 && return x
    k = 1
    while n >= fullshift(U)
        n -= fullshift(U)
        unshift!(x, zero(U))
        k += 1
    end
    carry = zero(U)
    @inbounds while k <= length(x)
        # nextcarry = x[k] >> (fullshift(U)-n)
        # x[k] = (x[k] << n) + carry
        # carry = nextcarry
        carry, x[k] = (x[k] >> (fullshift(U)-n)) , ((x[k] << n) + carry)
        k += 1
    end
    !iszero(carry) && push!(x, carry)
    return x
end
function _vrightshift!(x::Vector{U}, n::UInt) where {U<:BitUnsigned}
    n==0 && return x
    k = 1
    while n >= fullshift(U) && length(x) > 0
        n -= fullshift(U)
        shift!(x)
    end
    length(x) == 0 && return x
    while k < length(x)
        x[k] = (x[k] >> n) | (x[k+1] << (fullshift(U)-n))
        k += 1
    end
    x[k] = (x[k] >> n)

    iszero(x[k]) && pop!(x)
    return x
end
function _vand!(y::Vector{U}, x::Vector{U}) where {U<:BitUnsigned}
    while length(y) > length(x)
        pop!(y)
    end
    @inbounds for k in 1:length(y)
        y[k] &= x[k]
    end
    while length(y) > 0 && iszero(last(y))
        pop!(y)
    end
    return y
end
function _vor!(y::Vector{U}, x::Vector{U}) where {U<:BitUnsigned}
    while length(y) < length(x)
        push!(y, zero(U))
    end
    @inbounds for k = 1:length(x)
        y[k] |= x[k]
    end
    return y
end

function _vadd!(y::Vector{U}, x::Vector{U}) where {U<:BitUnsigned}
    while length(y) < length(x)
        push!(y, zero(U))
    end
    lx, ly = length(x), length(y)
    carry = zero(U)
    for k = 1:lx
        y[k], ncarry = _add(y[k], x[k])
        if iszero(carry)
            carry = ncarry
        elseif iszero(ncarry)
            y[k], carry = _add(y[k], carry)
        else # if ncarry is nonzero, y[k] < typemax and ...
            y[k] += carry # ... this cannot overflow since carry = 1
            carry = ncarry
        end
    end
    for k = (lx+1):ly
        iszero(carry) && break
        y[k], carry = _add(y[k], carry)
    end
    !iszero(carry) && push!(y, carry)
    return y
end

function _vsub!(y::Vector{U}, x::Vector{U}) where {U<:BitUnsigned}
    _vcmp(y,x) < 0 && throw(OverflowError())
    lx, ly = length(x), length(y)
    borrow = zero(U)
    for k = 1:lx
        y[k], nborrow = _sub(y[k], x[k])
        if iszero(borrow)
            borrow = nborrow
        elseif iszero(nborrow)
            y[k], borrow = _sub(y[k], borrow)
        else # if nborrow is nonzero, y[k] > 0 and ...
            y[k] -= borrow # ... this cannot overflow since borrow = 1
            borrow = nborrow
        end
    end
    for k=(lx+1):ly
        iszero(borrow) && break
        y[k], borrow = _sub(y[k], borrow)
    end
    while length(y) > 0 && iszero(last(y)) # remove zero words at beginning
        pop!(y)
    end
    return y
end

function _vmuladd!(y::Vector{U}, x1::Vector{U}, x2::Vector{U}) where {U<:BitUnsigned}
    while length(y) < length(x1)+length(x2)
        push!(y, zero(U))
    end
    l1, l2 = length(x1), length(x2)
    for l = 1:l2
        # if !iszero(x2[l]) # not worthwhile according to Knuth
            carry = zero(U)
            for k = 1:l1
                prod, carry1 = _mul(x1[k], x2[l]) # carry1 < typemax - 1
                y[k+l-1], carry2 = _add(y[k+l-1], prod) # carry2 <= 1
                y[k+l-1], carry3 = _add(y[k+l-1], carry) # carry3 <= 1 - carry2
                carry = carry1 + carry2 + carry3 # this should not overflow
            end
            for k = (l1+l):(l1+l2)
                iszero(carry) && break
                y[k], carry = _add(y[k], carry)
            end
        # end
    end
    while iszero(last(y))
        pop!(y)
    end
    return y
end
function _vdivrem!(num::Vector{U}, den::U) where {U<:BitUnsigned} # n will be destroyed in the process
    iszero(den) && throw(DivideError())
    length(num) == 0 && return(num, num)
    length(num) == 1 && return divrem(num[1], den)

    last(num) >= den && push!(num, zero(U))
    q = Vector{U}(length(num)-1)
    k = length(q)
    while length(num) >= 2
        q[k], r = _divrem(num[end-1], num[end], den)
        pop!(num)
        num[end] = r
        k -= 1
    end
    iszero(last(num)) && pop!(num)
    return q, num
end

function _vdivrem!(num::Vector{U}, den::Vector{U}) where {U<:BitUnsigned} # n will be destroyed in the process
    isempty(den) && throw(diverr)
    c = _vcmp(num,den)
    c < 0 && return (Vector{U}(0), num)
    c == 0 && return ([one(U)], Vector{U}(0))
    length(den) == 1 && return _vdivrem!(num, den[1])

    # we now know that length(den) >= 2, length(num) >= length(den)
    #normalize
    s = unsigned(leading_zeros(last(den)))
    u = _vleftshift!(num, s)
    v = _vleftshift!(den, s) # d is modified but will be restored at the end

    n = length(v)
    # compare v with first n words of u
    m = length(u) - n
    k = n
    while k > 0 && u[m+k] == v[k]
        k -= 1
    end
    if k > 0 && u[m+k] < v[k] # smaller, but u in total larger
        m -= 1
    else
        push!(u, zero(U))
    end
    # now m = length(u) - n -1
    q = fill(zero(U), m+1)

    for j = m:-1:0
        q̂ = _estimateq(u[j+n-1],u[j+n],u[j+n+1], v[n-1], v[n])
        borrow = zero(U)
        for k=1:n # combination of _vmul! and _vadd!
            u[k+j], borrow1 = _sub(u[k+j], borrow) # borrow1 <= 1
            p, borrow2 = _mul(v[k], q̂) # borrow2 <= typemax - 2
            u[k+j], borrow3 = _sub(u[k+j], p) # borrow3 <= 1 - borrow1
            borrow = borrow1 + borrow2 + borrow3 # cannot overflow
        end
        u[j+n+1], borrow = _sub(u[n+j+1], borrow)
        if !iszero(borrow) # should happen rarely according to Knuth
            q̂ -= one(q̂)
            carry = zero(U)
            for k=1:n # add back: see _vadd!
                u[k+j], ncarry = _add(u[k+j], v[k])
                if iszero(carry)
                    carry = ncarry
                elseif iszero(ncarry)
                    u[k+j], carry = _add(u[k+j], carry)
                else
                    u[k+j] += carry
                    carry = ncarry
                end
            end
            u[n+j+1], carry = _add(u[n+j+1], carry)
            # @assert borrow == carry
        end
        # @assert u[n+j+1] == zero(u[n+j+1])
        q[j+1] = q̂
    end
    while length(u) > 0 && iszero(last(u))
        pop!(u)
    end
    # restore shift
    den = _vrightshift!(v, s) # d is modified but will be restored at the end
    r = _vrightshift!(u, s) # d is modified but will be restored at the end
    return q, r
end
end


text:
using Librsync
@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end

@test Librsync.MD4_SUM_LENGTH > 0
@test Librsync.BLAKE2_SUM_LENGTH > 0


text:
# benchmark suite execution and codespeed submission

using CUDA

using BenchmarkTools

using StableRNGs
rng = StableRNG(123)

# we only submit results when running on the master branch
real_run = get(ENV, "CODESPEED_BRANCH", nothing) == "master"
if real_run
    # to find untuned benchmarks
    BenchmarkTools.DEFAULT_PARAMETERS.evals = 0
end

# convenience macro to create a benchmark that requires synchronizing the GPU
macro async_benchmarkable(ex...)
    quote
        # use non-blocking sync to reduce overhead
        @benchmarkable CUDA.@sync $(ex...)
    end
end

# before anything else, run latency benchmarks. these spawn subprocesses, so we don't want
# to do so after regular benchmarks have caused the memory allocator to reserve memory.
@info "Running latency benchmarks"
latency_results = include("latency.jl")

SUITE = BenchmarkGroup()

# NOTE: don't use spaces in benchmark names (tobami/codespeed#256)

include("kernel.jl")
include("array.jl")

if real_run
    @info "Preparing main benchmarks"
    warmup(SUITE; verbose=false)
    tune!(SUITE)

    # reclaim memory that might have been used by the tuning process
    GC.gc(true)
    CUDA.reclaim()
end

# benchmark groups that aren't part of the suite
addgroup!(SUITE, "integration")

@info "Running main benchmarks"
results = run(SUITE, verbose=true)

# integration tests (that do nasty things, so need to be run last)
@info "Running integration benchmarks"
integration_results = BenchmarkGroup()
integration_results["volumerhs"] = include("volumerhs.jl")
integration_results["byval"] = include("byval.jl")
integration_results["cudadevrt"] = include("cudadevrt.jl")

results["latency"] = latency_results
results["integration"] = integration_results

println(results)

## comparison

# write out the results
BenchmarkTools.save(joinpath(@__DIR__, "results.json"), results)

# compare against previous results
# TODO: store these results so that we can compare when benchmarking PRs
reference_path = joinpath(@__DIR__, "reference.json")
if ispath(reference_path)
    reference = BenchmarkTools.load(reference_path)[1]
    comparison = judge(minimum(results), minimum(reference))

    println("Improvements:")
    println(improvements(comparison))

    println("Regressions:")
    println(regressions(comparison))
end

## submission

using JSON, HTTP

if real_run
    @info "Submitting to Codespeed..."

    basedata = Dict(
        "branch"        => ENV["CODESPEED_BRANCH"],
        "commitid"      => ENV["CODESPEED_COMMIT"],
        "project"       => ENV["CODESPEED_PROJECT"],
        "environment"   => ENV["CODESPEED_ENVIRONMENT"],
        "executable"    => ENV["CODESPEED_EXECUTABLE"]
    )

    # convert nested groups of benchmark to flat dictionaries of results
    flat_results = []
    function flatten(results, prefix="")
        for (key,value) in results
            if value isa BenchmarkGroup
                flatten(value, "$prefix$key/")
            else
                @assert value isa BenchmarkTools.Trial

                # codespeed reports maxima, but those are often very noisy.
                # get rid of measurements that unnecessarily skew the distribution.
                rmskew!(value)

                push!(flat_results,
                    Dict(basedata...,
                        "benchmark" => "$prefix$key",
                        "result_value" => median(value).time / 1e9,
                        "min" => minimum(value).time / 1e9,
                        "max" => maximum(value).time / 1e9))
            end
        end
    end
    flatten(results)

    HTTP.post("$(ENV["CODESPEED_SERVER"])/result/add/json/",
                ["Content-Type" => "application/x-www-form-urlencoded"],
                HTTP.URIs.escapeuri(Dict("json" => JSON.json(flat_results))))
end


text:
filter_asm(io::IOBuffer) = filter(l -> !startswith(l,";") && !isempty(l) && !startswith(l,"\tnopw") && !startswith(l,"\t.text"), split(String(take!(io)),"\n"))
asm_stat(ss::Array{SubString{String},1}) =
  ( total = length(ss)
  , movs = length(filter(l -> occursin("mov",l),ss))
  , mov  = length(filter(l -> startswith(l,"\tmov"),ss))
  , vmov = length(filter(l -> startswith(l,"\tvmov"),ss))
  )
display_asm_stat_io(io::IOBuffer) = display(asm_stat(filter_asm(io)))
io = IOBuffer()
structinfo(T) = [(fieldoffset(T,i), fieldname(T,i), fieldtype(T,i)) for i = 1:fieldcount(T)];

function mkptr(x::T) where T
  r = Ref{T}(x)
  return (r, Ptr{T}(pointer_from_objref(r)))
end


text:
using DataStructures

abstract type Fixture end
abstract type AbstractGlobalFixture <: Fixture end

function try_collect(val)
    try
        return collect(val)[:]
    catch
        return nothing
    end
end

is_callable(val) = !isempty(methods(val))

"Constant fixture type"
struct ConstantFixture <: AbstractGlobalFixture
    name :: String
    lvals :: Array{LabeledValue, 1}

    """
        ConstantFixture(vals[, labels])

    Create a [`ConstantFixture`](@ref) object.
    Only called to convert a given value to a fixture internally,
    users can just supply iterables directly as testcase or fixture parameters.
    """
    function ConstantFixture(vals, labels=nothing)

        c_vals = try_collect(vals)
        if c_vals === nothing
            error("`vals` must be an iterable")
        end

        if labels === nothing
            labeled_vals = map(LabeledValue, c_vals)
        else
            c_labels = try_collect(labels)
            if c_labels === nothing
                error("`labels` must be an iterable")
            end

            labeled_vals = map(LabeledValue, c_vals, c_labels)
        end

        new(String(gensym("constant")), labeled_vals)
    end
end

setup(fx::ConstantFixture) = fx.lvals

"Global fixture type"
struct GlobalFixture <: AbstractGlobalFixture
    name :: String
    ff :: FixtureFactory
    parameters :: Array{Fixture, 1}
    dependencies :: OrderedSet{AbstractGlobalFixture}

    function GlobalFixture(producer, params...; name=nothing, instant_teardown=false)
        if name === nothing
            name = String(gensym("fixture"))
        end

        if !is_callable(producer)
            error("Producer must be a callable")
        end

        params = collect(map(normalize_fixture, params))
        # TODO: check that it does not depend on any local fixtures
        deps = union(map(dependencies, params)..., global_fixtures(params))
        ff = fixture_factory(producer; instant_teardown=instant_teardown)
        new(name, ff, params, deps)
    end
end

function setup(fx::GlobalFixture, args)
    setup(fx.ff, args)
end

struct RunOptionsFixture <: AbstractGlobalFixture
end

# Shown as non-covered because `__precompile__()` optimizes out the only call
# to this function in `builtin_fixtures.jl`.
function run_options_fixture()
    RunOptionsFixture()
end

"Local fixture type"
struct LocalFixture <: Fixture
    name :: String
    ff :: FixtureFactory
    parameters :: Array{Fixture, 1}
    dependencies :: OrderedSet{AbstractGlobalFixture}

    function LocalFixture(producer, params...; name=nothing)
        if name === nothing
            name = String(gensym("local_fixture"))
        end

        if !is_callable(producer)
            error("Producer must be a callable")
        end

        params = collect(map(normalize_fixture, params))
        deps = union(map(dependencies, params)..., global_fixtures(params))
        ff = fixture_factory(producer; instant_teardown=false)
        new(name, ff, params, deps)
    end

end

is_global_fixture(::GlobalFixture) = true
is_global_fixture(::Any) = false

dependencies(::Fixture) = OrderedSet{AbstractGlobalFixture}()
dependencies(fx::RunOptionsFixture) = OrderedSet{AbstractGlobalFixture}([fx])
dependencies(fx::GlobalFixture) = fx.dependencies
dependencies(fx::LocalFixture) = fx.dependencies

global_fixtures(fxs) = OrderedSet{AbstractGlobalFixture}(filter(is_global_fixture, fxs))

normalize_fixture(f::Fixture) = f
normalize_fixture(f::Pair) = ConstantFixture(f[1], f[2])
normalize_fixture(f) = ConstantFixture(f)

parameters(fixture::LocalFixture) = fixture.parameters
parameters(fixture::GlobalFixture) = fixture.parameters


text:
# Implement number base conversion

function all_your_base(digits, base_in, base_out)

end


text:
# index stuff here
const IndexList = Union{AbstractArray, Base.Generator}

struct IndexSet
    full::IndexList
    reduced::IndexList
    mappings::Dict
end

struct EquivalenceClasses
  pred::Predicate
  indl::IndexList
  classes::Array
end

"""
    EquivalenceClasses(pred::Predicate, indl::IndexList)

Constructs an EquivalenceClasses struct from a predicate `pred`
over `n` indices in the list `indl`.
This inputs requires `n^2` operations in order to create the
adjecency matrix. For large inputs a mapping is therefore
preferabel.

# Examples
```
julia> EquivalenceClasses(Predicate((x,y)->all(x .== -y)),
                          [(i,j) for i in -2:2 for j in 4:7])
```
"""
function EquivalenceClasses(pred::Predicate, indl::IndexList)
  return EquivalenceClasses(pred, indl, find_classes(pred, indl))
end

"""
    EquivalenceClasses(m::Mapping, indl::IndexList)

Constructs an EquivalenceClasses struct from a mapping `m`
over the indices in the list `indl`.
and transitivity will be enforced by the algorithm creating the
adjecency matrix.

# Examples
```
julia> EquivalenceClasses(Predicate((x,y)->all(x .== -y)),
                          [(i,j) for i in -2:2 for j in 4:7])
```
"""
function EquivalenceClasses(pred::Predicate, indl::IndexList)
  return EquivalenceClasses(pred, indl, find_classes(pred, indl))
end

# ====================  Auxiliary functions ====================

Base.show(io::IO, z::EquivalenceClasses) = print(io, "Equivalence classes of ",
                                                 z.pred,
                                                 " over index list ",
                                                 z.indl,
                                                 " with classes ",
                                                 z.classes)


text:
export generate, generate_accompaniment

using NoteSequences
using NoteSequences: DEFAULT_STEPS_PER_SECOND

using StatsBase: wsample

const MIN_PITCH = 21
const MAX_PITCH = 108
const VELOCITY_BINS = 32

"""
    default_melody_encoder()

Return a TextMelodyEncoder with the following parameters.
`minpitch::Int=21`
`maxpitch::Int=109`
`steps_per_second::Int=100`
"""
default_melody_encoder() = TextMelodyEncoder(MIN_PITCH, MAX_PITCH + 1, DEFAULT_STEPS_PER_SECOND)

"""
    default_performance_encoder(add_eos::Bool=true)

Return a MidiPerformanceEncoder with the following parameters.
`steps_per_second::Int`=100
`velocity_bins::Int`=32
`minpitch::Int`=21
`maxpitch::Int`=109

If `add_eos` is true, `EOS_TOKEN` is appended at the end of the sequence.
"""
function default_performance_encoder(add_eos::Bool=true)
    MidiPerformanceEncoder(DEFAULT_STEPS_PER_SECOND, VELOCITY_BINS,
                           MIN_PITCH, MAX_PITCH, add_eos)
end

"""
    generate(model::UnconditionalMusicTransformer;
             primer::Performance=Performance(DEFAULT_STEPS_PER_SECOND, velocity_bins=32),
             numsteps::Int=3000,
             as_notesequence::Bool=false,
             encoder::MidiPerformanceEncoder=default_performance_encoder())

Generate a musical performance from the `UnconditionalMusicTransformer` with the given
`primer`. If the `primer` is empty, a default event of `PerformanceEvent(TIME_SHIFT, 100)` will be
appended. The performance will be converted to and from one-hot indices using the `encoder`.

`numsteps` determines the amount of steps to generate for.
A step is 10 milliseconds or 0.01 seconds at the default steps per second rate of 100.
If `as_notesequence` is true, a `NoteSequence` is returned.
Otherwise, a midifile is returned.
"""
function generate(model::UnconditionalMusicTransformer;
                  primer::Performance=Performance(DEFAULT_STEPS_PER_SECOND, velocity_bins=VELOCITY_BINS),
                  numsteps::Int=3000,
                  as_notesequence::Bool=false,
                  encoder::MidiPerformanceEncoder=default_performance_encoder())

    performance = deepcopy(primer)

    targets = map(event -> encode_event(event, encoder), performance)
    pred = 0
    @info "Generating..."

    # Sample till EOS is encountered or numsteps is reached
    while pred != EOS_TOKEN && performance.numsteps <= numsteps
        t = targets |> device
        logits = model(t)
        probabilities = softmax(logits[:, end]) |> cpu
        pred = wsample(1:encoder.vocab_size, probabilities)

        # Add prediction to the model targets
        push!(targets, pred)

        # Add predicted event to performance
        push!(performance, decode_event(pred, encoder))
    end

    ns = getnotesequence(performance)
    as_notesequence == true && return ns
    midifile(ns)
end

"""
Generate a musical performance from the `UnconditionalMusicTransformer` with the given
`target_indices`. `target_indices` is the one-hot indices of performance events. If it is empty,
a zero matrix will be passed to the model (default token).

`decode_len` determines the length of the sequence to decode.
The generated performance will be decoded from the predicted indices using the `encoder`.
If `as_notesequence` is true, a `NoteSequence` is returned.
Otherwise, a midifile is returned.
"""
function generate(model::UnconditionalMusicTransformer,
                  target_indices::Vector{Int};
                  decode_len::Int=1024,
                  as_notesequence::Bool=false,
                  encoder::MidiPerformanceEncoder=default_performance_encoder())

    targets = deepcopy(target_indices)
    pred = 0
    @info "Generating..."

    # Sample till EOS is encountered or decoding length is reached
    while pred != EOS_TOKEN && length(targets) <= decode_len
        t = targets |> device
        logits = model(t)
        probabilities = softmax(logits[:, end]) |> cpu
        pred = wsample(1:encoder.vocab_size, probabilities)
        push!(targets, pred)
    end

    ns = decode_to_notesequence(targets, encoder)
    as_notesequence == true && return ns
    midifile(ns)
end

"""
    generate_accompaniment(model::MelodyConditionedMusicTransformer,
                           melody::Vector{Int};
                           decode_len::Int=2048,
                           as_notesequence::Bool=false,
                           target_encoder::MidiPerformanceEncoder=default_performance_encoder())

Generate an accompaniment for the given melody using the `MelodyConditionedMusicTransformer`.

`decode_len` determines the length of the sequence to decode.
The generated performance will be decoded from the predicted indices using `target_encoder`.
If `as_notesequence` is true, a `NoteSequence` is returned.
Otherwise, a midifile is returned.
"""
function generate_accompaniment(model::MelodyConditionedMusicTransformer,
                                melody::Vector{Int};
                                decode_len::Int=(4096 - length(inputs)),
                                as_notesequence::Bool=false,
                                target_encoder::MidiPerformanceEncoder=default_performance_encoder())

    if isempty(melody)
        throw(error("Melody must not be empty"))
    end

    targets = Vector{Int}()
    @info "Generating..."

    model = model |> device

    melody, targets = device.([melody, targets])
    logits, encoder_context = model(melody, targets)
    pred = wsample(1:target_encoder.vocab_size, softmax(logits[:, end]) |> cpu)
    push!(targets, pred)

    while pred != EOS_TOKEN && length(targets) <= decode_len
        t = targets |> device
        logits = model(t, encoder_context)
        probabilities = softmax(logits[:, end]) |> cpu
        pred = wsample(1:target_encoder.vocab_size, probabilities)
        push!(targets, pred)
    end

    ns = decode_to_notesequence(targets, target_encoder)
    as_notesequence == true && return ns
    midifile(ns)
end


text:
__precompile__()
module Derivatives

using InplaceRealFFT
using LESFilter: rfftfreq, fftfreq

export dx, dy, dz, dx!, dy!, dz!

function loopdx!(fieldhat::AbstractArray{<:Complex,3}, kim::AbstractVector)
  nx,ny,nz = size(fieldhat)
  let nx = nx, ny = ny, nz = nz
  Threads.@threads for l = 1:nz
    for j = 1:ny
      @simd for i = 1:nx
      @inbounds fieldhat[i,j,l] = fieldhat[i,j,l]*kim[i]
      end
    end
  end
  end
end

function dx(field::AbstractArray{T,3},len::Real,n::Integer=1) where T<:Union{Float32,Float64}
  nx,ny,nz = size(field)

  fieldhat = rfft(field,1)
  kim = (2π .* rfftfreq(nx,len) .* im) .^ n

  loopdx!(fieldhat,kim)

  return irfft(fieldhat,nx,1)
end

function loopdy!(fieldhat::AbstractArray{<:Complex,3}, kim::AbstractVector)
  nx,ny,nz = size(fieldhat)
  let nx = nx, ny = ny, nz = nz
  Threads.@threads for l = 1:nz
    for j = 1:ny
      @simd for i = 1:nx
        @inbounds fieldhat[i,j,l] = fieldhat[i,j,l]*kim[j]
      end
    end
  end
  end
end

function dy(field::AbstractArray{T,3},len::Real,n::Integer=1) where T<:Union{Float32,Float64}
  nx,ny,nz = size(field)

  fieldhat = rfft(field,2)

  kim = (2π .* rfftfreq(ny,len) .* im) .^ n

  loopdy!(fieldhat,kim)

  return irfft(fieldhat,ny,2)
end

function loopdz!(fieldhat::AbstractArray{<:Complex,3}, kim::AbstractVector)
  nx,ny,nz = size(fieldhat)
  let nx = nx, ny = ny, nz = nz
  Threads.@threads for l = 1:nz
    for j = 1:ny
      @simd for i = 1:nx
        @inbounds fieldhat[i,j,l] = fieldhat[i,j,l]*kim[l]
      end
    end
  end
  end
end

function dz(field::AbstractArray{T,3},len::Real,n::Integer=1) where T<:Union{Float32,Float64}
  nx,ny,nz = size(field)

  fieldhat = rfft(field,3)

  kim = (2π .* rfftfreq(nz,len) .* im) .^ n

  loopdz!(fieldhat,kim)

  return irfft(fieldhat,nz,3)
end

function dx!(field::AbstractPaddedArray{T,3},len::Real,n::Integer=1) where {T<:Union{Float64,Float32}}

  nx,ny,nz = size(real(field))

  rfft!(field,1)

  kim = (2π .* rfftfreq(nx,len) .* im) .^ n

  loopdx!(InplaceRealFFT.complex_view(field),kim)

  return irfft!(field,1)
end

function dx(field::AbstractPaddedArray{T,3},len::Real,n::Integer=1) where {T<:Union{Float64,Float32}}
  fieldhat = copy(field)
  dx!(fieldhat,len,n)
  return fieldhat
end

function dy!(field::AbstractPaddedArray{T,3},len::Real,n::Integer=1) where {T<:Union{Float64,Float32}}

  nx,ny,nz = size(real(field))
  rfft!(field,1:2)

  kim = (2π .* fftfreq(ny,len) .* im) .^ n

  loopdy!(InplaceRealFFT.complex_view(field),kim)

  return irfft!(field,1:2)
end

function dy(field::AbstractPaddedArray{T,3},len::Real,n::Integer=1) where {T<:Union{Float64,Float32}}
  fieldhat = copy(field)
  dy!(fieldhat,len,n)
  return fieldhat
end

function dz!(field::AbstractPaddedArray{T,3},len::Real,n::Integer=1) where {T<:Union{Float64,Float32}}

  nx,ny,nz = size(real(field))
  rfft!(field,(1,3))

  kim = (2π .* fftfreq(nz,len) .* im) .^ n

  loopdz!(InplaceRealFFT.complex_view(field),kim)

  return irfft!(field,(1,3))
end

function dz(field::AbstractPaddedArray{T,3},len::Real,n::Integer=1) where {T<:Union{Float64,Float32}}
  fieldhat = copy(field)
  dz!(fieldhat,len,n)
  return fieldhat
end

end # module


text:
function pseudos(server, pseudoset, atsyms::Vector{Symbol}, fuzzy = "")
    s = Servers.maybe_start_server(server)
    pseudo_paths = list_pseudoset(s, pseudoset, fuzzy)
    ps = Dict{Symbol,String}()
    for a in atsyms
        path = get(pseudo_paths, a, nothing)
        if path === nothing
            error("No pseudo for atom $a found in set $pseudoset.")
        end
        t = tempname()
        Servers.pull(s, path, t)
        ps[a] = read(t, String)
        rm(t)
    end
    return ps
end

function list_pseudoset(server, pseudoset, fuzzy="")
    s = Servers.maybe_start_server(server)
    resp = HTTP.get(s, "/pseudos/$pseudoset", [], JSON3.write(fuzzy))
    if resp.status == 204
        error("No pseudoset $pseudoset found on Server $(s.name). Please first configure it using configure_pseudoset.")
    end
    return JSON3.read(resp.body, Dict{Symbol,String})
end

"""
    list_pseudosets(server = "localhost")

Lists the pseudosets that have previously been set up.
"""
function list_pseudosets(server = "localhost")
    s = Servers.maybe_start_server(server)
    return JSON3.read(HTTP.get(s, "/pseudos").body, Vector{String})
end

"""
    configure_pseudoset(set_name::String, dir::String, server = "localhost")

Reads the specified `dir` and sets up the pseudos for `set`.
"""
function configure_pseudoset(set_name::String, dir::String, server = "localhost")
    s = Servers.maybe_start_server(server)
    p = isabspath(dir) ? dir : joinpath(s, dir)
    n_pseudos = JSON3.read(HTTP.post(s, "/configure_pseudoset/" * p, [],
                                     JSON3.write(set_name)).body, Int)
    @info "Configured $n_pseudos pseudos on Server $(s.name), found in dir $p."
end

"""
    rm_pseudoset!(set_name::String, server = "localhost")

Removes the pseudo set from the server.
"""
function rm_pseudoset!(set_name::String, server = "localhost")
    s = Servers.maybe_start_server(server)
    return HTTP.put(s, "/rm_pseudos", [], JSON3.write(set_name))
end
#---#

"""
    set_pseudos!(job::Job, set::Symbol, specifier::String=""; kwargs...)
    set_pseudos!(structure::Structure, set::Symbol, specifier::String=""; kwargs...)
    set_pseudos!(job::Job, atsym::Symbol, set::Symbol, specifier::String=""; kwargs...)
    set_pseudos!(structure::Structure, atsym::Symbol, set::Symbol, specifier::String=""; kwargs...)

Sets the pseudopotentials of the atoms inside the `structure` (or `job.structure`) to the ones of `set`.
`specifier` can be specified to select a specific pseudo if multiple pseudopotentials
for a given element exist in the set.
Example:
```
set_pseudos!(job, :pbesol, "rrkjus")
```
will select the pseudo file that contains "rrkjus" in the filename.

If `atsym` is used, only the pseudos of the atoms with that name will be set.

    set_pseudos!(job::Job, at_pseudos::Pair{Symbol, Pseudo}...; kwargs...)
    set_pseudos!(structure::Structure, at_pseudos::Pair{Symbol, Pseudo}...; kwargs...)

Convenience function that allows to set pseudopotentials for multiple atom types at the same time.
e.g. `set_pseudos!(job, :Si => getdefault_pseudo(:Si, :sssp)
"""
function Structures.set_pseudos!(job::Job, set, server=job.server; specifier::String = "", kwargs...)
    atsyms = unique(map(x -> x.element.symbol, job.structure.atoms))
    return Structures.set_pseudos!(job.structure, pseudos(server, set, atsyms, specifier);
                                   kwargs...)
end