text:
#--------------------
#--------------------
#### Symbolic
#--------------------
abstract type Symbolic{T} end

### Interface to be defined for `simplify` to work:

"""
    istree(x::T)

Check if `x` represents an expression tree. If returns true,
it will be assumed that `operation(::T)` and `arguments(::T)`
methods are defined. Definining these three should allow use
of `simplify` on custom types. Optionally `symtype(x)` can be
defined to return the expected type of the symbolic expression.
"""
istree(x) = false

"""
    operation(x::T)

Returns the operation (a function object) performed by an expression
tree. Called only if `istree(::T)` is true. Part of the API required
for `simplify` to work. Other required methods are `arguments` and `istree`
"""
function operation end

"""
    arguments(x::T)

Returns the arguments (a `Vector`) for an expression tree.
Called only if `istree(x)` is `true`. Part of the API required
for `simplify` to work. Other required methods are `operation` and `istree`
"""
function arguments end

"""
    symtype(x)

The supposed type of values in the domain of x. Tracing tools can use this type to
pick the right method to run or analyse code.

This defaults to `typeof(x)` if `x` is numeric, or `Any` otherwise.
For the types defined in this package, namely `T<:Symbolic{S}` it is `S`.

Define this for your symbolic types if you want `simplify` to apply rules
specific to numbers (such as commutativity of multiplication). Or such
rules that may be implemented in the future.
"""
function symtype end

symtype(x::Number) = typeof(x)

symtype(x) = Any

symtype(::Symbolic{T}) where {T} = T

Base.isequal(s::Symbolic, x) = false
Base.isequal(x, s::Symbolic) = false
Base.isequal(x::Symbolic, y::Symbolic) = false
### End of interface

"""
    to_symbolic(x)

Convert `x` to a `Symbolic` type, using the `istree`, `operation`, `arguments`,
and optionally `symtype` if available.
"""
to_symbolic(x::Symbolic) = x

function to_symbolic(x)
    if !istree(x)
        return x
    end

    op = to_symbolic(operation(x))

    if symtype(x) === Any
        Term(op, map(to_symbolic, arguments(x)))
    else
        Term{symtype(x)}(op, map(to_symbolic, arguments(x)))
    end
end

Base.one( s::Symbolic) = one( symtype(s))

Base.zero(s::Symbolic) = zero(symtype(s))

"""
    promote_symtype(f, Ts...)

The result of applying `f` to arguments of [`symtype`](#symtype) `Ts...`

```julia
julia> promote_symtype(+, Real, Real)
Real

julia> promote_symtype(+, Complex, Real)
Number

julia> @syms f(x)::Complex
(f(::Number)::Complex,)

julia> promote_symtype(f, Number)
Complex
```

When constructing [`Term`](#Term)s without an explicit symtype,
`promote_symtype` is used to figure out the symtype of the Term.
"""
promote_symtype(f, Ts...) = Any

#--------------------
#--------------------
#### Syms
#--------------------
"""
    Sym{T}(name::Symbol)

A named variable of type `T`. Type `T` can be `FnType{X,Y}` which
means the variable is a function with the type signature X -> Y where
`X` is a tuple type of arguments and `Y` is any type.
"""
struct Sym{T} <: Symbolic{T}
    name::Symbol
end

const Variable = Sym # old name
Sym(x) = Sym{symtype(x)}(x)

Base.nameof(v::Sym) = v.name

Base.isequal(v1::Sym{T}, v2::Sym{T}) where {T} = v1 === v2

Base.show(io::IO, v::Sym) = print(io, v.name)

#---------------------------
#---------------------------
#### Function-like variables
#---------------------------

# Maybe don't even need a new type, can just use Sym{FnType}
struct FnType{X<:Tuple,Y} end

(f::Sym{<:FnType})(args...) = Term{promote_symtype(f, symtype.(args)...)}(f, [args...])

function (f::Sym)(args...)
    error("Sym $f is not callable. " *
          "Use @syms $f(var1, var2,...) to create it as a callable. " *
          "See ?@fun for more options")
end

"""
`promote_symtype(f::Sym{FnType{X,Y}}, arg_symtypes...)`

The output symtype of applying variable `f` to arugments of symtype `arg_symtypes...`.
if the arguments are of the wrong type then this function will error.
"""
function promote_symtype(f::Sym{FnType{X,Y}}, args...) where {X, Y}
    nrequired = fieldcount(X)
    ngiven    = nfields(args)

    if nrequired !== ngiven
        error("$f takes $nrequired arguments; $ngiven arguments given")
    end

    for i in 1:ngiven
        t = X.parameters[i]
        if !(args[i] <: t)
            error("Argument to $f at position $i must be of symbolic type $t")
        end
    end
    return Y
end

"""
    @syms <lhs_expr>[::T1] <lhs_expr>[::T2]...

For instance:

    @syms foo::Real bar baz(x, y::Real)::Complex

Create one or more variables. `<lhs_expr>` can be just a symbol in which case
it will be the name of the variable, or a function call in which case a function-like
variable which has the same name as the function being called. The Sym type, or
in the case of a function-like Sym, the output type of calling the function
can be set using the `::T` syntax.

# Examples:

- `@syms foo bar::Real baz::Int` will create
variable `foo` of symtype `Number` (the default), `bar` of symtype `Real`
and `baz` of symtype `Int`
- `@syms f(x) g(y::Real, x)::Int h(a::Int, f(b))` creates 1-arg `f` 2-arg `g`
and 2 arg `f`. The second argument to `h` must be a one argument function-like
variable. So, `h(1, g)` will fail and `h(1, f)` will work.
"""
macro syms(xs...)
    defs = map(xs) do x
        n, t = _name_type(x)
        :($(esc(n)) = Sym{$(esc(t))}($(Expr(:quote, n))))
    end

    Expr(:block, defs...,
         :(tuple($(map(x->esc(_name_type(x).name), xs)...))))
end

function syms_syntax_error()
    error("Incorrect @syms syntax. Try `@syms x::Real y::Complex g(a) f(::Real)::Real` for instance.")
end

function _name_type(x)
    if x isa Symbol
        return (name=x, type=Number)
    elseif x isa Expr && x.head === :(::)
        if length(x.args) == 1
            return (name=nothing, type=x.args[1])
        end
        lhs, rhs = x.args[1:2]
        if lhs isa Expr && lhs.head === :call
            # e.g. f(::Real)::Unreal
            type = map(x->_name_type(x).type, lhs.args[2:end])
            return (name=lhs.args[1], type=:($FnType{Tuple{$(type...)}, $rhs}))
        else
            return (name=lhs, type=rhs)
        end
    elseif x isa Expr && x.head === :call
        return _name_type(:($x::Number))
    else
        syms_syntax_error()
    end
end

function Base.show(io::IO, f::Sym{<:FnType{X,Y}}) where {X,Y}
    print(io, f.name)
    argrepr = join(map(t->"::"*string(t), X.parameters), ", ")
    print(io, "(", argrepr, ")")
    print(io, "::", Y)
end

#--------------------
#--------------------
#### Terms
#--------------------
"""
    Term{T}(f, args::AbstractArray)

or
    Term(f, args::AbstractArray)

Symbolic expression representing the result of calling `f(args...)`.

- `operation(t::Term)` returns `f`
- `arguments(t::Term)` returns `args`
- `symtype(t::Term)` returns `T`

If `T` is not provided during construction, it is queried by calling
`SymbolicUtils.promote_symtype(f, map(symtype, args)...)`.

See [promote_symtype](#promote_symtype)
"""
struct Term{T} <: Symbolic{T}
    f::Any
    arguments::Any
end

istree(t::Term) = true

Term(f, args) = Term{rec_promote_symtype(f, map(symtype, args)...)}(f, args)

operation(x::Term) = x.f

arguments(x::Term) = x.arguments

function Base.hash(t::Term{T}, salt::UInt) where {T}
    hash(arguments(t), hash(operation(t), hash(T, salt)))
end

function Base.isequal(t1::Term, t2::Term)
    t1 === t2 && return true
    a1 = arguments(t1)
    a2 = arguments(t2)

    isequal(operation(t1), operation(t2)) && length(a1) == length(a2) &&
        all(isequal(l,r) for (l, r) in zip(a1,a2))
end

function term(f, args...; type = nothing)
    if type === nothing
        T = rec_promote_symtype(f, symtype.(args)...)
    else
        T = type
    end
    Term{T}(f, [args...])
end

node_count(t) = istree(t) ? reduce(+, node_count(x) for x in  arguments(t), init=0) + 1 : 1

#--------------------
#--------------------
####  Pretty printing
#--------------------
const show_simplified = Ref(false)

function Base.show(io::IO, t::Term)
    if get(io, :simplify, show_simplified[])
        s = simplify(t)
        color = isequal(s, t) ? :white : :yellow

        Base.printstyled(IOContext(io, :simplify=>false, :withcolor=>color),
                         s, color=color)
    else
        f = operation(t)
        args = arguments(t)
        fname = nameof(f)
        binary = Base.isbinaryoperator(fname)
        color = get(io, :withcolor, :white)
        if binary
            get(io, :paren, false) && Base.printstyled(io, "(",color=color)
            for i = 1:length(args)
                length(args) == 1 && Base.printstyled(io, fname, color=color)

                args[i] isa Complex && Base.printstyled(io, "(",color=color)
                Base.printstyled(IOContext(io, :paren => true),
                                 args[i], color=color)
                args[i] isa Complex && Base.printstyled(io, ")",color=color)

                i != length(args) && Base.printstyled(io, " $fname ", color=color)
            end
            get(io, :paren, false) && Base.printstyled(io, ")", color=color)
        else
            Base.printstyled(io, "$fname(", color=color)
            for i=1:length(args)
                Base.printstyled(IOContext(io, :paren => false),
                                 args[i], color=color)
                i != length(args) && Base.printstyled(io, ", ", color=color)
            end
            Base.printstyled(io, ")", color=color)
        end
    end
    return nothing
end

showraw(io, t) = Base.show(IOContext(io, :simplify=>false), t)
showraw(t) = showraw(stdout, t)


text:
"""
    NumaAllocators.Windows

NUMA support for Windows.

See also https://docs.microsoft.com/en-us/windows/win32/api/memoryapi/nf-memoryapi-virtualallocexnuma
"""
module Windows

using ..NumaAllocators: AbstractNumaAllocator
import ArrayAllocators: allocate, iszeroinit
using ArrayAllocators.ByteCalculators: nbytes
using ArrayAllocators.Windows: wrap_virtual, hCurrentProcess, MEM_COMMIT_RESERVE, PAGE_READWRITE, kernel32

function GetNumaProcessorNode(processor = GetCurrentProcessorNumber())
    node_number = Ref{Cuchar}(0)
    status = ccall((:GetNumaProcessorNode, kernel32), Cint, (Cuchar, Ptr{Cuchar}), processor, node_number)
    if status == 0
        error("Could not retrieve NUMA node for processor $processor.")
    end
    return node_number[]
end

# https://docs.microsoft.com/en-us/windows/win32/api/systemtopologyapi/nf-systemtopologyapi-getnumahighestnodenumber
function GetNumaHighestNodeNumber()
    node_number = Ref{Culong}(0)
    status = ccall((:GetNumaHighestNodeNumber, kernel32), Cint, (Ptr{Culong},), node_number)
    if status == 0
        error("Could not retrieve highest NUMA node.")
    end
    return node_number[]
end

function GetCurrentProcessorNumber()
    return ccall((:GetCurrentProcessorNumber, "kernel32"), Cint, ())
end

#=
LPVOID VirtualAllocExNuma(
  [in]           HANDLE hProcess,
  [in, optional] LPVOID lpAddress,
  [in]           SIZE_T dwSize,
  [in]           DWORD  flAllocationType,
  [in]           DWORD  flProtect,
  [in]           DWORD  nndPreferred
);
=#

function VirtualAllocExNuma(hProcess, lpAddress, dwSize, flAllocationType, flProtect, nndPreferred)
    ccall((:VirtualAllocExNuma, kernel32),
        Ptr{Nothing}, (Ptr{Nothing}, Ptr{Nothing}, Csize_t, Culong, Culong, Culong),
        hProcess, lpAddress, dwSize, flAllocationType, flProtect, nndPreferred)
end
function VirtualAllocExNuma(dest_size, numa_node)
    VirtualAllocExNuma(hCurrentProcess, C_NULL, dest_size, MEM_COMMIT_RESERVE, PAGE_READWRITE, numa_node)
end

abstract type AbstractWinNumaAllocator{B} <: AbstractNumaAllocator{B} end

"""
    WinNumaAllocator

Allocate memory on a specific NUMA node with `VirtualAllocExNuma`.

See also https://docs.microsoft.com/en-us/windows/win32/api/memoryapi/nf-memoryapi-virtualallocexnuma
"""
struct WinNumaAllocator{B} <: AbstractWinNumaAllocator{B}
    node::Int
end

function allocate(n::WinNumaAllocator, num_bytes)
    return VirtualAllocExNuma(num_bytes, n.node)
end
Base.unsafe_wrap(::WinNumaAllocator, args...) = wrap_virtual(args...)
iszeroinit(::Type{A}) where A <: WinNumaAllocator = true

end # module Windows

ext: jl


text:
printstyled("  scan_seed\n", color=:light_green)

fname = string(path, "/SampleFiles/SEED/test.mseed")
test_sac_file = string(path, "/SampleFiles/SAC/test_le.sac")
@test_throws ErrorException scan_seed(test_sac_file)

redirect_stdout(out) do
  soh = scan_seed(fname, v=3)
  S = read_data("mseed", fname)
  h = String.(strip.(split(soh[1], ",")))

  nx = parse(Int64, split(h[2], "=")[2])
  @test nx == length(S.x[1])

  ng = parse(Int64, split(h[3], "=")[2])
  @test ng == 0

  nfs = parse(Int64, split(h[4], "=")[2])
  @test nfs == 1
end

if has_restricted
  redirect_stdout(out) do
    fname = path * "/SampleFiles/Restricted/Steim2-AllDifferences-BE.mseed"
    soh = scan_seed(fname, quiet=true)
    h = String.(strip.(split(soh[1], ",")))

    nx = parse(Int64, split(h[2], "=")[2])
    @test nx == 3096

    ng = parse(Int64, split(h[3], "=")[2])
    @test ng == 0

    nfs = parse(Int64, split(h[4], "=")[2])
    @test nfs == 1

    fname = path * "/SampleFiles/Restricted/SHW.UW.mseed"
    scan_seed(fname, fs_times=true)
    scan_seed(fname, seg_times=true)
    soh = scan_seed(fname, quiet=true)
    S = read_data("mseed", fname)
    nfs_expect = [143, 7]

    for i in 1:S.n
      h = String.(strip.(split(soh[i], ",")))
      nx = parse(Int64, split(h[2], "=")[2])
      @test nx == length(S.x[i])

      # This is occasionally off-by-one from the true total
      ng = parse(Int64, split(h[3], "=")[2])
      @test abs(ng - (size(S.t[i], 1) - 2)) ≤ 1

      nfs = parse(Int64, split(h[4], "=")[2])
      @test nfs == nfs_expect[i]
    end
  end
end


text:
module NeuroGraphs

using ArrayInterface

export AdjMap, AdjDiMap, WAdjMap, WAdjDiMap

include("AdjacencyGraphs/AdjacencyGraphs.jl")
using .AdjacencyGraphs

end


text:
using BlackBoxOptim

function fitness_for_opt(problem, numDimensions, populationSize, numSteps, method)
    problem = instantiate(problem, numDimensions)

    println("\n$(problem.name), n = $(numdims(problem)), optimizer = $(string(method))")

    res = bboptimize(problem; Method = method,
        NumDimensions = numDimensions,
        PopulationSize = populationSize,
        MaxSteps = numSteps
    )

    return best_fitness(res)
end


text:
using NiLangCore, Test

@testset "stack" begin
    for (stack, x) in [
                (FLOAT64_STACK, 0.3), (FLOAT32_STACK, 0f4),
                (INT64_STACK, 3), (INT32_STACK, Int32(3)),
                (COMPLEXF64_STACK, 4.0+0.3im), (COMPLEXF32_STACK, 4f0+0f3im),
                (BOOL_STACK, true),
                ]
        println(stack)
        push!(stack, x)
        @test pop!(stack) === x
    end
end

@testset "stack operations" begin
    z = 1.0
    NiLangCore.empty_global_stacks!()
    @test_throws ArgumentError (@instr GLOBAL_STACK[end] ↔ y::∅)
    y = 4.0
    @test_throws ArgumentError (@instr GLOBAL_STACK[end] → y)
    @test_throws BoundsError (@instr @invcheckoff GLOBAL_STACK[end] ↔ y)
    @test_throws ArgumentError (@instr @invcheckoff GLOBAL_STACK[end] → y)
    x = 0.3
    NiLangCore.empty_global_stacks!()
    @instr GLOBAL_STACK[end+1] ↔ x
    @instr GLOBAL_STACK[end] ↔ x::∅
    @test x === 0.3
    @instr @invcheckoff GLOBAL_STACK[end+1] ↔ x
    y = 0.5
    @instr GLOBAL_STACK[end+1] ↔ y
    @instr @invcheckoff GLOBAL_STACK[end] ↔ x::∅
    @test x == 0.5

    x =0.3
    st = Float64[]
    @instr st[end+1] ↔ x
    @test length(st) == 1
    @instr st[end] ↔ x::∅
    @test length(st) == 0
    @test x === 0.3
    @instr st[end+1] ↔ x
    @test length(st) == 1

    y = 0.5
    @instr st[end+1] ↔ y
    @instr @invcheckoff st[end] ↔ x::∅
    @test x == 0.5

    @i function test(x)
        x2 ← zero(x)
        x2 += x^2
        GLOBAL_STACK[end+1] ↔ x
        x::∅ ↔ x2
    end
    @test test(3.0) == 9.0
    l = length(NiLangCore.GLOBAL_STACK)
    @test check_inv(test, (3.0,))
    @test length(NiLangCore.GLOBAL_STACK) == l

    @i function test2(x)
        x2 ← zero(x)
        x2 += x^2
        @invcheckoff GLOBAL_STACK[end+1] ↔ x
        x::∅ ↔ x2
    end
    @test test2(3.0) == 9.0
    l = length(NiLangCore.GLOBAL_STACK)
    @test check_inv(test2, (3.0,))
    @test length(NiLangCore.GLOBAL_STACK) == l

    x = 3.0
    @instr GLOBAL_STACK[end+1] ↔ x
    NiLangCore.empty_global_stacks!()
    l = length(NiLangCore.GLOBAL_STACK)
    @test l == 0
end

@testset "copied push/pop stack operations" begin
    NiLangCore.empty_global_stacks!()
    x =0.3
    @instr GLOBAL_STACK[end+1] ← x
    @test x === 0.3
    @instr GLOBAL_STACK[end] → x
    @test x === 0.3
    @instr GLOBAL_STACK[end+1] ← x
    x = 0.4
    @test_throws InvertibilityError @instr GLOBAL_STACK[end] → x
    y = 0.5
    @instr GLOBAL_STACK[end+1] ← y
    @instr @invcheckoff GLOBAL_STACK[end] → x
    @test x == 0.5

    st = []
    x = [0.3]
    @instr st[end+1] ← x
    @test st[1] !== [0.3]
    @test st[1] ≈ [0.3]

    x =0.3
    st = Float64[]
    @instr ~(st[end] → x)
    @test x === 0.3
    @test length(st) == 1
    @instr ~(st[end+1] ← x)
    @test length(st) == 0
    @test x === 0.3
    @instr @invcheckoff st[end+1] ← x
    @test length(st) == 1
    x = 0.4
    @test_throws InvertibilityError @instr st[end] → x
    @test length(st) == 0

    y = 0.5
    @instr st[end+1] ← y
    @instr @invcheckoff st[end] → x
    @test x == 0.5

    @i function test(x, x2)
        x2 += x^2
        GLOBAL_STACK[end+1] ← x
        x ↔ x2
    end
    @test test(3.0, 0.0) == (9.0, 3.0)
    l = length(NiLangCore.GLOBAL_STACK)
    @test check_inv(test, (3.0, 0.0))
    @test length(NiLangCore.GLOBAL_STACK) == l
end

@testset "dictionary & vector" begin
    # allocate and deallocate
    @i function f1(d, y)
        d["y"] ← y
    end
    d = Dict("x" => 34)
    @test f1(d, 3) == (Dict("x"=>34, "y"=>3), 3)
    @test_throws InvertibilityError f1(d, 3)
    d = Dict("x" => 34)
    @test check_inv(f1, (d, 3))
    # not available on vectors
    @i function f2(d, y)
        d[2] ← y
    end
    @test_throws MethodError f2([1,2,3], 3)

    # swap
    @i function f3(d, y)
        d["y"] ↔ y
    end
    d = Dict("y" => 34)
    @test f3(d, 3) == (Dict("y"=>3), 34)
    d = Dict("z" => 34)
    @test_throws KeyError f3(d, 3)
    d = Dict("y" => 34)
    @test check_inv(f3, (d, 3))
    # swap on vector
    @i function f4(d, y, x)
        d[2] ↔ y
        d[end] ↔ x
    end
    d = [11,12,13]
    @test f4(d, 1,2) == ([11,1,2],12,13)
    d = [11,12,13]
    @test check_inv(f4, (d, 1,2))

    # swap to empty
    @i function f5(d, x::T) where T
        d["x"]::∅ ↔ x   # swap in
        d["y"] ↔ x::∅  # swap out
    end
    d = Dict("y" => 34)
    @test f5(d, 3) == (Dict("x"=>3), 34)
    d = Dict("y" => 34)
    @test check_inv(f5, (d, 3))
    d = Dict("x" => 34)
    @test_throws InvertibilityError f5(d, 3)
    # not available on vectors
    @i function f6(d, y)
        d[2]::∅ ↔ y
    end
    @test_throws MethodError f6([1,2,3], 3)
end

@testset "inverse stack" begin
    @i function f(x)
        x[end+1] ← 1
    end

    x = FastStack{Int}(3)
    @test check_inv(f, (x,))
end


text:
@testset "Basic ser/deser example" begin
    using ONNXNaiveNASflux, Test, Statistics

    l1 = Conv((3,3), 2=>3, relu)
    l2 = Dense(3, 4, elu)

    f = function(x,y)
        x = l1(x)
        # Home-brewed global average pool
        x = dropdims(mean(x, dims=(1,2)), dims=(1,2))
        x = l2(x)
        return x + y
    end

    # Serialize f. Use a string to save to file instead
    io = PipeBuffer()
    x_shape = (:W, :H, 2, :Batch)
    y_shape = (4, :Batch)
    save(io, f, x_shape, y_shape)

    # Deserialize as a NaiveNASflux CompGraph (or use the deserialization method from ONNX.jl)
    g = load(io)

    x = ones(Float32, 5,4,2,3)
    y = ones(Float32, 4, 3)
    @test g(x,y) ≈ f(x,y)

    # Serialization of CompGraphs does not require input shapes to be provided as they can be inferred.
    io = PipeBuffer()
    save(io, g)

    g = load(io)
    @test g(x,y) ≈ f(x,y)

    @test name.(vertices(g)) == ["data_0", "conv_0", "reducemean_0", "squeeze_0", "dense_0", "data_1", "add_0"]
    @test nout.(vertices(g)) == [2, 3, 3, 3, 4, 4, 4]
end

@testset "Check supported OPs list" begin
    # I'm sure there is a better way to do this...
    sosec = "## Supported Operations"
    bts = "```"
    want = sosec
    add = false
    found = String[]
    for line in eachline(joinpath(dirname(pathof(ONNXNaiveNASflux)),"..", "README.md"))
        add && line == bts && break
        add && !isempty(line) && push!(found, line)
        if line == sosec
            want = bts
        end
        if want == bts && line == bts
            add =true
        end
    end

    io = IOBuffer()
    ONNXNaiveNASflux.list_supported_ops(io)

    actual = split(String(take!(io)), "\n"; keepempty=false);

    @test found == actual
end

@testset "Custom serialization example" begin

    import ONNXNaiveNASflux: AbstractProbe, recursename, nextname, newfrom, add!, name, ONNX
    function myfun(probes::AbstractProbe...)
        p = probes[1] # select any probe
        optype = "MyOpType"
        # Naming strategy (e.g. how to avoid duplicate names) is provided by the probe
        # Not strictly needed, but the onnx model is basically corrupt if duplicates exist
        nodename = recursename(optype, nextname(p))

        # Add ONNX node info
        add!(p, ONNX.NodeProto(
        # Names of input is provided by probes. This is why new probes need to be provided as output
        input = collect(name.(probes)),
        # Name of output from this node
        output = [nodename],
        op_type = optype))

        # Probes can procreate like this
        return newfrom(p, nodename, s -> s)
    end

    gp = ONNXNaiveNASflux.graphproto()
    pps = ONNXNaiveNASflux.inputprotoprobe!.(Ref(gp), ("in1", "in2"), ((1,2), (3,4)), s -> "out")
    pp = myfun(pps...)

    @test name(pp) == "out"
    @test length(gp.node) == 1
    @test gp.node[1].input == ["in1", "in2"]
    @test gp.node[1].output == ["out"]
    @test gp.node[1].op_type == "MyOpType"

end

@testset "Custom deserialization example" begin
    import ONNXNaiveNASflux: actfuns

    # All inputs which are not output from another node in the graph are provided in the method call
    actfuns[:MyActFun1] = (params, α, β) -> x -> x^α + β
    # Params contains a Dict with attributes.
    actfuns[:MyActFun2] = function(params)
        α = get(params, :alpha, 1)
        return x -> α / x
    end
    ONNXNaiveNASflux.refresh()

    af1 = ONNXNaiveNASflux.verts[:MyActFun1]("af1", [inputvertex("in", 3)], Dict(), 2, 3)
    @test af1(4) == 4^2 + 3

    af2 = ONNXNaiveNASflux.verts[:MyActFun2]("af2", [inputvertex("in", 3)], Dict(:alpha => 3))
    @test af2(4) == 3/4

    delete!(ONNXNaiveNASflux.verts, :MyActFun1)
    delete!(ONNXNaiveNASflux.verts, :MyActFun2)
    delete!(ONNXNaiveNASflux.invariantops, :MyActFun1)
    delete!(ONNXNaiveNASflux.invariantops, :MyActFun2)
    delete!(actfuns, :MyActFun1)
    delete!(actfuns, :MyActFun2)
    ONNXNaiveNASflux.refresh()
end


text:
module GtkJulia

function julia_main()
    try
        real_main()
    catch
        Base.invokelatest(Base.display_error, Base.catch_stack())
        return 1
    end
    return 0
end

function real_main()
    println("Hello World!")
end

end #module GtkJulia

max_stars_repo_name: Lecrapouille/BacASable

max_issues_repo_name: Lecrapouille/BacASable

max_forks_repo_name: Lecrapouille/BacASable


text:
###################################################
#
# A script to generate a list of all packages in reST format
#
# The JSON and Calendar packages are prerequisites for this script, and should be installed
#
# This script also needs a Github OAuth token, since un-authenticated requests are limited to 60 per hour
# Therefore, you need to specify the token as an environment variable GH_AUTH before running this program
#
# Once the auth token is specified, run `make listpkg` in this directory to create the latest listing.
# Note that this will update your packages.
#
# Checkin the resultant `packagelist.rst`, to allow the docs to be built locally without running this script.
#
# **Note**:
#
# To create a github oauth token, you will need to create a registered github application. Once you have an application, you can then call
#
#     curl -i -u <gh-user> -d @gh_auth https://api.github.com/authorizations
#
# This will prompt for your gh password using basic authentication. `gh_auth` is a local file containing:
#
# ```JSON
# {
#   "scopes": ["public_repo"],
#   "client_id":"<application client id>",
#   "client_secret":"<application client secret>",
#   "note": "For Julia Metadata Listings"
# }
# ```
# Registered applications and authorised token are visible within your github profile page.
###################################################

require("JSON")
require("Calendar")
using Calendar

try
	global gh_auth
	gh_auth = ENV["GH_AUTH"]
catch e
	error ("Please provide a Github OAuth Token as environment variable GH_AUTH")
end

function gen_listpkg()

	Pkg.update()
	io=open("packages/packagelist.rst","w+");
	print(io, ".. _available-packages:\n\n")
	print(io, "********************\n Available Packages  \n********************\n\n")
	cd(Pkg.dir()) do
	for pkg in Pkg.available()
		print(" Processing $(pkg)\n")
		url = (Pkg.Read.url(pkg))
		maxv = Pkg.available(pkg)[end]
		url_reg = r"^(([^:/?#]+):)?(//([^/?#]*))?([^?#]*)(\?([^#]*))?(#(.*))?"
		gh_path_reg_git=r"^/(.*)?/(.*)?.git$"
		gh_path_reg_http=r"^/(.*)?/(.*)?$"
		m=match(url_reg, url)
		host=m.captures[4]
		path=m.captures[5]
		scheme=m.captures[2]
		u=get_user_details_default()

		if ismatch(r"github\.com", host)
			m2 = match(gh_path_reg_git, path)
			user=m2.captures[1]
			repo=m2.captures[2]
			u=get_user_details_gh(user)
			gh_repo_url = "https://api.github.com/repos/$(user)/$(repo)?access_token=$(gh_auth)"
			gh_contrib_url = "https://api.github.com/repos/$(user)/$(repo)/contributors?access_token=$(gh_auth)"
			travis_url= "https://api.travis-ci.org/repositories/$(user)/$(repo).json"
			#print("processing $gh_repo_url")
			gh_repo=JSON.parse(readall(download(gh_repo_url)))
			#print("processing $gh_user_url")
			gh_contrib=JSON.parse(readall(download(gh_contrib_url)))
			travis = Dict()
			try
				travis = JSON.parse(readall(download(travis_url)))
			catch
				print("Error reading travis url for $(repo)")
			end

			desc = get(gh_repo, "description", "No description provided")
			homepage = get(gh_repo, "homepage", nothing)
			html_url = gh_repo["html_url"]

			sha1_file = "METADATA/$pkg/versions/$(maxv)/sha1"
			if isfile(sha1_file)
				sha1 = readchomp(sha1_file)
				gh_update_url = "https://api.github.com/repos/$(user)/$(repo)/git/commits/$(sha1)?access_token=$(gh_auth)"
				gh_update = JSON.parse(readall(download(gh_update_url)))
				author = get(gh_update, "author", nothing)
				date = author != nothing ? get(author, "date", nothing) : nothing
			end
		end
		print(io, "`$(pkg) <$(html_url)>`_");
		if (get(travis, "file", nothing) != "not found" && get(travis, "last_build_status", nothing) != nothing)
			print(io, " |$(pkg)_build|\n")
			print(io, "_"^(length("`$(pkg) <$(html_url)>`_")) * "\n\n")
			print(io, "  .. |$(pkg)_build| image:: ../images/travis-icon.png\n")
			print(io, "     :height: 19\n")
			print(io, "     :width: 19\n")
			print(io, "     :target: https://travis-ci.org/$(user)/$(repo)\n\n")
		else
			print(io, "\n")
			print(io, "_"^(length("`$(pkg) <$(html_url)>`_")) * "\n\n")
		end

		print(io, "  .. image:: $(u[:avatar])\n     :height: 80px\n     :width: 80px\n     :align: right\n     :alt: $(u[:fullname])\n     :target: $(u[:url])\n\n")
		print(io, "  $(desc) \n\n")
		print(io, "  Current Version: ``$(maxv)``");
		if date != nothing
			print(io, "  (updated: $(date[1:10])) \n\n")
		end

		print(io, "  Maintainer: `$(u[:fullname]) <$(u[:url])>`_\n\n")

		if homepage != nothing && length(chomp(homepage)) > 0
			print(io, "  Documentation: `<$(homepage)>`_ \n\n")
		end
		print(io, "  Dependencies::\n\n" )
		ver_dir = "METADATA/$pkg/versions/$(maxv)/requires"
		if isfile(ver_dir)
			vset = Pkg.Reqs.parse(ver_dir)
			if length(vset) > 0
				for (dpkg, dver) in vset
					print(io, "      $(dpkg)"); print(io, " "^(25-length(dpkg))); print(io, "$(dver)\n")
				end
			else
				print(io, "      None\n")
			end
		else
			print(io, "      None\n")
		end
		print(io, "\n")

		if ismatch(r"github\.com", host)
			print(io, "  Contributors:\n\n")
			for contributor in gh_contrib
				c_user = get(contributor, "login", "")
				u=get_user_details_gh(c_user)
				print(io, "    .. image:: $(u[:avatar])\n        :height: 40px\n        :width: 40px\n")
				print(io, "        :alt: $(u[:fullname])\n        :target: $(u[:url])\n\n")
			end  #for contributor
		end

		print(io, "----\n\n")
	end  #for pkg
	print(io, ".. footer: $(length(Pkg.available())) packages, generated $(now()) \n\n")
	end  #cd

	close(io)
end #function

global user_cache = Dict{String, Dict}()

function get_user_details_default()
	u=Dict{Symbol, String}()
	u[:login] = "Unknown"
	u[:fullname] = "Unknown"
	u[:avatar] = "Not provided"
	u[:url] = "Not provided"

end

function get_user_details_gh(user)

	if !has(user_cache, user)
		gh_user_url = "https://api.github.com/users/$(user)?access_token=$(gh_auth)"
		gh_user=JSON.parse(readall(download(gh_user_url)))
		fullname = get(gh_user, "name", user)
		if fullname == nothing; fullname = user; end
		avatar =
		user_url = gh_user["html_url"]

	  u=Dict{Symbol, String}()

	  u[:login] = user
	  u[:avatar] = gh_user["avatar_url"]
	  #Sometimes name is missing, sometimes it is null in the JSON
	  if get(gh_user, "name", user) == nothing
	  	u[:fullname] = user
	  else
	  	u[:fullname] = get(gh_user, "name", user)
	  end
	  u[:url] = gh_user["html_url"]

	  user_cache[user] = u
	end

	return user_cache[user]

end #function

gen_listpkg()

max_stars_repo_name: porterjamesj/julia

max_issues_repo_name: porterjamesj/julia

max_forks_repo_name: porterjamesj/julia


text:
module Rice

export decode_aec, decode_rice, decode_rice_inv

function decode_aec(d, J = 16, n = 11; start_byte = 1)
    done = false
    i_byte = start_byte
    consumed_bits = 0
    samples = []
    while !done

        if consumed_bits >= 8
            i_byte += 1
            consumed_bits -= 8
        end

        if i_byte === length(d)
            b = (d[i_byte] << consumed_bits)
        else
            b = (d[i_byte] << consumed_bits) | (d[i_byte + 1] >> (8 - consumed_bits))
        end
        @info string(b, base = 2, pad = 8)

        if (b >> 3) === 0b00000
            # Zero Block
            @info "Zero Block" i_byte consumed_bits
            consumed_bits += 5

            counter = 0
            while true
                if consumed_bits >= 8
                    i_byte += 1
                    consumed_bits -= 8
                end

                if i_byte === length(d)
                    b = (d[i_byte] << consumed_bits)
                else
                    b = (d[i_byte] << consumed_bits) | (d[i_byte + 1] >> (8 - consumed_bits))
                end

                if counter > 63
                    error("impossible FS $counter")
                end

                if (b >> 7) === 0b1
                    consumed_bits += 1
                    break
                elseif (b >> 7) === 0b0
                    counter += 1
                    consumed_bits += 1
                else
                    error("WTF")
                end
            end

            if counter < 4
                counter += 1
            elseif counter === 4
                counter = 64
            elseif counter > 4
            else
                error("weird counter $counter")
            end

            @info repeat([0], counter * J) counter
            push!(samples, repeat([0], counter * J)...)

        elseif (b >> 3) === 0b00001
            # Second Extension
            @info "Second Extension" i_byte consumed_bits
            consumed_bits += 5
            error("not implemented")

        elseif (b >> 4) === 0b1111
            # No Compression
            @info "No Compression" i_byte consumed_bits
            consumed_bits += 4

            if consumed_bits >= 8
                i_byte += 1
                consumed_bits -= 8
            end

            if i_byte === length(d)
                b = (d[i_byte] << consumed_bits)
            else
                b = (d[i_byte] << consumed_bits) | (d[i_byte + 1] >> (8 - consumed_bits))
            end

            sample = Int(b) << (n - 8)
            consumed_bits += 8

            if consumed_bits >= 8
                i_byte += 1
                consumed_bits -= 8
            end

            if i_byte === length(d)
                b = (d[i_byte] << consumed_bits)
            else
                b = (d[i_byte] << consumed_bits) | (d[i_byte + 1] >> (8 - consumed_bits))
            end

            sample += Int(b) >> (16 - n)
            consumed_bits += 16 - n

            @info sample
            push!(samples, sample)

        else
            k = Int(b >> 4) - 1
            @info "Split Sample" i_byte consumed_bits k
            consumed_bits += 4

            for j in 1:J

                counter = 0
                while true
                    if consumed_bits >= 8
                        i_byte += 1
                        consumed_bits -= 8
                    end

                    if i_byte === length(d)
                        b = (d[i_byte] << consumed_bits)
                    else
                        b = (d[i_byte] << consumed_bits) | (d[i_byte + 1] >> (8 - consumed_bits))
                    end

                    if counter > 2^n - 1
                        error("impossible FS $counter")
                    end

                    if (b >> 7) === 0b1
                        consumed_bits += 1
                        break
                    elseif (b >> 7) === 0b0
                        counter += 1
                        consumed_bits += 1
                    else
                        error("WTF")
                    end
                end

                val = 0
                for i in 1:k
                    if consumed_bits >= 8
                        i_byte += 1
                        consumed_bits -= 8
                    end

                    if i_byte === length(d)
                        b = (d[i_byte] << consumed_bits)
                    else
                        b = (d[i_byte] << consumed_bits) | (d[i_byte + 1] >> (8 - consumed_bits))
                    end

                    val = (val << 1) + (b >> 7)
                    consumed_bits += 1
                end

                sample = counter * (2^k) + val
                @info sample k counter val
                push!(samples, sample)

                if sample > 2^n
                    error("impossible value $sample")
                end
            end

        end

        if i_byte === length(d)
            @info "end almost reached" i_byte consumed_bits
        end

    end
    samples
end

function decode_rice(d, k = 4, n = 11; start_byte = 1)
    done = false
    i_byte = start_byte
    consumed_bits = 0
    samples = []
    while !done
        counter = 0
        while true
            if consumed_bits >= 8
                i_byte += 1
                consumed_bits -= 8
            end

            if i_byte === length(d)
                b = (d[i_byte] << consumed_bits)
            elseif i_byte > length(d)
                @info "end overshot" i_byte consumed_bits string(d[i_byte-1], base = 2, pad = 8)
                return samples
            else
                b = (d[i_byte] << consumed_bits) | (d[i_byte + 1] >> (8 - consumed_bits))
            end

            if counter > 2^n - 1
                error("impossible FS $counter")
            end

            if (b >> 7) === 0b1
                consumed_bits += 1
                break
            elseif (b >> 7) === 0b0
                counter += 1
                consumed_bits += 1
            else
                error("WTF")
            end
        end

        val = 0
        for i in 1:k
            if consumed_bits >= 8
                i_byte += 1
                consumed_bits -= 8
            end

            if i_byte === length(d)
                b = (d[i_byte] << consumed_bits)
            elseif i_byte > length(d)
                @info "end overshot" i_byte consumed_bits string(d[i_byte-1], base = 2, pad = 8)
                return samples
            else
                b = (d[i_byte] << consumed_bits) | (d[i_byte + 1] >> (8 - consumed_bits))
            end

            val = (val << 1) + (b >> 7)
            consumed_bits += 1
        end

        sample = counter * (2^k) + val
        @info sample k counter val
        push!(samples, sample)

        if sample > 2^n
            error("impossible value $sample")
        end

        if i_byte === length(d)
            @info "end almost reached" i_byte consumed_bits string(d[i_byte], base = 2, pad = 8)
            break
        end

    end

    @info "end reached" i_byte consumed_bits string(d[i_byte], base = 2, pad = 8)
    return samples
end

function decode_rice_inv(d, k = 4, n = 11; start_byte = 1)
    done = false
    i_byte = start_byte
    consumed_bits = 0
    samples = []
    while !done
        while i_byte <= length(d)
            counter = 0
            while true
                if consumed_bits >= 8
                    i_byte += 1
                    consumed_bits -= 8
                end

                if i_byte === length(d)
                    b = (d[i_byte] << consumed_bits)
                else
                    b = (d[i_byte] << consumed_bits) | (d[i_byte + 1] >> (8 - consumed_bits))
                end

                if counter > 2^n - 1
                    error("impossible FS $counter")
                end

                if (b >> 7) === 0b0
                    consumed_bits += 1
                    break
                elseif (b >> 7) === 0b1
                    counter += 1
                    consumed_bits += 1
                else
                    error("WTF")
                end
            end

            val = 0
            for i in 1:k
                if consumed_bits >= 8
                    i_byte += 1
                    consumed_bits -= 8
                end

                if i_byte === length(d)
                    b = (d[i_byte] << consumed_bits)
                else
                    b = (d[i_byte] << consumed_bits) | (d[i_byte + 1] >> (8 - consumed_bits))
                end

                val = (val << 1) + (b >> 7)
                consumed_bits += 1
            end

            sample = counter * (2^k) + val
            @info sample k counter val
            push!(samples, sample)

            if sample > 2^n
                error("impossible value $sample")
            end
        end

        if i_byte === length(d)
            @info "end almost reached" i_byte consumed_bits
        end

    end
    samples
end

end  # module Rice


text:
export SOSDecomposition, SOSDecompositionWithDomain, sos_decomposition

"""
    struct SOSDecomposition{T, PT}

Represend SOSDecomposition without domain
"""
struct SOSDecomposition{T, PT <: MP.APL{T}, U} <: MP.AbstractPolynomialLike{U}
    ps::Vector{PT}
    function SOSDecomposition{T, PT, U}(ps::Vector{PT}) where {T, PT, U}
        new(ps)
    end
end

SOSDecomposition(ps::Vector{PT}) where {T, PT <: MP.APL{T}} = SOSDecomposition{T, PT, _promote_add_mul(T)}(ps)
MP.polynomialtype(::Union{SOSDecomposition{T, PT, U}, Type{SOSDecomposition{T, PT, U}}}) where {T, PT, U} = MP.polynomialtype(PT, U)

#function SOSDecomposition(ps::Vector)
#    T = reduce(promote_type, Int, map(eltype, ps))
#    SOSDecomposition{T}(ps)
#end

function GramMatrix(p::SOSDecomposition{T}) where {T}
    X = MP.mergemonovec(map(MP.monomials, p))
    m = length(p)
    n = length(X)
    Q = zeros(T, m, n)
    for i in 1:m
        j = 1
        for t in MP.terms(p[i])
            while X[j] != MP.monomial(t)
                j += 1
            end
            Q[i, j] = MP.coefficient(t)
            j += 1
        end
    end
    GramMatrix(Q' * Q, X)
end

function SOSDecomposition(p::GramMatrix, ranktol=0.0,
                          dec::MultivariateMoments.LowRankChol=SVDChol())
    n = length(p.basis)
    # TODO LDL^T factorization for SDP is missing in Julia
    # it would be nice to have though
    nM, cM, Q = MultivariateMoments.lowrankchol(Matrix(getmat(p)), dec, ranktol)
    ps = [MP.polynomial(Q[i,:], p.basis) for i in 1:size(Q, 1)]
    return SOSDecomposition(ps)
end
# Without LDL^T, we need to do float(T)
#SOSDecomposition(p::GramMatrix{C, T}) where {C, T} = SOSDecomposition{C, float(T)}(p)

Base.length(p::SOSDecomposition) = length(p.ps)
Base.isempty(p::SOSDecomposition) = isempty(p.ps)
Base.iterate(p::SOSDecomposition, args...) = Base.iterate(p.ps, args...)
Base.getindex(p::SOSDecomposition, i::Int) = p.ps[i]

(p::GramMatrix)(s::MP.AbstractSubstitution...) = MP.polynomial(p)(s...)

function Base.show(io::IO, p::SOSDecomposition)
    for (i, q) in enumerate(p)
        print(io, "(")
        print(io, q)
        print(io, ")^2")
        if i != length(p)
            print(io, " + ")
        end
    end
end

function Base.isapprox(p::SOSDecomposition, q::SOSDecomposition; kwargs...)
    m = length(p.ps)
    if length(q.ps) != m
        false
    else
        MultivariateMoments.permcomp((i, j) -> isapprox(p.ps[i], q.ps[j]; kwargs...), m)
    end
end

function Base.promote_rule(::Type{SOSDecomposition{T1, PT1, U1}}, ::Type{SOSDecomposition{T2, PT2, U2}}) where {T1, T2, PT1<:MP.APL{T1}, PT2<:MP.APL{T2}, U1, U2}
    T = promote_type(T1, T2)
    return SOSDecomposition{T, promote_type(PT1, PT2), _promote_add_mul(T)}
end

function Base.convert(::Type{SOSDecomposition{T, PT, U}}, p::SOSDecomposition) where {T, PT, U}
    return SOSDecomposition(convert(Vector{PT}, p.ps))
end

function MP.polynomial(decomp::SOSDecomposition)
    return sum(decomp.ps.^2)
end
function MP.polynomial(decomp::SOSDecomposition, T::Type)
    return MP.polynomial(MP.polynomial(decomp), T)
end

"""
    function sos_decomposition(cref::JuMP.ConstraintRef, ranktol::Float64, dec::MultivariateMoments.lowrankchol)

Return representation as a sum of squares.
"""
function sos_decomposition(cref::JuMP.ConstraintRef, args...)
    return SOSDecomposition(gram_matrix(cref), args...)
end

"""
    struct SOSDecompositionWithDomain{T, PT, S}

Represend SOSDecomposition on a basic semi-algebraic domain.
"""
struct SOSDecompositionWithDomain{T, PT <: MP.APL{T}, U, S <: AbstractSemialgebraicSet}
    sos::SOSDecomposition{T, PT, U}
    sosj::Vector{SOSDecomposition{T, PT, U}}
    domain::S
end

function SOSDecompositionWithDomain(ps::SOSDecomposition{T1, PT1, U1}, vps::Vector{SOSDecomposition{T2, PT2, U2}}, set::AbstractSemialgebraicSet ) where {T1, T2, PT1, PT2, U1, U2}
    ptype = promote_type(SOSDecomposition{T1, PT1, U1}, SOSDecomposition{T2, PT2, U2})
    return SOSDecompositionWithDomain(convert(ptype, ps), convert(Vector{ptype}, vps), set)
end

function Base.show(io::IO, decomp::SOSDecompositionWithDomain)
    print(io, decomp.sos)
    for (sos, g) in zip(decomp.sosj, inequalities(decomp.domain))
        print(io, " + ")
        print(io, sos)
        print(io, " * ")
        print(io, "(")
        print(io, g)
        print(io, ")")
    end
end

function MP.polynomial(decomp::SOSDecompositionWithDomain)
    p = MP.polynomial(decomp.sos)
    if !(isempty(equalities(decomp.domain)))
        @error "Semialgebraic set has equality constraints"
    end
    for (Gj, gj) in zip(decomp.sosj, inequalities(decomp.domain))
        p += MP.polynomial(Gj)*gj
    end
    return p
end

function Base.isapprox(p::SOSDecompositionWithDomain, q::SOSDecompositionWithDomain; kwargs...)
    return isapprox(p.sos, q.sos) && all(isapprox.(p.sosj, q.sosj)) && p.domain == q.domain
end

"""
    sos_decomposition(cref::JuMP.ConstraintRef, K<:AbstractBasicSemialgebraicSet)

Return representation in the quadraic module associated with K.
"""
function sos_decomposition(cref::JuMP.ConstraintRef, K::AbstractBasicSemialgebraicSet, args...)
    lm = SOSDecomposition.(lagrangian_multipliers(cref), args...)
    gm = sos_decomposition(cref, args...)
    return SOSDecompositionWithDomain(gm, lm, K)
end


text:
using LoopVectorization
using ExplainableAI
using ExplainableAI: modify_params
import ExplainableAI: modify_layer, lrp!
using Flux
using LinearAlgebra: I
using ReferenceTests
using Random

const RULES = Dict(
    "ZeroRule" => ZeroRule(),
    "EpsilonRule" => EpsilonRule(),
    "GammaRule" => GammaRule(),
    "ZBoxRule" => ZBoxRule(0.0f0, 1.0f0),
)

## Hand-written tests
@testset "ZeroRule analytic" begin
    rule = ZeroRule()

    ## Simple dense layer
    Rₖ₊₁ = reshape([1 / 3 2 / 3], 2, 1)
    aₖ = reshape([1.0 2.0;], 2, 1)
    W = [3.0 4.0; 5.0 6.0]
    b = [7.0, 8.0]
    Rₖ = reshape([17 / 90, 316 / 675], 2, 1) # expected output

    layer = Dense(W, b, relu)
    R̂ₖ = similar(aₖ) # will be inplace updated
    @inferred lrp!(R̂ₖ, rule, layer, aₖ, Rₖ₊₁)
    @test R̂ₖ ≈ Rₖ

    ## Pooling layer
    Rₖ₊₁ = Float32.([1 2; 3 4]//30)
    aₖ = Float32.([1 2 3; 10 5 6; 7 8 9])
    Rₖ = Float32.([0 0 0; 4 0 2; 0 0 4]//30) # expected output

    # Repeat in color channel dim and add batch dim
    Rₖ₊₁ = reshape(repeat(Rₖ₊₁, 1, 3), 2, 2, 3, 1)
    aₖ = reshape(repeat(aₖ, 1, 3), 3, 3, 3, 1)
    Rₖ = reshape(repeat(Rₖ, 1, 3), 3, 3, 3, 1)

    layer = MaxPool((2, 2); stride=(1, 1))
    R̂ₖ = similar(aₖ) # will be inplace updated
    @inferred lrp!(R̂ₖ, rule, layer, aₖ, Rₖ₊₁)
    @test R̂ₖ ≈ Rₖ
end

# Fixed pseudo-random numbers
T = Float32
pseudorandn(dims...) = randn(MersenneTwister(123), T, dims...)

## Test individual rules
@testset "modify_params" begin
    W, b = [1.0 -1.0; 2.0 0.0], [-1.0, 1.0]
    ρW, ρb = @inferred modify_params(GammaRule(0.42), W, b)
    @test ρW ≈ [1.42 -1.0; 2.84 0.0]
    @test ρb ≈ [-1.0, 1.42]
end

## Test Dense layer
# Define Dense test input
ins_dense = 20 # input dimension
outs_dense = 10 # output dimension
batchsize = 2
aₖ_dense = pseudorandn(ins_dense, batchsize)

layers = Dict(
    "Dense_relu" => Dense(ins_dense, outs_dense, relu; init=pseudorandn),
    "Dense_identity" => Dense(Matrix{Float32}(I, outs_dense, ins_dense), false, identity),
)
@testset "Dense" begin
    for (rulename, rule) in RULES
        @testset "$rulename" begin
            for (layername, layer) in layers
                @testset "$layername" begin
                    Rₖ₊₁ = layer(aₖ_dense)
                    Rₖ = similar(aₖ_dense)
                    @inferred lrp!(Rₖ, rule, layer, aₖ_dense, Rₖ₊₁)

                    @test typeof(Rₖ) == typeof(aₖ_dense)
                    @test size(Rₖ) == size(aₖ_dense)

                    if rulename == "Dense_identity"
                        # First `outs_dense` dimensions should propagate
                        # activations as relevances, rest should be ≈ 0.
                        @test Rₖ[1:outs_dense] ≈ aₖ_dense[1:outs_dense]
                        @test all(Rₖ[outs_dense:end] .< 1e-8)
                    end

                    @test_reference "references/rules/$rulename/$layername.jld2" Dict(
                        "R" => Rₖ
                    ) by = (r, a) -> isapprox(r["R"], a["R"]; rtol=0.02)
                end
            end
        end
    end
end

## Test PoolingLayers
insize = (6, 6, 2, batchsize)
aₖ = pseudorandn(insize)

equalpairs = Dict( # these pairs of layers are all equal
    "AdaptiveMaxPool" => (AdaptiveMaxPool((3, 3)), MaxPool((2, 2); pad=0)),
    "GlobalMaxPool" => (GlobalMaxPool(), MaxPool((6, 6); pad=0)),
    "AdaptiveMeanPool" => (AdaptiveMeanPool((3, 3)), MeanPool((2, 2); pad=0)),
    "GlobalMeanPool" => (GlobalMeanPool(), MeanPool((6, 6); pad=0)),
)

@testset "PoolingLayers" begin
    for (rulename, rule) in RULES
        @testset "$rulename" begin
            for (layername, layers) in equalpairs
                @testset "$layername" begin
                    l1, l2 = layers
                    Rₖ₊₁ = l1(aₖ)
                    @test Rₖ₊₁ == l2(aₖ)
                    Rₖ1 = similar(aₖ)
                    Rₖ2 = similar(aₖ)
                    @inferred lrp!(Rₖ1, rule, l1, aₖ, Rₖ₊₁)
                    @inferred lrp!(Rₖ2, rule, l2, aₖ, Rₖ₊₁)
                    @test Rₖ1 == Rₖ2

                    @test typeof(Rₖ1) == typeof(aₖ)
                    @test size(Rₖ1) == size(aₖ)

                    @test_reference "references/rules/$rulename/$layername.jld2" Dict(
                        "R" => Rₖ1
                    ) by = (r, a) -> isapprox(r["R"], a["R"]; rtol=0.02)
                end
            end
        end
    end
end

## Test ConvLayers and others
layers = Dict(
    "Conv" => Conv((3, 3), 2 => 4; init=pseudorandn),
    "MaxPool" => MaxPool((3, 3)),
    "MeanPool" => MaxPool((3, 3)),
    "ConvTranspose" => ConvTranspose((3, 3), 2 => 4; init=pseudorandn),
    "CrossCor" => CrossCor((3, 3), 2 => 4; init=pseudorandn),
    "flatten" => flatten,
    "Dropout" => Dropout(0.2),
    "AlphaDropout" => AlphaDropout(0.2),
)
@testset "Other Layers" begin
    for (rulename, rule) in RULES
        @testset "$rulename" begin
            for (layername, layer) in layers
                @testset "$layername" begin
                    Rₖ₊₁ = layer(aₖ)
                    Rₖ = similar(aₖ)
                    @inferred lrp!(Rₖ, rule, layer, aₖ, Rₖ₊₁)

                    @test typeof(Rₖ) == typeof(aₖ)
                    @test size(Rₖ) == size(aₖ)

                    @test_reference "references/rules/$rulename/$layername.jld2" Dict(
                        "R" => Rₖ
                    ) by = (r, a) -> isapprox(r["R"], a["R"]; rtol=0.02)
                end
            end
        end
    end
end

## Test custom layers & default AD fallback using the ZeroRule
# Compare with references of non-wrapped layers
struct TestWrapper{T}
    layer::T
end
(w::TestWrapper)(x) = w.layer(x)
modify_layer(r::AbstractLRPRule, w::TestWrapper) = modify_layer(r, w.layer)
lrp!(Rₖ, rule::ZBoxRule, w::TestWrapper, aₖ, Rₖ₊₁) = lrp!(Rₖ, rule, w.layer, aₖ, Rₖ₊₁)

layers = Dict(
    "Conv" => (Conv((3, 3), 2 => 4; init=pseudorandn), aₖ),
    "Dense_relu" => (Dense(ins_dense, outs_dense, relu; init=pseudorandn), aₖ_dense),
    "flatten" => (flatten, aₖ),
)
@testset "Custom layers" begin
    for (rulename, rule) in RULES
        @testset "$rulename" begin
            for (layername, (layer, aₖ)) in layers
                @testset "$layername" begin
                    wrapped_layer = TestWrapper(layer)
                    Rₖ₊₁ = wrapped_layer(aₖ)
                    Rₖ = similar(aₖ)
                    @inferred lrp!(Rₖ, rule, wrapped_layer, aₖ, Rₖ₊₁)

                    @test typeof(Rₖ) == typeof(aₖ)
                    @test size(Rₖ) == size(aₖ)

                    @test_reference "references/rules/$rulename/$layername.jld2" Dict(
                        "R" => Rₖ
                    ) by = (r, a) -> isapprox(r["R"], a["R"]; rtol=0.02)
                end
            end
        end
    end
end


text:
using Test
using P4est

const P4EST_TEST = get(ENV, "P4EST_TEST", "")

@testset "P4est.jl tests" begin
  @testset "p4est_connectivity_new_periodic" begin
    @test p4est_connectivity_new_periodic() isa Ptr{p4est_connectivity}
  end

  connectivity = p4est_connectivity_new_periodic()
  @testset "p4est_connectivity_is_valid" begin
    @test p4est_connectivity_is_valid(connectivity) == 1
  end

  @testset "unsafe_wrap" begin
    obj = unsafe_wrap(connectivity)
    @test obj.num_vertices == 4
  end

  @testset "uses_mpi" begin
    if P4EST_TEST == "P4EST_JLL_USES_MPI_PRE_GENERATED_BINDINGS"
      @test P4est.uses_mpi() == true
    elseif P4EST_TEST == "P4EST_JLL_USES_MPI"
      @test P4est.uses_mpi() == true
    elseif P4EST_TEST == "P4EST_CUSTOM_NON_MPI"
      @test P4est.uses_mpi() == false
    elseif P4EST_TEST == "P4EST_CUSTOM_USES_MPI"
      @test P4est.uses_mpi() == true
    end
  end
end


text:
# This will be called by __init__
function init_us_census()
    register(DataDep(
        "US Census - names",
        """US Census data, 1880-2017.
            https://catalog.data.gov/dataset/baby-names-from-social-security-card-applications-national-level-data""",
        "https://www.ssa.gov/oact/babynames/names.zip",
        post_fetch_method=unpack
        ))
end


text:
int main() {
        x = 14;
	return 0 ;
}

max_stars_repo_name: carlostome/Javalette-Compiler

max_issues_repo_name: carlostome/Javalette-Compiler

max_forks_repo_name: carlostome/Javalette-Compiler


text:
using EffectiveWaves, Test

@testset "test different numerical integration methods" begin

    using ApproxFun
    import SpecialFunctions: hankelh1

    # function imitates kernal in Fredholm equation that is solved for AverageWave
    n = 1; X = 0.9; θin = 0.2;
    num_coefs = 10000;
    tol = (1.0/num_coefs)^2

    K(Y) = cos(Y*sin(θin) + n*atan(Y,X))*hankelh1(n,sqrt(X^2+Y^2))

    # approximate function with Chebyshev polynomial (to high precision) then integrate from Y0 to Y1
    Y0 = 0.0; Y1 = 1000.0;
    s = sum(Fun(K,Y0..Y1, num_coefs))

    # compare with basic discrete integration
    Ys = LinRange(Y0,Y1,2*num_coefs+1);
    σ = integration_scheme(Ys);
    @test abs(sum(K.(Ys).*σ) - s)/abs(s) < 500*tol # trapezoidal scheme not so accurate

    σ = integration_scheme(Ys; scheme=:simpson);
    @test abs(sum(K.(Ys).*σ) - s)/abs(s) < tol

    integration_scheme(Ys; scheme=:fish);
    @test true
end


text:
"""
    CoboundaryMatrix{I}

This `struct` is used to compute cohomology. The `I` parameter sets whether the implicit
algoritm is used or not.
"""
struct CoboundaryMatrix{
    I,T<:Number,F,S<:AbstractCell,R<:ReducedMatrix,B<:Chain{T},C<:Chain{T}
}
    filtration::F
    reduced::R
    buffer::B
    chain::C
    columns_to_reduce::Vector{S}
    columns_to_skip::Vector{S}
end

function CoboundaryMatrix{I}(
    ::Type{T}, filtration, columns_to_reduce, columns_to_skip
) where {I,T}
    S = eltype(columns_to_reduce)
    C = simplex_type(filtration, dim(S) + 1)
    if I
        reduced = ReducedMatrix{C,T,S}()
        buffer = Chain{T,S}()
    else
        reduced = ReducedMatrix{C,T,C}()
        buffer = Chain{T,C}()
    end
    sizehint!(reduced, length(columns_to_reduce))
    chain = Chain{T,C}()

    return CoboundaryMatrix{
        I,T,typeof(filtration),S,typeof(reduced),typeof(buffer),typeof(chain)
    }(
        filtration, reduced, buffer, chain, columns_to_reduce, columns_to_skip
    )
end

field_type(::CoboundaryMatrix{<:Any,T}) where {T} = T
dim(cm::CoboundaryMatrix{<:Any,<:Any,<:Any,S}) where {S} = dim(S)

is_implicit(::CoboundaryMatrix{I}) where {I} = I
is_cohomology(::CoboundaryMatrix) = true
ordering(::CoboundaryMatrix) = Base.Order.Forward

function coboundary(matrix::CoboundaryMatrix, simplex::AbstractCell)
    return coboundary(matrix.filtration, simplex)
end

function next_matrix(matrix::CoboundaryMatrix{I}, verbose) where {I}
    C = simplex_type(matrix.filtration, dim(matrix) + 1)
    new_to_reduce = C[]
    new_to_skip = C[]
    sizehint!(new_to_skip, length(matrix.reduced))

    if verbose
        progbar = ProgressUnknown("Assembling columns:")
    end
    for simplex in columns_to_reduce(
        matrix.filtration,
        Iterators.flatten((matrix.columns_to_reduce, matrix.columns_to_skip)),
    )
        if haskey(matrix.reduced, simplex)
            push!(new_to_skip, abs(simplex))
        else
            push!(new_to_reduce, abs(simplex))
        end
        verbose && next!(
            progbar;
            showvalues=(
                ("cleared", length(new_to_skip)), ("to reduce", length(new_to_reduce))
            ),
        )
    end
    @prog_print verbose '\r'

    return CoboundaryMatrix{I}(
        field_type(matrix), matrix.filtration, new_to_reduce, new_to_skip
    )
end

"""
    BoundaryMatrix{I}

This `struct` is used to compute homology. The `I` parameter sets whether the implicit
algoritm is used or not.
"""
struct BoundaryMatrix{
    I,T<:Number,F,S<:AbstractCell,R<:ReducedMatrix,B<:Chain{T},C<:Chain{T}
}
    filtration::F
    reduced::R
    buffer::B
    chain::C
    columns_to_reduce::Vector{S}
end

function BoundaryMatrix{I}(::Type{T}, filtration, columns_to_reduce) where {I,T}
    if eltype(columns_to_reduce) === Any
        S = typeof(first(columns_to_reduce))
    else
        S = eltype(columns_to_reduce)
    end
    C = simplex_type(filtration, dim(S) - 1)

    columns = S[]
    foreach(columns_to_reduce) do c
        push!(columns, abs(c))
    end

    if !I
        reduced = ReducedMatrix{C,T,C}()
        buffer = Chain{T,C}()
    else
        reduced = ReducedMatrix{C,T,S}()
        buffer = Chain{T,S}()
    end
    sizehint!(reduced, length(columns))
    chain = Chain{T,C}()

    return BoundaryMatrix{
        I,T,typeof(filtration),S,typeof(reduced),typeof(buffer),typeof(chain)
    }(
        filtration, reduced, buffer, chain, columns
    )
end

field_type(::BoundaryMatrix{<:Any,T}) where {T} = T
dim(::BoundaryMatrix{<:Any,<:Any,<:Any,S}) where {S} = dim(S) - 1
ordering(::BoundaryMatrix) = Base.Order.Reverse

is_implicit(::BoundaryMatrix{I}) where {I} = I
is_cohomology(::BoundaryMatrix) = false

# The naming here is not ideal...
function coboundary(matrix::BoundaryMatrix, simplex::AbstractCell)
    return boundary(matrix.filtration, simplex)
end

"""
    initialize_coboundary!(matrix, column)

Initialize the column indexed by `column` by adding its (co)boundary to `matrix.chain`. This
is where the emergent pairs optimization gets triggered for implicit versions of the
algorithm.

"""
function initialize_coboundary!(matrix, column)
    empty!(matrix.chain)
    # Emergent pairs: we are looking for pairs of simplices (σ, τ) where σ is the youngest
    # facet of τ and τ is the oldest cofacet of σ. These pairs give birth to persistence
    # intervals with zero length and can be skipped.

    # This implementation of this optimization only works if (co)boundary simplices are
    # returned in the correct order and if the birth times of σ and τ are the same.
    emergent_check = emergent_pairs(matrix.filtration) && is_implicit(matrix)
    for cofacet in coboundary(matrix, column)
        if emergent_check && birth(cofacet) == birth(column)
            emergent_check = false
            if !haskey(matrix.reduced, cofacet)
                return ChainElement{typeof(cofacet),field_type(matrix)}(cofacet)
            end
        end
        push!(matrix.chain, cofacet)
    end
    if isempty(matrix.chain)
        return nothing
    else
        heapify!(matrix.chain, ordering(matrix))
        return heappop!(matrix.chain, ordering(matrix))
    end
end

"""
    add!(matrix, column, pivot)

Add already reduced column indexed by `column` multiplied by `-coefficient(pivot)` to
`matrix.chain`.

"""
function add!(matrix, column, pivot)
    return add!(Val(is_implicit(matrix)), matrix, column, pivot)
end
# Implicit version
function add!(::Val{true}, matrix, column, pivot)
    factor = -coefficient(pivot)
    for element in column
        for cofacet in coboundary(matrix, simplex(element))
            simplex(pivot) == cofacet && continue
            heappush!(
                matrix.chain, (cofacet, coefficient(element) * factor), ordering(matrix)
            )
        end
    end
    return append!(matrix.buffer, (s, c * factor) for (s, c) in column)
end
# Explicit version
function add!(::Val{false}, matrix, column, pivot)
    factor = -coefficient(pivot)
    for element in column
        # The pivot is not stored in the column, so there is no need to check for it.
        heappush!(matrix.chain, element * factor, ordering(matrix))
    end
end

"""
    finalize!(matrix, column, pivot)

After reduction is done, finalize the column by adding it to the reduced matrix.
"""
function finalize!(matrix, column, pivot)
    return finalize!(Val(is_implicit(matrix)), matrix, column, pivot)
end
# Implicit version
function finalize!(::Val{true}, matrix, column, pivot)
    push!(matrix.buffer, column)
    return matrix.reduced[simplex(pivot)] = clean!(
        matrix.buffer, ordering(matrix), inv(coefficient(pivot))
    )
end
# Explicit version
function finalize!(::Val{false}, matrix, column, pivot)
    return matrix.reduced[simplex(pivot)] = clean!(
        matrix.chain, ordering(matrix), inv(coefficient(pivot))
    )
end

"""
    reduce_column!(matrix, column_to_reduce)

Reduce `column_to_reduce` by repeatedly adding other columns to it. Once nothing more can be
added, `finalize!` the column.

"""
function reduce_column!(matrix, column_to_reduce)
    empty!(matrix.buffer)
    pivot = initialize_coboundary!(matrix, column_to_reduce)
    while !isnothing(pivot)
        column = matrix.reduced[pivot]
        isempty(column) && break

        add!(matrix, column, pivot)
        pivot = heappop!(matrix.chain, ordering(matrix))
    end
    if !isnothing(pivot)
        finalize!(matrix, column_to_reduce, pivot)
    end

    return pivot
end

"""
    collect_cocycle!(matrix, pivot, reps)

Collect the representative (co)cycle.
"""
function collect_cocycle!(matrix, column, pivot)
    if is_cohomology(matrix)
        if !is_implicit(matrix)
            error("representative cocycles for explicit cohomology not supported")
        elseif isnothing(pivot)
            push!(matrix.buffer, column)
            return copy(clean!(matrix.buffer, ordering(matrix)))
        else
            return Chain(matrix.reduced[pivot])
        end
    else
        if is_implicit(matrix)
            heappush!(matrix.chain, pivot, ordering(matrix))
            return heapmove!(matrix.chain, ordering(matrix))
        else
            rep = Chain(matrix.reduced[pivot])
            push!(rep, pivot)
            return clean!(rep, ordering(matrix))
        end
    end
end

"""
    interval(matrix, column, pivot, cutoff, reps)

Construct a persistence interval.
"""
function interval(matrix, column, pivot, cutoff, reps)
    if is_cohomology(matrix)
        birth_simplex = column
        death_simplex = isnothing(pivot) ? nothing : simplex(pivot)
    elseif isnothing(pivot)
        # In homology, birth simplex is nothing when column is fully reduced.
        return nothing
    else
        birth_simplex, death_simplex = simplex(pivot), column
    end
    birth_time = Float64(birth(birth_simplex))
    death_time = isnothing(death_simplex) ? Inf : Float64(birth(death_simplex))
    if death_time - birth_time > cutoff
        if reps
            rep = (; representative=collect_cocycle!(matrix, column, pivot))
        else
            rep = NamedTuple()
        end
        meta = (; birth_simplex=birth_simplex, death_simplex=death_simplex, rep...)
        return PersistenceInterval(birth_time, death_time, meta)
    else
        return nothing
    end
end

"""
    handle_apparent_pairs!(matrix, intervals, cutoff, verbose, reps)

Handle apparent pair precomputation, if defined for `matrix.filtration`. Only does anything
for implicit cohomology. Resulting intervals (if any) are `push!`ed to `intervals`.
"""
function handle_apparent_pairs!(matrix, intervals, cutoff, verbose, reps)
    coho = Val(is_cohomology(matrix))
    impl = Val(is_implicit(matrix))
    return handle_apparent_pairs!(coho, impl, matrix, intervals, cutoff, verbose, reps)
end
# Implicit cohomology version
function handle_apparent_pairs!(
    ::Val{true}, ::Val{true}, matrix, intervals, cutoff, verbose, reps
)
    columns, apparent = find_apparent_pairs(
        matrix.filtration, matrix.columns_to_reduce, verbose
    )
    bulk_insert!(matrix.reduced, apparent)
    for (σ, τ) in apparent
        τ_elem = ChainElement{typeof(τ),field_type(matrix)}(τ)
        int = interval(matrix, σ, τ_elem, cutoff, reps)
        !isnothing(int) && push!(intervals, int)
    end
    return columns
end
# Other versions don't support this
function handle_apparent_pairs!(::Val, ::Val, matrix, _, _, _, _)
    return matrix.columns_to_reduce
end

"""
    compute_intervals!(matrix, cutoff, verbose, reps)

Compute all intervals by fully reducing `matrix`.
"""
function compute_intervals!(matrix, cutoff, verbose, reps)
    ###
    ### Set up output.
    ###
    intervals = PersistenceInterval[]

    ###
    ### Apparent pair stuff.
    ###
    columns = handle_apparent_pairs!(matrix, intervals, cutoff, verbose, reps)

    ###
    ### Interval computation.
    ###
    @prog_print(
        verbose,
        fmt_number(length(columns)),
        " ",
        simplex_name(eltype(columns)),
        " to reduce.",
    )
    # One-dimensional columns in cohomology are already sorted.
    if !is_cohomology(matrix) || dim(matrix) > 1
        sort_t = time_ns()
        sort!(columns; rev=is_cohomology(matrix))
        elapsed = round((time_ns() - sort_t) / 1e9; digits=3)
        @prog_println verbose " Sorted in " elapsed "s."
    else
        @prog_println verbose
    end

    if verbose
        progbar = Progress(length(columns); desc="Computing $(dim(matrix))d intervals... ")
    end
    for column in columns
        pivot = reduce_column!(matrix, column)
        int = interval(matrix, column, pivot, cutoff, reps)
        !isnothing(int) && push!(intervals, int)

        verbose && next!(progbar; showvalues=((:intervals, length(intervals)),))
    end

    diagram = PersistenceDiagram(
        sort!(intervals; by=persistence);
        threshold=Float64(threshold(matrix.filtration)),
        dim=dim(matrix),
        field=field_type(matrix),
        filtration=matrix.filtration,
    )
    return postprocess_diagram(matrix.filtration, diagram)
end

"""
    compute_death_simplices!(matrix::CoboundaryMatrix{true}, verbose, cutoff)

Fully reduce `matrix`, but only compute (homological) death simplices. Return all death
simplices up to the last that produces an interval with persistence greater than `cutoff`.

Used for assisted homology.
"""
function compute_death_simplices!(matrix::CoboundaryMatrix{true}, verbose, cutoff)
    columns, apparent = find_apparent_pairs(
        matrix.filtration, matrix.columns_to_reduce, verbose
    )
    bulk_insert!(matrix.reduced, apparent)
    deaths = simplex_type(matrix.filtration, dim(matrix) + 1)[]
    inf_births = simplex_type(matrix.filtration, dim(matrix))[]
    if isempty(columns)
        return deaths, inf_births
    else
        dim(matrix) > 1 && sort!(columns; rev=true)
        thresh = typemin(birth(first(columns)))
        for pair in apparent
            if birth(pair[2]) - birth(pair[1]) > cutoff
                thresh = max(thresh, birth(pair[2]))
            end
            push!(deaths, pair[2])
        end
        if verbose
            progbar = Progress(length(columns); desc="Precomputing columns...   ")
        end
        for column in columns
            pivot = reduce_column!(matrix, column)
            if !isnothing(pivot)
                if birth(pivot) - birth(column) > cutoff
                    thresh = max(thresh, birth(pivot))
                end
                push!(deaths, simplex(pivot))
            else
                push!(inf_births, column)
            end
            verbose && next!(progbar; showvalues=((:simplices, length(deaths)),))
        end
        return filter!(x -> birth(x) ≤ thresh, deaths), inf_births
    end
end


text:
const PlotTypes{Typ, Arg} = Union{Atomic{Typ, Arg}, Combined{Typ, Arg}}

argtypes(x::PlotTypes{T, A}) where {T, A} = A
argtypes(x) = Any

function data_limits(x)
    error("No datalimits for $(typeof(x)) and $(argtypes(x))")
end
"""
Data limits calculate a minimal boundingbox from the data points in a plot.
This doesn't include any transformations, markers etc.
"""
data_limits(x::Atomic{Typ, <: Tuple{Arg1}}) where {Typ, Arg1} = FRect3D(value(x[1]))
# TODO makes this generically work
data_limits(x::Atomic{Typ, <: Tuple{<: AbstractVector{<: NTuple{N, <: Number}}}}) where {Typ, N} = FRect3D(Point{N, Float32}.(value(x[1])))
function data_limits(x::Atomic{Typ, <: Tuple{<: AbstractVector{<: NTuple{2, T}}}}) where {Typ, T <: VecTypes}
    FRect3D(reinterpret(T, value(x[1])))
end

function data_limits(x::Atomic{Typ, <: Tuple{X, Y, Z}}) where {Typ, X, Y, Z}
    _boundingbox(value.(x[1:3])...)
end

function data_limits(x::Atomic{Typ, <: Tuple{X, Y}}) where {Typ, X, Y}
    _boundingbox(value.(x[1:2])...)
end

_isfinite(x) = isfinite(x)
_isfinite(x::VecTypes) = all(isfinite, x)
scalarmax(x::AbstractArray, y::AbstractArray) = max.(x, y)
scalarmax(x, y) = max(x, y)
scalarmin(x::AbstractArray, y::AbstractArray) = min.(x, y)
scalarmin(x, y) = min(x, y)

extrema_nan(itr::Pair) = (itr[1], itr[2])
extrema_nan(itr::ClosedInterval) = (minimum(itr), maximum(itr))

function extrema_nan(itr)
    vs = iterate(itr)
    vs === nothing && return (NaN, NaN)
    v, s = vs
    vmin = vmax = v
    # find first finite value
    while vs !== nothing && !_isfinite(v)
        v, s = vs
        vmin = vmax = v
        vs = iterate(itr, s)
    end
    while vs !== nothing
        x, s = vs
        vs = iterate(itr, s)
        _isfinite(x) || continue
        vmax = scalarmax(x, vmax)
        vmin = scalarmin(x, vmin)
    end
    return (vmin, vmax)
end

function _boundingbox(x, y, z = (0 => 0))
    minmax = extrema_nan.((x, y, z))
    mini, maxi = first.(minmax), last.(minmax)
    FRect3D(mini, maxi .- mini)
end

const ImageLike{Arg} = Union{Heatmap{Arg}, Image{Arg}}
function data_limits(x::ImageLike{<: Tuple{X, Y, Z}}) where {X, Y, Z}
    _boundingbox(value.((x[1], x[2]))...)
end

function data_limits(x::Volume)
    _boundingbox(value.((x[1], x[2], x[3]))...)
end

function text_limits(x::VecTypes)
    p = to_ndim(Vec3f0, x, 0.0)
    FRect3D(p, p)
end
function text_limits(x::AbstractVector)
    FRect3D(x)
end
function data_limits(x::Text)
    text_limits(value(x[:position]))
end

function data_limits(x::Annotations)
    # data limits is supposed to not include any transformation.
    # for the annotation, we use the model matrix directly, so we need to
    # to inverse that transformation for the correct limits
    bb = data_limits(x.plots[1])
    inv(value(x[:model])) * bb
end

Base.isfinite(x::Rect) = all(isfinite.(minimum(x))) &&  all(isfinite.(maximum(x)))

function data_limits(plots::Vector)
    isempty(plots) && return FRect3D(Vec3f0(0), Vec3f0(0))
    bb = FRect3D()
    plot_idx = iterate(plots)
    while plot_idx !== nothing
        plot, idx = plot_idx
        plot_idx = iterate(plots, idx)
        # axis shouldn't be part of the data limit
        isaxis(plot) && continue
        isa(plot, Legend) && continue
        bb2 = data_limits(plot)
        isfinite(bb) || (bb = bb2)
        isfinite(bb2) || continue
        bb = union(bb, bb2)
    end
    bb
end

data_limits(s::Scene) = data_limits(plots_from_camera(s))
data_limits(plot::Combined) = data_limits(plot.plots)


text:
module test_9
using AdventOfCode2018.Day9
using Test

@testset "day 8, part 1" begin
    @test parse_input("10 players; last marble is worth 1618 points") == (10, 1618)
    @test part_1(9, 25) == 32
    @test part_1(10, 1618) == 8317
    @test part_1(13, 7999) == 146373
    @test part_1(17, 1104) == 2764
    @test part_1(21, 6111) == 54718
    @test part_1(30, 5807) == 37305
end

end


text:
using FEMfunctions

import Gmsh: gmsh
using LinearAlgebra

using AbstractPlotting
using SparseArrays
using SpecialMatrices
using Polynomials
using Printf
using WriteVTK

include("makeMesh.jl")
include("twodOseen.jl")

μ = 1.0
ω =-1.0

### define the N parameters that describe the inner boundary
N = 160      # number of BSplines used to represent the inner boundary

### describe the inner boundary using BSplines
r = ones(N,1)

x,eConn,eConn2, innerNodes,innerX, outerNodes,outerX = makeMesh(r)

sort!(innerNodes);
sort!(outerNodes);

#   Let's look at the mesh...
nNodes = size(x,2)
nElements = size(eConn,2)

xT = zeros(Float64,nNodes,2)   # xT = transpose(x[1:2,:])
eC = zeros(Int64,nElements,6)  # eC = transpose(eConn), converted to Int64
for i=1:nNodes
  xT[i,1] = x[1,i]
  xT[i,2] = x[2,i]
end
for i=1:nElements#-nElementsSensor # the sensor elements are orientated correctly (why?)
  eC[i,1] = convert(Int64,eConn[1,i])
  eC[i,2] = convert(Int64,eConn[3,i])
  eC[i,3] = convert(Int64,eConn[2,i])
  eC[i,4] = convert(Int64,eConn[6,i])
  eC[i,5] = convert(Int64,eConn[5,i])
  eC[i,6] = convert(Int64,eConn[4,i])
end

advection = Array{Float64,2}(undef,nNodes,2)
for i=1:nNodes
  advection[i,1] = -ω*xT[i,2]
  advection[i,2] = ω*xT[i,1]
end

velocity,pressure,exactVelocity,exactPressure = twodOseen(xT,eC,innerNodes,outerNodes,advection,μ,ω)

#  To manually check the function saveFEMasVTK, enter the following line
#  and use ParaView:
# saveFEMasVTK("Oseen",xT,eC,["pressure"],pressure,["velocity"],velocity)
# saveFEMasVTK("OseenEx",xT,eC,["pressure"],exactPressure,["velocity"],exactVelocity)

M = twodMassMatrix(xT,eC)
errorU = velocity[:,1]-exactVelocity[:,1]
errorV = velocity[:,2]-exactVelocity[:,2]
errorP = pressure-exactPressure

#  The inner boundary isn't exactly a circle of radius 1, so we have a weaker error
@test sqrt( dot(errorU,(M*errorU)) + dot(errorV,(M*errorV)) ) < 1e-4
#         sqrt( dot(errorP,(M*errorP)) ) < 1e-2


text:
# SPD-License-Identifier: MIT

using ArgParse
using BayesianModels
using CUDA
using HDF5
using JLD2
using LinearAlgebra
using Random

CUDA.allowscalar(false)

include("utils.jl")

function parse_commandline()
    s = ArgParseSettings()
    @add_arg_table s begin
        "model"
            required = true
            help = "input GMM"
        "features"
            required = true
            help = "HDF5 archive containing the features"
    end
    s.description = """
    Compute the GMM most likely components. The list of utterances
    to use is read from the standard input and the results is output
    on the standard output.
    """
    parse_args(s)
end

function main(args)
    model, _ = loadmodel(args["model"])
    h5open(args["features"], "r") do ffea
        for utt in eachline(stdin)
            X = read(ffea[utt])
            println("$utt\t$(join(predict(model, X), " "))")
        end
    end
end

args = parse_commandline()
main(args)

max_stars_repo_name: MartinKocour/SpeechLab

max_issues_repo_name: MartinKocour/SpeechLab

max_forks_repo_name: MartinKocour/SpeechLab


text:
using Revise
using NNFEM
using PyPlot

using Revise
using ADCME
using PyCall
using Test
using MAT
using NNFEM
using PyPlot
using JLD2
using LinearAlgebra
using Distributions, Random
using ForwardDiff
using DelimitedFiles

np = pyimport("numpy")
"""
Force load on one edge of the plate,
The edge is [0,L], which is discretized to ne elements, with porder-th order polynomials
type: a string of load type, which can be "constant", or "Gaussian"
args: an array,
      for Constant load, it has p1 and p2, in tangential and normal direction
      for "Gaussian" load, it has p, x0, and σ, the force is in the normal direction

"""
function ComputeLoad(L, ne, porder, ngp, type,  args)
    @assert(ngp <= 4)
    dx = L/ne

    xx = np.linspace(0, L, ne*porder + 1)

    Ft, Fn = zeros(Float64, ne*porder + 1), zeros(Float64, ne*porder + 1)

    elem_coords = zeros(Float64, porder + 1, 2)

    # construct pressure load function
    if type == "Constant"
        pt, pn = args
        f_px, f_py = x->pt, x->pn
    elseif type == "Gaussian"
        pn, x0, σ = args
        f_px, f_py = x->0, x-> pn * 1.0/(sqrt(2*pi*σ^2)) * exp.(-0.5*(x .- x0).^2/σ^2)
    else
        error("Force load type is not recogonized ", type)
    end

    # assemble the force

    @show ne
    for e = 1:ne
        if porder == 1
            loc_id = [e, e+1]
        elseif porder == 2
            loc_id = [2*(e-1)+1, 2*(e-1)+3, 2*(e-1)+2]
        else
            error("porder error porder ==", porder)
        end

        loc_xx = xx[loc_id]         #x-coordinates of element nodes
        elem_coords[:,1] .= loc_xx

        #return list
        weights, hs = get1DElemShapeData( elem_coords, ngp)

        for igp = 1:ngp
            gp_xx = loc_xx' * hs[igp]   #x-coordinates of element Gaussian points
            Ft[loc_id] += f_px(gp_xx) * hs[igp] * weights[igp]
            Fn[loc_id] += f_py(gp_xx) * hs[igp] * weights[igp]

        end
    end

    return Ft, Fn
end

function BoundaryCondition(tid, nx, ny, porder=2, Lx = 1.0, Ly = 0.5; force_scale=5.0)
    nnodes, neles = (nx*porder + 1)*(ny*porder + 1), nx*ny

    x = np.linspace(0.0, Lx, nx*porder + 1)
    y = np.linspace(0.0, Ly, ny*porder + 1)

    X, Y = np.meshgrid(x, y)
    nodes = zeros(nnodes,2)
    nodes[:,1], nodes[:,2] = X'[:], Y'[:]
    ndofs = 2

    EBC, g = zeros(Int64, nnodes, ndofs), zeros(nnodes, ndofs)
    FBC, fext = zeros(Int64, nnodes, ndofs), zeros(nnodes, ndofs)
    gt = nothing

    #todo only change 203

    # setting EBC:
    # tid = 1XX   |==|-->
    # tid = 2XX
    if div(tid,100)==1
        EBC[collect(1:nx*porder+1), :] .= -1 # fix bottom
        FBC[collect((nx*porder+1)*ny*porder+1:(nx*porder+1)*(ny*porder+1)), :] .= -2 # force on the top
    elseif div(tid,100)==2
        EBC[collect(1:nx*porder+1:(nx*porder+1)*(ny*porder+1)), :] .= -1 # fix left
        FBC[collect(nx*porder+1:nx*porder+1:(nx*porder+1)*(ny*porder+1)), :] .= -2 # force on the right
    elseif div(tid,100)==3
        EBC[collect(1:nx*porder+1:(nx*porder+1)*(ny*porder+1)), :] .= -1 # fix left
        FBC[collect(2:nx*porder+1), :] .= -2 # force on the bottom
    end

    #
    if tid==100
        F1, F2 = ComputeLoad(Lx, nx, porder, ngp, "Constant",  [0, P1])
        fext[collect((nx*porder+1)*ny*porder+1:(nx*porder+1)*(ny*porder+1)), 1] .= F1
        fext[collect((nx*porder+1)*ny*porder+1:(nx*porder+1)*(ny*porder+1)), 2] .= F2

    elseif tid == 200
        F1, F2 = ComputeLoad(Ly, ny, porder, ngp, "Constant",  [P1, 0])
        fext[collect(nx*porder+1:nx*porder+1:(nx*porder+1)*(ny*porder+1)), 1] .= F1
        fext[collect(nx*porder+1:nx*porder+1:(nx*porder+1)*(ny*porder+1)), 2] .= F2

    elseif tid == 300

        F1, F2 = ComputeLoad(Lx, nx, porder, ngp, "Gaussian",  [Lx*P2, Lx*5.0/6.0, 0.2*Lx])
        fext[collect(2:nx*porder+1), 1] .= F1[2:end]
        fext[collect(2:nx*porder+1), 2] .= F2[2:end]
    else
        error("tid = $tid is not understood")
    end

    dof_to_active = findall(FBC[:].==-2)
    ft = t->fext[:][dof_to_active]*sin(π*t/(T))

    npoints = (nx+1)*(ny+1)
    node_to_point = zeros(Int64, nnodes) .- 1
    pointid = 1
    for j = 1:2:(ny*porder + 1)
        for i = 1:2:(nx*porder + 1)
            node_to_point[(j-1)*(nx*porder + 1) + i] = pointid
            pointid += 1
        end
    end

    return nodes, EBC, g, gt, FBC, fext, ft, npoints, node_to_point
end

force_scale = 2.0 #50
tid =200
time = 0.03
fiber_size = 1
porder = 2
T = 0.2
NT = 200

P1 = 44800 /force_scale #gcm/ms^2 compress/pull
P2 =  4480 /force_scale #gcm/ms^2 bend
P3 = 16800 /force_scale

@show P1, P2, force_scale
ngp = 3

testtype = "PlaneStress"
prop = Dict("name"=> testtype, "rho"=> 0.800, "E"=>2.0e6, "nu"=>0.2)

nxc, nyc = 10,5
nx, ny =  nxc*fiber_size, nyc*fiber_size
Lx, Ly = 0.1, 0.05 #m
nodes, EBC, g, gt, FBC, fext, ft, npoints, node_to_point = BoundaryCondition(tid, nx, ny, porder, Lx, Ly; force_scale=force_scale)

function Edge_func(x, y, t, tid)
    zero = similar(x) * 0.0
    if (tid == 100)
        f1 = @. zero
        f2 = @. P1*sin(π*t/(T)) + zero
    end
    if (tid == 200)
        f1 = @. P1*sin(π*t/(T)) + zero
        f2 = @. zero
    end
    if (tid == 300)
        pn, x0, σ = [Lx*P2, Lx*5.0/6.0, 0.2*Lx]
        f1 = @. zero
        f2 = @. pn * 1.0/(sqrt(2*pi*σ^2)) * exp.(-0.5*(x .- x0).^2/σ^2)*sin(π*t/(T))
    end

    return [f1 f2]
end

if tid == 100
    Edge_Traction_Data = zeros(Int64, nx, 3)
    for i = 1:nx
        Edge_Traction_Data[i, :] .= i+(ny-1)*nx, 3, 100
    end
elseif tid == 200
    Edge_Traction_Data = zeros(Int64, ny, 3)
    for i = 1:ny
        Edge_Traction_Data[i, :] .= nx+(i-1)*nx, 2, 200
    end
elseif tid == 300
    Edge_Traction_Data = zeros(Int64, nx, 3)
    for i = 1:nx
        Edge_Traction_Data[i, :] .= i, 1, 300
    end
end

elements = []
for j = 1:ny
    for i = 1:nx
        n = (nx*porder+1)*(j-1)*porder + (i-1)porder+1
        #element (i,j)
        if porder == 1
            #   4 ---- 3
            #
            #   1 ---- 2

            elnodes = [n, n + 1, n + 1 + (nx + 1), n + (nx + 1)]
        elseif porder == 2
            #   4 --7-- 3
            #   8   9   6
            #   1 --5-- 2
            elnodes = [n, n + 2, n + 2 + 2*(2*nx+1),  n + 2*(2*nx+1), n+1, n + 2 + (2*nx+1), n + 1 + 2*(2*nx+1), n + (2*nx+1), n+1+(2*nx+1)]
        else
            error("polynomial order error, porder= ", porder)
        end

        coords = nodes[elnodes,:]
        push!(elements,FiniteStrainContinuum(coords,elnodes, prop, 3))
    end
end

ndofs = 2
domain = Domain(nodes, elements, ndofs, EBC, g, FBC, fext, Edge_Traction_Data)
setGeometryPoints!(domain, npoints, node_to_point)

state = zeros(domain.neqs)
∂u = zeros(domain.neqs)
globdat = GlobalData(state,zeros(domain.neqs), zeros(domain.neqs),∂u, domain.neqs, gt, ft, nothing, Edge_func)

assembleMassMatrix!(globdat, domain)
updateStates!(domain, globdat)

fold = zeros(Float64, size(domain.fext))
ID = domain.ID
nodal_force = ft(time) # user defined time-dependent boundary
dof_id = 0
#update fext for active nodes (length of neqs)
for idof = 1:domain.ndims
    for inode = 1:domain.nnodes
        if (domain.FBC[inode, idof] == -2)
            global dof_id
            dof_id += 1
            fold[ID[inode, idof]] = nodal_force[dof_id]
        end
    end
end

fnew = getEdgeForce(domain, globdat, time)

@show norm(fold), norm(fnew)
@show norm(fold - fnew)


text:
function ObservationDataSet(arg0::RinexHeader, arg1::SatelliteSystem, arg2::jint, arg3::AbsoluteDate, arg4::jdouble, arg5::List)
    return ObservationDataSet((RinexHeader, SatelliteSystem, jint, AbsoluteDate, jdouble, List), arg0, arg1, arg2, arg3, arg4, arg5)
end

function get_date(obj::ObservationDataSet)
    return jcall(obj, "getDate", AbsoluteDate, ())
end

function get_header(obj::ObservationDataSet)
    return jcall(obj, "getHeader", RinexHeader, ())
end

function get_observation_data(obj::ObservationDataSet)
    return jcall(obj, "getObservationData", List, ())
end

function get_prn_number(obj::ObservationDataSet)
    return jcall(obj, "getPrnNumber", jint, ())
end

function get_rcvr_clk_offset(obj::ObservationDataSet)
    return jcall(obj, "getRcvrClkOffset", jdouble, ())
end

function get_satellite_system(obj::ObservationDataSet)
    return jcall(obj, "getSatelliteSystem", SatelliteSystem, ())
end


text:
mutable struct DriverBehavior
    desired_speed        :: Float64
    desired_acceleraion  :: Float64
    desired_deceleration :: Float64
    curve_speed          :: Float64
    observe_speed_limits :: Float64
    distance_keeping     :: Float64
    lane_keeping         :: Float64
    speed_keeping        :: Float64
    lane_change_dynamic  :: Float64
    urge_to_overtake     :: Float64
    keep_right_rule      :: Float64
    respond_to_tailgate  :: Float64
    foresight_distance   :: Float64
    steering_distance    :: Float64
    use_of_indicator     :: Float64
    obey_traffic_signs   :: Bool
    obey_traffic_lights  :: Bool
end

const HASTY_DRIVER       = DriverBehavior(0.85,0.80,0.80,0.80,0.35,0.35,0.7,0.7,0.75,0.85,0.5,0.5,0.50,0.7,0.30,true,true)
const BRISK_DRIVER       = DriverBehavior(0.70,0.65,0.65,0.60,0.45,0.40,0.5,0.5,0.60,0.60,0.5,0.5,0.55,0.5,0.65,true,true)
const DEFAULT_DRIVER     = DriverBehavior(0.50,0.50,0.50,0.50,0.50,0.50,0.5,0.5,0.50,0.50,0.5,0.5,0.50,0.5,0.50,true,true)
const COMFORTABLE_DRIVER = DriverBehavior(0.45,0.40,0.40,0.40,0.60,0.55,0.5,0.5,0.45,0.40,0.5,0.5,0.60,0.4,0.75,true,true)
const INSECURE_DRIVER    = DriverBehavior(0.40,0.35,0.65,0.35,0.90,0.65,0.3,0.3,0.35,0.20,0.5,0.5,0.30,0.3,0.20,true,true)

function Base.copy(d::DriverBehavior)
    DriverBehavior(
        d.desired_speed,
        d.desired_acceleraion,
        d.desired_deceleration,
        d.curve_speed,
        d.observe_speed_limits,
        d.distance_keeping,
        d.lane_keeping,
        d.speed_keeping,
        d.lane_change_dynamic,
        d.urge_to_overtake,
        d.keep_right_rule,
        d.respond_to_tailgate,
        d.foresight_distance,
        d. steering_distance,
        d. use_of_indicator,
        d.obey_traffic_signs,
        d.obey_traffic_lights
    )
end

max_stars_repo_name: sisl/VirtualTestDrive

max_issues_repo_name: sisl/VirtualTestDrive

max_forks_repo_name: sisl/VirtualTestDrive


text:
import ChainRulesCore: frule, rrule
using LinearAlgebra

function compute_gn_jac(u_n::AbstractVector{T}) where T
    n = length(u_n) - 1
    g_n = zero(u_n)
    # using `similar` in case u_n is e.g., `CuArray`
    dgdu = fill!(similar(u_n, n + 1, n + 1), zero(T))

    ## First: calculate the a_n terms (eqn. 10)
    a_n = zero(g_n)
    dadu = copy(dgdu)
    # constant term
    a_n[begin] = one(T)
    for i in 1:n
        bcoeff = one(T) # binomial coefficient
        for j in 0:i
            r = (-1)^j * bcoeff
            a_n[begin + j] -= r * u_n[begin + i]
            dadu[begin + i, begin + j] -= r
            bcoeff *= (i - j) / (j + 1)
        end
    end

    ## Second: calculate the g_n terms (eqn. 15)
    for j in n:-1:2
        if j ≥ n - 1
            g_n[begin + j] = a_n[begin + j] / (j + 2)
            for i in 1:n
                dgdu[begin + i, begin + j] = dadu[begin + i, begin + j] / (j + 2)
            end
        else
            g_n[begin + j] = a_n[begin + j] / (j + 2) + g_n[begin + j + 2]
            for i in 1:n
                dgdu[begin + i, begin + j] = dadu[begin + i, begin + j] / (j + 2) + dgdu[begin + i, begin + j + 2]
            end
        end
    end
    if n ≥ 3
        g_n[begin + 1] = a_n[begin + 1] + 3 * g_n[begin + 3]
        for i in 1:n
            dgdu[begin + i, begin + 1] = dadu[begin + i, begin + 1] + 3 * dgdu[begin + i, begin + 3]
        end
    elseif n ≥ 1
        g_n[begin + 1] = a_n[begin + 1]
        for i in 1:n
            dgdu[begin + i, begin + 1] = dadu[begin + i, begin + 1]
        end
    end
    if n ≥ 2
        g_n[begin] = a_n[begin] + 2 * g_n[begin + 2]
        for i in 1:n
            dgdu[begin + i, begin] = dadu[begin + i, begin] + 2 * dgdu[begin + i, begin + 2]
        end
    else
        g_n[begin] = a_n[begin]
        for i in 1:n
            dgdu[begin + i, begin] = dadu[begin + i, begin]
        end
    end

    return g_n, dgdu
end

function compute_uniform_grad(b::T, r; r2, b2, sqarea, fourbrinv) where T
    if b ≤ 1 - r
        flux = π * (1 - r2)
        kap0 = convert(T, π)
        kck = zero(T)
        kite_area2 = zero(T)
        ∇s0 = SA[zero(T), -2 * π * r]
    else
        kite_area2 = sqrt(sqarea)
        r2m1 = (r - 1) * (r + 1)
        kap0 = atan(kite_area2, r2m1 + b2)
        Πmkap1 = atan(kite_area2, r2m1 - b2)
        kck = kite_area2 * fourbrinv
        flux = Πmkap1 - r2 * kap0 + 0.5 * kite_area2
        ∇s0 = SA[kite_area2 / b, -2 * r * kap0]
    end

    return (flux, kap0, kite_area2, kck), ∇s0
end

function compute_linear_grad(b::T, r; k2, kc, kc2, onembmr2, onembpr2, onembmr2inv, sqonembmr2, r2=r^2, b2=b^2, br=b * r, fourbr=4 * br, sqbr=sqrt(br)) where T
    if iszero(b) # case 10
        Λ1 = -2 * π * sqrt(1 - r2)^3
        Eofk = 0.5 * π
        Em1mKdm = 0.25 * π

        ∇s1 = SA[zero(T), r * Λ1]
    elseif b == r
        if r == 0.5 # case 6
            Λ1 = π - 4 / 3
            Eofk = one(T)
            Em1mKdm = one(T)
            ∇s1 = SA[2 / 3, -2]
        elseif r < 0.5 # case 5
            m = 4 * r2
            Eofk = cel(m, one(T), one(T), 1 - m)
            Em1mKdm = cel(m, one(T), one(T), zero(T))
            Λ1 = π + 2 / 3 * ((2 * m - 3) * Eofk - m * Em1mKdm)
                # + (b - r) * 4 * r * (Eofk - 2 * Em1mKdm)
            ∇s1 = SA[-4 / 3 * r * (Eofk - 2 * Em1mKdm), -4 * r * Eofk]
        else # case 7
            m = 4 * r2
            minv = inv(m)
            Eofk = cel(minv, one(T), one(T), 1 - minv)
            Em1mKdm = cel(minv, one(T), one(T), zero(T))
            Λ1 = π + ((2 * m - 3) * Em1mKdm - m * Eofk) / (3 * r)
            #  -(b - r) * 2 * (2 * Eofk - Em1mKdm)
            ∇s1 = SA[2 / 3 * (2 * Eofk - Em1mKdm), -2 * Em1mKdm]
        end
    else
        if b + r > 1 # case 2, case 8
            Πofk, Eofk, Em1mKdm = cel(k2, kc, (b - r)^2 * kc2, zero(T), one(T), one(T), 3 * kc2 * (b - r) * (b + r), kc2, zero(T))
            Λ1 = onembmr2 * (Πofk + (-3 + 6 * r2 + 2 * br) * Em1mKdm - fourbr * Eofk) / (3 * sqbr)
            ∇s1 = SA[2 * r * onembmr2 * (-Em1mKdm + 2 * Eofk) / (sqbr * 3), -2 * r * onembmr2 * Em1mKdm / sqbr]
        elseif b + r < 1 # case 3 case 9
            bmrdbpr = (b - r) / (b + r)
            μ = 3 * bmrdbpr * onembmr2inv
            p = bmrdbpr^2 * onembpr2 * onembmr2inv
            Πofk, Eofk, Em1mKdm = cel(inv(k2), kc, p, 1 + μ, one(T), one(T), p + μ, kc2, zero(T))
            Λ1 = 2 * sqonembmr2 * (onembpr2 * Πofk - (4 - 7 * r2 - b2) * Eofk) / 3
            ∇s1 = SA[-4 / 3 * r * sqonembmr2 * (Eofk - 2 * Em1mKdm), -4 * r * sqonembmr2 * Eofk]
        else # case 4
            Λ1 = 2 * acos(1 - 2 * r) - 2 * π * (r > 0.5) - (4 / 3 * (3 + 2 * r - 8 * r2) + 8 * (r + b - 1) * r) * sqrt(r * (1 - r))
            Eofk = one(T)
            Em1mKdm = one(T)
            g = -8 * r * sqrt(r * (1 - r))
            ∇s1 = SA[-g / 3, g]
        end
    end

    flux = ((1 - T(r > b)) * 2π - Λ1) / 3

    return (flux, Eofk, Em1mKdm), ∇s1
end

function compute_quadratic_grad(b, r; s0, r2, b2, kap0, kite_area2, k2, ∇s0)
    r2pb2 = r2 + b2
    η2 = r2 * (r2pb2 + b2)
    if k2 > 1
        four_pi_eta = 2 * π * (η2 - 1)
        ∇η = 4 * π * SA[2 * b * r2, 2 * r * r2pb2]
    else
        Πmkap1 = atan(kite_area2, (r - 1) * (r + 1) - b2)
        four_pi_eta = 2 * (-Πmkap1 + η2 * kap0 - 0.25 * kite_area2 * (1 + 5 * r2 + b2))
        ∇η = SA[2 / b * (4 * b2 * r2 * kap0 - (1 + r2pb2) * kite_area2),
                8 * r * (r2pb2 * kap0 - kite_area2)]
    end
    s2 = 2 * s0 + four_pi_eta
    ∇s2 = 2 * ∇s0 + ∇η
    return  s2, ∇s2
end

function compute_grad(ld::PolynomialLimbDark, b::S, r) where S
    T = float(S)
    bcut = 1e-3
    n = ld.n_max
    dfdg = zero(ld.g_n)

    ## check for trivial cases
    if b ≥ 1 + r || iszero(r)
        # completely unobscured
        return one(T), dfdg, zero(T), zero(T)
    elseif r ≥ 1 + b
        # completely obscured
        return zero(T), dfdg, zero(T), zero(T)
    end

    r2 = r^2
    b2 = b^2

    if iszero(b)
        onemr2 = 1 - r2
        sqrt1mr2 = sqrt(onemr2)

        # annular ellipse
        facd = -2 * r
        if ld.n_max == 0
            flux = ld.g_n[begin] * onemr2
            dfdr = ld.g_n[begin] * facd
        else
            flux = ld.g_n[begin] * onemr2 + 2 / 3 * ld.g_n[begin + 1] * sqrt1mr2^3
            dfdr = ld.g_n[begin] * facd + ld.g_n[begin + 1] * facd * sqrt1mr2
        end
        fac = 2 * r2 * onemr2
        for i in 2:ld.n_max
            flux -= ld.g_n[begin + i] * fac
            dfdr += ld.g_n[begin + i] * facd * (2 * onemr2 - i * r2)
            dfdg[begin + i] -= fac
            fac *= sqrt1mr2
            facd *= sqrt1mr2
        end
        dfdr *= π * ld.norm
        dfdg[1] = (onemr2 - flux) * π * ld.norm
        dfdg[2] = 2 / 3 * (sqrt1mr2^3 - flux) * π * ld.norm
        return flux * π * ld.norm, dfdg * ld.norm, zero(T), dfdr
    end

    # take a moment to calculate values used repeatedly in the
    # numerical routines ahead
    onembmr2 = (r - b + 1) * (1 - r + b)
    onembmr2inv = inv(onembmr2)
    sqonembmr2 = sqrt(onembmr2)
    br = b * r
    fourbr = 4 * br
    fourbrinv = inv(fourbr)
    sqbr = sqrt(br)
    sqbrinv = inv(sqbr)
    onembpr2 = (1 - r - b) * (1 + b + r)
    sqarea = sqarea_triangle(one(T), r, b)
    k2 = max(zero(T), onembmr2 * fourbrinv)
    k = sqrt(k2)
    omemr2 = (1 - r) * (1 + r)
    onemr2mb2 = omemr2 - b2
    onemr2pb2 = omemr2 + b2
    if k2 ≥ 1
        if k2 > 2
            kc2 = 1 - inv(k2)
        else
            kc2 = onembpr2 * onembmr2inv
        end
    else
        if k2 > 0.5
            kc2 = (r - 1 + b) * (b + r + 1) * fourbrinv
        else
            kc2 = 1 - k2
        end
    end
    kc = sqrt(abs(kc2))

    ## Compute uniform term

    (s0, kap0, kite_area2, kck), ∇s0 = compute_uniform_grad(b, r; b2, r2, sqarea, fourbrinv)

    flux = ld.g_n[begin] * s0
    ∇flux = ld.g_n[begin] * ∇s0
    dfdg[begin] = s0

    if ld.n_max == 0
        dfdg[begin] -= flux * ld.norm * π
        dfdb, dfdr = ∇flux * ld.norm
        return flux * ld.norm, dfdg * ld.norm, dfdb, dfdr
    end

    ## compute linear term
    (s1, Eofk, Em1mKdm), ∇s1 = compute_linear_grad(b, r; k2, kc, kc2, r2, b2, br, fourbr, sqbr, onembmr2, onembpr2, onembmr2inv, sqonembmr2)

    flux += ld.g_n[begin + 1] * s1
    ∇flux += ld.g_n[begin + 1] * ∇s1
    dfdg[begin + 1] = s1

    if ld.n_max == 1
        dfdg[begin] -= flux * ld.norm * π
        dfdg[begin + 1] -= flux * ld.norm * π * 2 / 3
        dfdb, dfdr = ∇flux * ld.norm
        return flux * ld.norm, dfdg * ld.norm, dfdb, dfdr
    end

    ## calculate quadratic term
    s2, ∇s2 = compute_quadratic_grad(b, r; s0, r2, b2, kap0, kite_area2, k2, ∇s0)
    flux += ld.g_n[begin + 2] * s2
    ∇flux = ∇flux + ld.g_n[begin + 2] * ∇s2
    dfdg[begin + 2] = s2
    if ld.n_max == 2
        dfdg[begin] -= flux * ld.norm * π
        dfdg[begin + 1] -= flux * ld.norm * π * 2 / 3
        dfdb, dfdr = ∇flux * ld.norm
        return flux * ld.norm, dfdg * ld.norm, dfdb, dfdr
    end

    ## higher order (N > 3) terms require solving Mn and Nn integrals
    if k2 < 0.5 && ld.n_max > 3
        downwardM!(ld.Mn, ld.Mn_coeff; n_max=ld.n_max, sqonembmr2, sqbr, onemr2mb2, k, k2, sqarea, kap0, Eofk, Em1mKdm, kite_area2)
        if b < bcut
            Nn_lower!(ld.Nn, ld.Nn_coeff; Mn=ld.Mn, n_max=ld.n_max, onembpr2, sqonembmr2, kap0, kc, k2, sqbr, k, Eofk, Em1mKdm)
        end
    else
        upwardM!(ld.Mn; sqbr, n_max=ld.n_max, sqonembmr2, onemr2mb2, sqarea, k2, kap0, Eofk, Em1mKdm, kite_area2)
        if b < bcut
            Nn_raise!(ld.Nn, ld.Mn; n_max=ld.n_max, onembpr2, kap0, kc, k2, sqbr, k, Eofk, Em1mKdm)
        end
    end

    # compute remaining terms
    for n in 3:ld.n_max
        pofgn_M = 2 * r2 * ld.Mn[begin + n] - n / (n + 2) *
                       (onemr2mb2 * ld.Mn[begin + n] + sqarea * ld.Mn[begin + n - 2])
        flux -= ld.g_n[begin + n] * pofgn_M
        dfdg[begin + n] = -pofgn_M
        dpdr_M = 2 * r * ((n + 2) * ld.Mn[begin + n] - n * ld.Mn[begin + n - 2])

        if b < bcut
            dpdb_M = n * (ld.Mn[begin + n - 2] * (2 * r^3 + b^3 - b - 3 * r2 * b) + b * ld.Mn[begin + n] - 4 * r^3 * ld.Nn[begin + n - 2])
        else
            dpdb_M = n / b * ((ld.Mn[begin + n] - ld.Mn[begin + n - 2]) * (r2 + b2) + (b2 - r2)^2 * ld.Mn[begin + n - 2])
        end
        ∇flux -= ld.g_n[begin + n] * SA[dpdb_M, dpdr_M]

    end
    dfdg[begin] -= flux * ld.norm * π
    dfdg[begin + 1] -= flux * ld.norm * π * 2 / 3
    dfdb, dfdr = ∇flux * ld.norm
    return flux * ld.norm, dfdg * ld.norm, dfdb, dfdr
end

function frule((_, Δld, Δb, Δr), ::typeof(compute), ld::LD, b, r) where {LD <: PolynomialLimbDark}
    f, dfdg, dfdb, dfdr = compute_grad(ld, b, r)
    ∂g_n = dot(dfdg, Δld.g_n)
    return f, ∂g_n + dfdb * Δb + dfdr * Δr
end

function rrule(::typeof(compute), ld::LD, b, r) where {LD <: PolynomialLimbDark}
    f, dfdg, dfdb, dfdr = compute_grad(ld, b, r)
    function compute_pullback(Δf)
        ∂ld = Tangent{LD}(g_n=dfdg * Δf)
        ∂b = dfdb * Δf
        ∂r = dfdr * Δf
        return NoTangent(), ∂ld, ∂b, ∂r
    end
    return f, compute_pullback
end

function frule((_, Δu_n), ::Type{<:PolynomialLimbDark}, u::AbstractVector{S}; maxiter=100) where S
    T = float(S)
    # add constant u_0 term
    n_max = length(u)
    u_n = pushfirst!(copy(u), -one(S))
    # get Green's basis coefficients
    g_n, ∇g_n = compute_gn_jac(u_n)
    # compute series expansion for M_n and N_n
    Mn_coeff = compute_Mn_coeffs(T, n_max; maxiter=maxiter)
    Nn_coeff = compute_Nn_coeffs(T, n_max; maxiter=maxiter)

    # calculate flux normalization factor, which only depends on first two terms
    if length(g_n) > 1
        norm = inv(π * (g_n[1] + 2 / 3 * g_n[2]))
    else
        norm = inv(π * g_n[1])
    end

    # pre-allocate temp arrays
    Mn = similar(g_n)
    Nn = similar(g_n)

    Ω = PolynomialLimbDark(n_max, u_n, g_n, Mn_coeff, Nn_coeff, norm, Mn, Nn)
    ∂g_n = @views ∇g_n[begin + 1:end, :]' * Δu_n
    ∂Ω = Tangent{typeof(Ω)}(g_n=∂g_n)
    return Ω, ∂Ω
end

function rrule(::Type{<:PolynomialLimbDark}, u::AbstractVector{S}; maxiter=100) where S
    T = float(S)
    # add constant u_0 term
    n_max = length(u)
    u_n = pushfirst!(copy(u), -one(S))
    # get Green's basis coefficients
    g_n, ∇g_n = compute_gn_jac(u_n)
    # compute series expansion for M_n and N_n
    Mn_coeff = compute_Mn_coeffs(T, n_max; maxiter=maxiter)
    Nn_coeff = compute_Nn_coeffs(T, n_max; maxiter=maxiter)

    # calculate flux normalization factor, which only depends on first two terms
    if length(g_n) > 1
        norm = inv(π * (g_n[1] + 2 / 3 * g_n[2]))
    else
        norm = inv(π * g_n[1])
    end

    # pre-allocate temp arrays
    Mn = similar(g_n)
    Nn = similar(g_n)

    Ω = PolynomialLimbDark(n_max, u_n, g_n, Mn_coeff, Nn_coeff, norm, Mn, Nn)
    ∂g_n = @view ∇g_n[begin + 1:end, :]

    proj = ProjectTo(u)

    function PolynomialLimbDark_pullback(Δld)
        ∂u = ∂g_n * Δld.g_n
        return NoTangent(), proj(∂u)
    end
    return Ω, PolynomialLimbDark_pullback
end