text:
module Elements

export Element, Singleton, Group, kind, @kind, circulation, flux, streamfunction,
       complexpotential, conftransform, inverse_conftransform, jacobian, image

using ..Properties

abstract type Element end

#== Trait definitions ==#

abstract type Singleton end
abstract type Group end

function kind end

macro kind(element, k)
    esc(quote
        Elements.kind(::$element) = $k
        Elements.kind(::Type{$element}) = $k
        end)
end

@kind ComplexF64 Singleton
kind(::AbstractArray{T}) where {T <: Union{Element, ComplexF64}} = Group
kind(::Tuple) = Group

# Convenience functions to define wrapper types
# e.g. vortex sheets as a wrapper around vortex blobs
unwrap_src(e) = unwrap(e)
unwrap_targ(e) = unwrap(e)
unwrap(e) = e

## Actual Definitions of Properties

"""
    Elements.position(src::Element)

Returns the complex position of a potential flow element.
This is a required method for all `Element` types.

# Example

```jldoctest
julia> point = Vortex.Point(1.0 + 0.0im, 1.0);

julia> Elements.position(point)
1.0 + 0.0im

julia> points = Vortex.Point.([1.0im, 2.0im], 1.0);

julia> Elements.position.(points)
2-element Array{Complex{Float64},1}:
 0.0 + 1.0im
 0.0 + 2.0im
```
"""
@property begin
    signature = position(src::Source)
    stype = ComplexF64
end

"""
    Elements.circulation(src)

Returns the total circulation contained in `src`.

# Example

```jldoctest
julia> points = Vortex.Point.([1.0im, 2.0im], [1.0, 2.0]);

julia> blobs = Vortex.Blob.([1.0im, 2.0im], [1.0, 2.0], 0.1);

julia> Elements.circulation(points[1])
1.0

julia> Elements.circulation(points)
3.0

julia> Elements.circulation((points, blobs))
6.0

julia> Elements.circulation.(points)
2-element Array{Float64,1}:
 1.0
 2.0

julia> Elements.circulation.((points, blobs))
(3.0, 3.0)

julia> Elements.circulation(Source.Point(rand(), rand()))
0.0

julia> Elements.circulation(Source.Blob(rand(), rand(), rand()))
0.0
```
"""
@property begin
    signature = circulation(src::Source)
    reduce = (+)
    stype = Float64
end

"""
    Elements.flux(src)

Returns the flux through a unit circle induced by `src`.

# Example

```jldoctest
julia> points = Source.Point.([1.0im, 2.0im], [1.0, 2.0]);

julia> blobs = Source.Blob.([1.0im, 2.0im], [1.0, 2.0], 0.1);

julia> Elements.flux(points[1])
1.0

julia> Elements.flux((points, blobs))
6.0

julia> Elements.flux.(points)
2-element Array{Float64,1}:
 1.0
 2.0

julia> Elements.flux.((points, blobs))
(3.0, 3.0)

julia> Elements.flux(Vortex.Point(rand(), rand()))
0.0

julia> Elements.flux(Vortex.Blob(rand(), rand(), rand()))
0.0
```
"""
@property begin
    signature = flux(src::Source)
    reduce = (+)
    stype = Float64
end

@doc raw"""
    Elements.impulse(src)

Return the aerodynamic impulse of `src` about (0,0):
```math
P := \int \boldsymbol{x} \times \boldsymbol{\omega}\,\mathrm{d}A.
```
This is a required method for all vortex types.

# Example

```jldoctest
julia> sys = (Vortex.Point(1.0im, π), Vortex.Blob(1.0im, -π, 0.1));

julia> Elements.impulse(sys[1])
3.141592653589793 + 0.0im

julia> Elements.impulse(sys)
0.0 + 0.0im
```
"""
@property begin
    signature = impulse(src::Source)
    reduce = (+)
    stype = ComplexF64
end

@doc raw"""
    Elements.angularimpulse(src)

Return the aerodynamic angular impulse of `src` about (0,0):
```math
\Pi := \frac{1}{2}\int \boldsymbol{x} \times (\boldsymbol{x}\times\boldsymbol{\omega})\,\mathrm{d}A.
```
This is a required method for all vortex types.

# Example

```jldoctest
julia> sys = (Vortex.Point(1.0im, π), Vortex.Blob(1.0im, -π, 0.1));

julia> Elements.angularimpulse(sys[1])
3.141592653589793 + 0.0im

julia> Elements.angularimpulse(sys)
0.0 + 0.0im
```
"""
@property begin
    signature = angularimpulse(src::Source)
    reduce = (+)
    stype = ComplexF64
end

@property begin
    signature = streamfunction(targ::Target, src::Source)
    preallocator = allocate_streamfunction
    stype = Float64
end

@property begin
    signature = complexpotential(targ::Target, src::Source)
    preallocator = allocate_complexpotential
    stype = ComplexF64
end

@doc raw"""
    Elements.conftransform(src,body)

Return the conformally transformed position of `src` via the transform defined
by `body`. In this context, the position(s) in `src` is interpreted as
coordinates in the circle plane.

# Example

```jldoctest
julia> sys = (Vortex.Point(1.0im, π), Vortex.Blob(2.0im, -π, 0.1));

julia> b = ConformalBody([1/4,0,1/4]);

julia> Elements.conftransform(sys,b)
(0.0 + 0.0im, 0.0 + 0.375im)
```
"""
@property begin
    signature = conftransform(src::Source,b)
    preallocator = allocate_conftransform
    stype = ComplexF64
end

@doc raw"""
    Elements.inverse_conftransform(src,body)

Return the inverse conformally transformed position of `src` via the transform defined
by `body`. In this context, the position(s) in `src` is interpreted as
coordinates in the physical plane.

# Example

```jldoctest
julia> sys = (Vortex.Point(1.0im, π), Vortex.Blob(2.0im, -π, 0.1));

julia> b = ConformalBody([1/4,0,1/4]);

julia> Elements.inverse_conftransform(sys,b)
(0.0 + 0.0im, 0.0 + 0.375im)
```
"""
@property begin
    signature = inverse_conftransform(src::Source,b)
    preallocator = allocate_inv_conftransform
    stype = ComplexF64
end

@doc raw"""
    Elements.jacobian(src,body)

Return the Jacobian of the conformal transform at the position of `src` via the transform defined
by `body`. In this context, the position(s) in `src` is interpreted as
coordinates in the circle plane.

# Example

```jldoctest
julia> sys = (Vortex.Point(1.0im, π), Vortex.Blob(2.0im, -π, 0.1));

julia> b = ConformalBody([1/4,0,1/4]);

julia> Elements.jacobian(sys,b)
(0.5 + 0.0im, 0.0 + 0.5im)

```
"""
@property begin
    signature = jacobian(src::Source,b)
    preallocator = allocate_jacobian
    stype = ComplexF64
end

@doc raw"""
    Elements.image(src,body)

Return the image position of `src` in the circle-plane representation of `body`.
In this context, the position(s) in `src` is interpreted as
coordinates in the circle plane.

# Example

```jldoctest
julia> sys = (Vortex.Point(1.0im, π), Vortex.Blob(2.0im, -π, 0.1));

julia> b = PowerBody([1/4,0,1/4],zero(ComplexF64),0.0)

julia> Elements.image(sys,b)
(0.0 + 1.0im, 0.0 + 0.5im)

```
"""
@property begin
    signature = image(src::Source,b)
    preallocator = allocate_image
    stype = ComplexF64
end

end


text:
"The simplest viable reverse mode a AD, only supports `Float64`"
module ReverseDiffZero
using ChainRulesCore
using Test

#########################################
# Initial rule setup
@scalar_rule x + y (1, 1)
@scalar_rule x - y (1, -1)
##########################
#Define the AD

struct Tracked{F} <: Real
    propagate::F
    primal::Float64
    tape::Vector{Tracked}  # a reference to a shared tape
    partial::Base.RefValue{Float64} # current accumulated sensitivity
end

"An intermediate value, a Branch in Nabla terms."
function Tracked(propagate, primal, tape)
    v = Tracked(propagate, primal, tape, Ref(zero(primal)))
    push!(tape, v)
    return v
end

"Marker for inputs (leaves) that don't need to propagate."
struct NoPropagate end

"An input, a Leaf in Nabla terms. No inputs of its own to propagate to."
function Tracked(primal, tape)
    # don't actually need to put these on the tape, since they don't need to propagate
    return Tracked(NoPropagate(), primal, tape, Ref(zero(primal)))
end

primal(d::Tracked) = d.primal
primal(d) = d

partial(d::Tracked) = d.partial[]
partial(d) = nothing

tape(d::Tracked) = d.tape
tape(d) = nothing

"we have many inputs grab the tape from the first one that is tracked"
get_tape(ds) = something(tape.(ds)...)

"propagate the currently stored partial back to my inputs."
propagate!(d::Tracked) = d.propagate(d.partial[])

"Accumulate the sensitivity, if the value is being tracked."
accum!(d::Tracked, x̄) = d.partial[] += x̄
accum!(d, x̄) = nothing

# needed for `^` to work from having `*` defined
Base.to_power_type(x::Tracked) = x

"What to do when a new rrule is declared"
function define_tracked_overload(sig)
    sig = Base.unwrap_unionall(sig)  # not really handling most UnionAll
    opT, argTs = Iterators.peel(sig.parameters)
    opT isa Type{<:Type} && return  # not handling constructors
    fieldcount(opT) == 0 || return  # not handling functors
    all(argT isa Type && Float64 <: argT for argT in argTs) || return # only handling purely Float64 ops.

    N = length(sig.parameters) - 1  # skip the op
    fdef = quote
        # we use the function call overloading form as it lets us avoid namespacing issues
        # as we can directly interpolate the function type into to the AST.
        function (op::$opT)(tracked_args::Vararg{Union{Tracked, Float64}, $N}; kwargs...)
            args = (op, primal.(tracked_args)...)
            y, y_pullback = rrule(args...; kwargs...)
            the_tape = get_tape(tracked_args)
            y_tracked = Tracked(y, the_tape) do ȳ
                # pull this partial back and propagate it to the input's partial store
                _, ārgs = Iterators.peel(y_pullback(ȳ))
                accum!.(tracked_args, ārgs)
            end
            return y_tracked
        end
    end
    eval(fdef)
end

# !Important!: Attach the define function to the `on_new_rule` hook
on_new_rule(define_tracked_overload, rrule)

"Do a calculus. `f` should have a single output."
function derv(f, args::Vararg; kwargs...)
    the_tape = Vector{Tracked}()
    tracked_inputs = Tracked.(args, Ref(the_tape))
    tracked_output = f(tracked_inputs...; kwargs...)
    @assert tape(tracked_output) === the_tape

    # Now the backward pass
    out = primal(tracked_output)
    ōut = one(out)
    accum!(tracked_output, ōut)
    # By going down the tape backwards we know we will have fully accumulated partials
    # before propagating them onwards
    for op in reverse(the_tape)
        propagate!(op)
    end
    return partial.(tracked_inputs)
end

# End AD definition
################################

# add a rule later also
function ChainRulesCore.rrule(::typeof(*), x::Number, y::Number)
    function times_pullback(ΔΩ)
        # we will use thunks here to show we handle them fine.
        return (NO_FIELDS,  @thunk(ΔΩ * y'), @thunk(x' * ΔΩ))
    end
    return x * y, times_pullback
end

# Manual refresh needed as new rule added in same file as AD after the `on_new_rule` call
refresh_rules();

@testset "ReversedDiffZero" begin
    foo(x) = x + x
    @test derv(foo, 1.6) == (2.0,)

    bar(x) = x + 2.1 * x
    @test derv(bar, 1.2) == (3.1,)

    baz(x) = 2.0 * x^2 + 3.0*x + 1.2
    @test derv(baz, 1.7) == (2 * 2.0 * 1.7 + 3.0,)

    qux(x) = foo(x) + bar(x) + baz(x)
    @test derv(qux, 1.7) == ((2 * 2.0 * 1.7 + 3.0) + 3.1 + 2,)

    function quux(x)
        y = 2.0*x + 3.0*x
        return 4.0*y + 5.0*y
    end
    @test derv(quux, 11.1) == (4*(2+3) + 5*(2+3),)
end
end  # module

ext: jl


text:
module TestRelationalPrograms
using Test

using Catlab.CategoricalAlgebra.CSets
using Catlab.WiringDiagrams.UndirectedWiringDiagrams
using Catlab.Programs.RelationalPrograms

# Relation macro
################

# Untyped
#--------

parsed = @relation (x,z) where (x,y,z) begin
  R(x,y)
  S(y,z)
end

d = RelationDiagram{Symbol}(2)
add_box!(d, 2, name=:R); add_box!(d, 2, name=:S)
add_junctions!(d, 3, variable=[:x,:y,:z])
set_junction!(d, [1,2,2,3])
set_junction!(d, [1,3], outer=true)
@test parsed == d

parsed = @relation ((x,z) where (x,y,z)) -> (R(x,y); S(y,z))
@test parsed == d
parsed = @relation function (x,z) where (x,y,z); R(x,y); S(y,z) end
@test parsed == d

# Abbreviated syntax where context is inferred.
d1 = @relation (x,y,z) -> (R(x,y); S(y,z))
d2 = @relation ((x,y,z) where (x,y,z)) -> (R(x,y); S(y,z))
@test d1 == d2

d1 = @relation (x,z) -> (R(x,y); S(y,z))
d2 = @relation ((x,z) where (x,z,y)) -> (R(x,y); S(y,z))
@test d1 == d2

# Special case: closed diagram.
parsed = @relation (() where (a,)) -> R(a,a)
d = RelationDiagram{Symbol}(0)
add_box!(d, 2, name=:R)
add_junction!(d, variable=:a); set_junction!(d, [1,1])
@test parsed == d

parsed = @relation () -> A()
d = RelationDiagram{Symbol}(0)
add_box!(d, 0, name=:A)
@test parsed == d

# Typed
#------

parsed = @relation (x,y,z) where (x::X, y::Y, z::Z, w::W) begin
  R(x,w)
  S(y,w)
  T(z,w)
end
d = RelationDiagram{Symbol}([:X,:Y,:Z])
add_box!(d, [:X,:W], name=:R)
add_box!(d, [:Y,:W], name=:S)
add_box!(d, [:Z,:W], name=:T)
add_junctions!(d, [:X,:Y,:Z,:W], variable=[:x,:y,:z,:w])
set_junction!(d, [1,4,2,4,3,4])
set_junction!(d, [1,2,3], outer=true)
@test parsed == d

# Untyped, named ports
#---------------------

parsed = @relation (start=u, stop=w) where (u, v, w) begin
  E(src=u, tgt=v)
  E(src=v, tgt=w)
end
d = RelationDiagram{Symbol}(2, port_names=[:start, :stop])
add_box!(d, 2, name=:E)
add_box!(d, 2, name=:E)
add_junctions!(d, 3, variable=[:u,:v,:w])
set_junction!(d, [1,2,2,3])
set_junction!(d, [1,3], outer=true)
set_subpart!(d, :port_name, [:src, :tgt, :src, :tgt])
@test parsed == d

# Abbreviated syntax where context is inferred.
d1 = @relation (start=u, mid=v, stop=w) -> (E(src=u, tgt=v); E(src=v, tgt=w))
d2 = @relation ((start=u, mid=v, stop=w) where (u,v,w)) ->
  (E(src=u, tgt=v); E(src=v, tgt=w))
@test d1 == d2

d1 = @relation (start=u, stop=w) -> (E(src=u, tgt=v); E(src=v, tgt=w))
d2 = @relation ((start=u, stop=w) where (u,w,v)) ->
  (E(src=u, tgt=v); E(src=v, tgt=w))
@test d1 == d2

# Special case: unary named tuple syntax.
parsed = @relation ((src=v) where (v, w)) -> E(src=v, tgt=w)
@test parsed[:outer_port_name] == [:src]

# Special case: closed diagram.
if VERSION >= v"1.5"
  d1 = @relation ((;) where (v,)) -> E(src=v, tgt=v)
  @test subpart(d1, :port_name) == [:src, :tgt]

  d2 = @relation (;) -> E(src=v, tgt=v)
  @test d1 == d2
end

# Typed, named ports
#-------------------

parsed = @relation (e=e, e′=e′) where (e::Employee, e′::Employee,
                                       d::Department) begin
  Employee(id=e, department=d)
  Employee(id=e′, department=d)
end
d = RelationDiagram{Symbol}([:Employee, :Employee], port_names=[:e,:e′])
add_box!(d, [:Employee, :Department], name=:Employee)
add_box!(d, [:Employee, :Department], name=:Employee)
add_junctions!(d, [:Employee, :Employee, :Department], variable=[:e,:e′,:d])
set_junction!(d, [1,3,2,3])
set_junction!(d, [1,2], outer=true)
set_subpart!(d, :port_name, [:id, :department, :id, :department])
@test parsed == d

end


text:
module PLext

using Reexport

#include("simulator/statevector.jl")
include("simulator/statevector_ptr.jl")
#include("simulator/statevector_copy.jl")
#include("simulator/statevector_ptr_dist.jl")
#include("compile/sysimg.jl")
#include("compile/precompile.jl")
#include("compile/library.jl")

#@reexport using PLext.StateVectorP
#@reexport using PLext.StateVectorC

end

max_stars_repo_name: mlxd/PLext

max_issues_repo_name: mlxd/PLext

max_forks_repo_name: mlxd/PLext


text:
# Michelle Medina
# Finite Differencing
# DATE

export Params, Params2D
export dX, dX_2D, dY_2D, functionOnGrid, Deriv1D, DerivX_2D, DerivY_2D

struct Params
#structure contains # of points N and Domain Space L
    N::Int
    L::Float64
end

struct Params2D
    #contains number of points Nx and Ny and the length of the
    #2D grid Lx and Ly
    Nx::Int
    Ny::Int
    Lx::Float64
    Ly::Float64
end

# Defining the Grid Space
""" computing grid space h """
function dX(grid::Params)
    return (grid.L/(grid.N-1))
end

""" computing grid space hx """
function dX_2D(grid::Params2D)
    return (grid.Lx/(grid.Nx-1))
end

""" computing grid space hy """
function dY_2D(grid::Params2D)
    #computes grid space hy
    return (grid.Ly/(grid.Ny-1))
end

""" 1D Function is computed on the grid """
function functionOnGrid(u::Function, grid::Params)::Array{Float64,1}
    ux=zeros(grid.N)
    h=dX(grid)
    for i in 1:grid.N
        x=(i-1)*h
        ux[i]=u(x)
    end
    return ux
end

""" 2D Function is computed on the grid """
function functionOnGrid(u::Function, grid::Params2D)::Array{Float64,2}
    uxy=zeros(grid.Nx,grid.Ny)
    hx=dX_2D(grid)
    hy=dY_2D(grid)
    for i in 1:grid.Nx, j in 1:grid.Ny
        x = (i-1)*hx
        y= (j-1)*hy
        uxy[i,j]=u(x,y)
    end
    return uxy
end

""" Derivative of a function is computed on a 1D grid """
function Deriv1D(u::Array{Float64,1},grid::Params)::Array{Float64}
    dudx = zeros(grid.N)
    h = dX(grid)
    for i in 1:grid.N
        if i == 1
            dudx[i] = (u[i+1] - u[i])/h
        elseif i == grid.N
            dudx[i] = (u[i]-u[i-1])/h
        else
            dudx[i] = (u[i+1]-u[i-1])/2h
        end
    end
    return dudx
end

function DerivX_2D(u::Array{Float64,2}, grid::Params2D)::Array{Float64}
    dudx= zeros(grid.Nx,grid.Ny)
    hx = dX_2D(grid)

    for i in 1:grid.Nx, j in 1:grid.Ny
        if i==1
            dudx[i,j] = (u[i+1,j] - u[i,j])/hx
        elseif i==grid.Nx
            dudx[i,j] = (u[i,j] - u[i-1,j])/hx
        else
            dudx[i,j] = (u[i+1,j]-u[i-1,j])/2hx
        end
    end
    return dudx
end

function DerivY_2D(u::Array{Float64,2},grid::Params2D)::Array{Float64}
    dudy= zeros(grid.Nx,grid.Ny)
    hy = dY_2D(grid)

    for i in 1:grid.Nx, j in 1:grid.Ny
        if j == 1
            dudy[i,j] = (u[i,j+1] - u[i,j])/hy
        elseif j== grid.Ny
            dudy[i,j] = (u[i,j]-u[i,j-1])/hy
        else
            dudy[i,j] = (u[i,j+1]-u[i,j-1])/2hy
        end
    end
    return dudy
end


text:
# https://adrhill.github.io/ExplainableAI.jl/stable/generated/advanced_lrp/#How-it-works-internally
abstract type AbstractLRPRule end

# Generic LRP rule. Since it uses autodiff, it is used as a fallback for layer types without custom implementation.
function lrp!(Rₖ, rule::R, layer::L, aₖ, Rₖ₊₁) where {R<:AbstractLRPRule,L}
    lrp_autodiff!(Rₖ, rule, layer, aₖ, Rₖ₊₁)
    return nothing
end

function lrp_autodiff!(
    Rₖ::T1, rule::R, layer::L, aₖ::T1, Rₖ₊₁::T2
) where {R<:AbstractLRPRule,L,T1,T2}
    layerᵨ = modify_layer(rule, layer)
    c::T1 = only(
        gradient(aₖ) do a
            z::T2 = layerᵨ(a)
            s = Zygote.@ignore Rₖ₊₁ ./ modify_denominator(rule, z)
            z ⋅ s
        end,
    )
    Rₖ .= aₖ .* c
    return nothing
end

# For linear layer types such as Dense layers, using autodiff is overkill.
function lrp!(Rₖ, rule::R, layer::Dense, aₖ, Rₖ₊₁) where {R<:AbstractLRPRule}
    lrp_dense!(Rₖ, rule, layer, aₖ, Rₖ₊₁)
    return nothing
end

function lrp_dense!(Rₖ, rule::R, l, aₖ, Rₖ₊₁) where {R<:AbstractLRPRule}
    ρW, ρb = modify_params(rule, get_params(l)...)
    ãₖ₊₁ = modify_denominator(rule, ρW * aₖ .+ ρb)
    @tullio Rₖ[j, b] = aₖ[j, b] * ρW[k, j] / ãₖ₊₁[k, b] * Rₖ₊₁[k, b]
    return nothing
end

# Other special cases that are dispatched on layer type:
lrp!(Rₖ, ::AbstractLRPRule, ::DropoutLayer, aₖ, Rₖ₊₁) = (Rₖ .= Rₖ₊₁)
lrp!(Rₖ, ::AbstractLRPRule, ::ReshapingLayer, aₖ, Rₖ₊₁) = (Rₖ .= reshape(Rₖ₊₁, size(aₖ)))

# To implement new rules, we can define two custom functions `modify_params` and `modify_denominator`.
# If this isn't done, the following fallbacks are used by default:
"""
    modify_params(rule, W, b)

Function that modifies weights and biases before applying relevance propagation.
Returns modified weights and biases as a tuple `(ρW, ρb)`.
"""
modify_params(::AbstractLRPRule, W, b) = (W, b) # general fallback

"""
    modify_denominator(rule, d)

Function that modifies zₖ on the forward pass, e.g. for numerical stability.
"""
modify_denominator(::AbstractLRPRule, d) = stabilize_denom(d, 1.0f-9) # general fallback

"""
    modify_layer(rule, layer)

Function that modifies a layer before applying relevance propagation.
Returns a new, modified layer.
"""
modify_layer(::AbstractLRPRule, layer) = layer # skip layers without modify_params
function modify_layer(rule::R, layer::L) where {R<:AbstractLRPRule,L<:Union{Dense,Conv}}
    return set_params(layer, modify_params(rule, get_params(layer)...)...)
end

"""
    ZeroRule()

Constructor for LRP-0 rule. Commonly used on upper layers.
"""
struct ZeroRule <: AbstractLRPRule end

"""
    GammaRule([γ=0.25])

Constructor for LRP-``γ`` rule. Commonly used on lower layers.

Arguments:
- `γ`: Optional multiplier for added positive weights, defaults to 0.25.
"""
struct GammaRule{T} <: AbstractLRPRule
    γ::T
    GammaRule(γ=0.25f0) = new{Float32}(γ)
end
function modify_params(r::GammaRule, W, b)
    T = eltype(W)
    ρW = W + convert(T, r.γ) * relu.(W)
    ρb = b + convert(T, r.γ) * relu.(b)
    return ρW, ρb
end

"""
    EpsilonRule([ϵ=1.0f-6])

Constructor for LRP-``ϵ`` rule. Commonly used on middle layers.

Arguments:
- `ϵ`: Optional stabilization parameter, defaults to `1f-6`.
"""
struct EpsilonRule{T} <: AbstractLRPRule
    ϵ::T
    EpsilonRule(ϵ=1.0f-6) = new{Float32}(ϵ)
end
modify_denominator(r::EpsilonRule, d) = stabilize_denom(d, r.ϵ)

"""
    ZBoxRule(low, high)

Constructor for LRP-``z^{\\mathcal{B}}``-rule.
Commonly used on the first layer for pixel input.

The parameters `low` and `high` should be set to the lower and upper bounds of the input features,
e.g. `0.0` and `1.0` for raw image data.
It is also possible to provide two arrays of that match the input size.
"""
struct ZBoxRule{T} <: AbstractLRPRule
    low::T
    high::T
end

# The ZBoxRule requires its own implementation of relevance propagation.
lrp!(Rₖ, r::ZBoxRule, layer::Dense, aₖ, Rₖ₊₁) = lrp_zbox!(Rₖ, r, layer, aₖ, Rₖ₊₁)
lrp!(Rₖ, r::ZBoxRule, layer::Conv, aₖ, Rₖ₊₁) = lrp_zbox!(Rₖ, r, layer, aₖ, Rₖ₊₁)

_zbox_bound(T, c::Real, in_size) = fill(convert(T, c), in_size)
function _zbox_bound(T, A::AbstractArray, in_size)
    size(A) != in_size && throw(
        ArgumentError(
            "Bounds `low`, `high` of ZBoxRule should either be scalar or match input size.",
        ),
    )
    return convert.(T, A)
end

function lrp_zbox!(Rₖ::T1, r::ZBoxRule, layer::L, aₖ::T1, Rₖ₊₁::T2) where {L,T1,T2}
    T = eltype(aₖ)
    in_size = size(aₖ)
    l = _zbox_bound(T, r.low, in_size)
    h = _zbox_bound(T, r.high, in_size)

    W, b = get_params(layer)
    layer⁺ = set_params(layer, max.(0, W), max.(0, b)) # W⁺, b⁺
    layer⁻ = set_params(layer, min.(0, W), min.(0, b)) # W⁻, b⁻

    c::T1, cₗ::T1, cₕ::T1 = gradient(aₖ, l, h) do a, l, h
        f::T2 = layer(a)
        f⁺::T2 = layer⁺(l)
        f⁻::T2 = layer⁻(h)

        z = f - f⁺ - f⁻
        s = Zygote.@ignore safedivide(Rₖ₊₁, z, 1e-9)
        z ⋅ s
    end
    Rₖ .= aₖ .* c + l .* cₗ + h .* cₕ
    return nothing
end


text:
## Copyright (c) 2015, Júlio Hoffimann Mendes <juliohm@stanford.edu>
##
## Permission to use, copy, modify, and/or distribute this software for any
## purpose with or without fee is hereby granted, provided that the above
## copyright notice and this permission notice appear in all copies.
##
## THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
## WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
## MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
## ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
## WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
## ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
## OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.

@userplot VoxelReusePlot

@recipe function f(vr::VoxelReusePlot;
        tmin=nothing, tmax=nothing, nreal=10,
        gpu=false, soft=[], hard=HardData(), tol=.1,
        overlapx=1/6, overlapy=1/6, overlapz=1/6,
        cut=:boykov, simplex=false, threads=CPU_PHYSICAL_CORES)
  # get input image
  img = vr.args[1]

  @assert ndims(img) == 3 "image is not 3D (add ghost dimension for 2D)"

  # image extent
  extent = size(img)
  idx = [extent...] .> 1

  # default support
  tmin == nothing && (tmin = 7)
  tmax == nothing && (tmax = min(100, minimum(extent[idx])))

  @assert tmin > 0 "tmin must be positive"
  @assert tmin < tmax "tmin must be smaller than tmax"

  # save support as an array
  ts = collect(tmin:tmax)

  # -- compute voxel reuse
  extent = size(img)
  idx = [extent...] .> 1

  p = Progress(length(ts), color=:black)
  μs = Float64[]; σs = Float64[]
  for T in ts
    tplconfig = [1,1,1]
    tplconfig[idx] = T
    μ, σ = voxelreuse(img, tplconfig...,
            nreal=nreal, gpu=gpu, soft=soft, hard=hard,
			overlapx=overlapx, overlapy=overlapy, overlapz=overlapz,
			cut=cut, simplex=simplex, threads=threads)
    push!(μs, μ); push!(σs, σ)

    next!(p)
  end
  # --

  # highlight the optimum template range
  @series begin
    rank = sortperm(μs, rev=true)
    best = ts[rank[1:min(5,length(ts))]]
    xmin, xmax = minimum(best), maximum(best)

    reference = minimum(μs - σs)
    ymin = max(0., reference - .05)
    ymax = ymin + .008

    tplrange = [xmin, xmax]

    seriestype := :shape
    linewidth := 0
    fillalpha := .5
    label --> "Optimum range: $tplrange"

    [xmin, xmax, xmax, xmin], [ymin, ymin, ymax, ymax]
  end

  seriestype := :path
  primary := false
  ribbon := σs
  fillalpha := .5
  xlim := (0, tmax)
  xlabel --> "Template size"
  ylabel --> "Voxel reuse"

  ts, μs
end


text:
######################################################################
# ergodicity.jl
# Handles the stuff needed for ergodicity
######################################################################

# TODO: need initial robot positions to initialize ck
# Gamma		constant weights
# h			constant normalizing weights
# phik		Fourier coefficients for current EID
# phi
# ck
type ErgodicManager
	K::Int
	L::Int
	Gamma::Matrix{Float64}
	h::Matrix{Float64}
	phi::Matrix{Float64}
	phik::Matrix{Float64}
	ck::Matrix{Float64}
	F

	function ErgodicManager(m::SearchDomain, K::Int)
		L = m.num_cells
		F = m.F

		# create data structures
		Gamma = zeros(K+1, K+1)
		h     = zeros(K+1, K+1)
		phik  = zeros(K+1, K+1)
		ck    = zeros(K+1, K+1)
		phi	  = EID(m)

		pil2 = pi * pi / (L * L)
		for K1 = 0:K
			for K2 = 0:K
				Gamma[K1+1, K2+1] = 1 / ( (1 + pil2*(K1*K1 + K2*K2))^1.5 )
				h[K1+1, K2+1] = h_ij(K1, K2, L)
			end
		end
		return new(K, L, Gamma, h, phi, phik, ck, F)
	end
end

"""
`fk(em, k1, k2, x1, x2)`

Arguments:

* `em` = ergodic manager
* `k1, k2` = coefficients (0-indexed)
* `x1, x2` = vehicle location

"""
function fk(em::ErgodicManager, K1::Int, K2::Int, x1, x2)
	return cos(K1*pi*x1/em.L) * cos(K2*pi*x2/em.L) / em.h[K1+1, K2+1]
end

######################################################################
# Trajectory Coefficients c_k
######################################################################
# t = 0 is the start point ... this method would fail...
# t = 1 is when we can start using this function
function update_ck!(em::ErgodicManager, X::VehicleSet, t::Int)
	for K1 = 0:em.K
		for K2 = 0:em.K
			em.ck[K1+1,K2+1] = ck_ij(em, K1, K2, X, t)
		end
	end
end

# Computes update for c_{ij}
# TODO: handle mulitple vehicles
function ck_ij(em::ErgodicManager, K1::Int, K2::Int, X::VehicleSet, t::Int)
	N = length(X)
	new_sum = 0.0
	for i = 1:N
		new_sum += fk(em, K1, K2, X[i].x, X[i].y)
	end
	return (N*em.ck[K1+1,K2+1] + new_sum) / (N*t)
end

"""
`h_ij(k1::Int, k2::Int, L::Int)`

`L` is size of search domain.
"""
function h_ij(K1::Int, K2::Int, L::Int)
	# Limits
	xmin = (0.,0.)
	xmax = (L, L)

	# make integrand
	function h_int(x)
		cos(K1*pi*x[1]/L)^2 * cos(K2*pi*x[2]/L)^2
	end

	return sqrt(hcubature(h_int, xmin, xmax)[1])
end

# This actually works remarkably well
function hk2(m::SearchDomain, k::Int)
	L = m.num_cells
	hksum = 0.0
	for x1 = 0.5:(L-.5)
		for x2 = 0.5:(L-.5)
			hksum += cos(k*pi*x1/L)^2 * cos(k*pi*x2/L)^2
		end
	end
	return sqrt(hksum)
end

######################################################################
# wtf
######################################################################
function update_phik!(em::ErgodicManager)
	K = em.K
	for K1 = 0:K
		for K2 = 0:K
			em.phik[K1+1, K2+1] = 0.0
			for x = 1:em.L
				for y = 1:em.L
					em.phik[K1+1, K2+1] += em.phi[x,y]*fk(em,K1,K2,x,y)
				end
			end
		end
	end
end

"""
`phik(m::SearchDomain, phi::Matrix{Float64}, k::Int)`

Currently performs a rough numerical integration.
"""
function phi_ij(em::ErgodicManager, K1::Int, K2::Int)
	L = em.L
	phiksum = 0.0
	k1 = K1 * pi / L
	k2 = K2 * pi / L
	for x1 = 1:L
		for x2 = 1:L
			phiksum += em.phi[x1, x2] * cos(k1*(x1-.5)) * cos(k2*(x2-.5))
		end
	end
	return phiksum / em.h[K1+1, K2+1]
end

# Returns bx, by, the two components of Bj
function Bj(em::ErgodicManager, x1, x2, N, t)
	bx = 0.
	by = 0.
	for K1 = 0:em.K
		for K2 = 0:em.K
			k1 = K1 * pi / em.L
			k2 = K2 * pi / em.L

			Sk = N*t*(em.ck[K1+1,K2+1] - em.phik[K1+1,K2+1])
			m = em.Gamma[K1+1,K2+1] * Sk

			bx += -k1 * m * sin(k1*x1) * cos(k2*x2) / em.h[K1+1, K2+1]
			by += -k2 * m * cos(k1*x1) * sin(k2*x2) / em.h[K1+1, K2+1]
		end
	end
	return bx,by
end

function update_phi!(em::ErgodicManager, b::Belief)
	n = em.L
	F = em.F
	for vx = 1:n
		for vy = 1:n
			ma = 0.0; mb = 0.0; mc = 0.0; md = 0.0
			for theta_x = 1:n
				for theta_y = 1:n
					p = b[theta_x, theta_y]
					ma += F[vx, vy, theta_x, theta_y,1,1] * p
					mb += F[vx, vy, theta_x, theta_y,1,2] * p
					mc += F[vx, vy, theta_x, theta_y,2,1] * p
					md += F[vx, vy, theta_x, theta_y,2,2] * p
				end
			end
			em.phi[vx, vy] = ma*md - mb*mc
		end
	end
end


text:
using Logging

const DEFAULT_COMPILE_THRESHOLD = 1.04

abstract type Optimizer end

mutable struct TopNOptimizer{TCostModel} <: Optimizer
    n::Int
    compile_threshold::Float32
    costmodel::TCostModel
    last_fixtypes
    fixtypes_history
    TopNOptimizer(
        n = 10;
        compile_threshold = DEFAULT_COMPILE_THRESHOLD,
        costmodel = basemodel) = new{typeof(costmodel)}(n, compile_threshold, costmodel, EmptyTypeList, [])
end

fixtypes(opt::TopNOptimizer, ::NoProfiler) = opt.last_fixtypes

percent(x) = round((x) * 10_000) / 100

function find_best_historic_cost(opt, _typefreqs, costmodel)
    best_historic_cost = typemax(ClockCycle)
    best_idx = -1
    for idx in 1:length(opt.fixtypes_history) # TODO Do not calculate for all if the list is lengthy
        old = opt.fixtypes_history[idx]
        oldcost = costof(_typefreqs, old, costmodel)
        if oldcost < best_historic_cost
            best_historic_cost = oldcost
            best_idx = idx
        end
    end
    return best_idx, best_historic_cost
end

function fixtypes(opt::TopNOptimizer, prof::FullProfiler)
    tfs = typefreqs(prof)
    @debug "Profiled type freqs: $(tfs)"
    ideal = ideal_fixtypes(opt, prof)
    idealcost = costof(tfs, ideal, opt.costmodel)
    best_idx, best_historic_cost = find_best_historic_cost(opt, tfs, opt.costmodel)
    if best_historic_cost > opt.compile_threshold * idealcost
        @debug begin
            lastcost = costof(typefreqs(prof), opt.last_fixtypes, opt.costmodel)
            speedup_last = percent(1.0 - idealcost / lastcost)
            speedup_best = best_historic_cost == typemax(Int) ? speedup_last : percent(1.0 - idealcost / best_historic_cost)
            "Selected new fixtype list for $(speedup_last)% speedup vs last used, and $(speedup_best)% vs best: $(Catwalk.pretty(ideal))"
        end
        opt.last_fixtypes = ideal
        push!(opt.fixtypes_history, ideal)
    else
        opt.last_fixtypes = opt.fixtypes_history[best_idx]
        @debug "Selected previous fixtype list: $(Catwalk.pretty(opt.last_fixtypes))"
    end
    return opt.last_fixtypes
end

# The N most frequent first types encoded
function ideal_fixtypes(opt::TopNOptimizer, prof::FullProfiler)
    freqs = sort(collect(typefreqs(prof)); order=Base.Order.ReverseOrdering())
    topfreqs = freqs[1:min(opt.n, length(freqs))]
    return encode(map(f -> f.type, topfreqs)...)
end


text:
module Example

using Glimmer, Glimmer.FlexUI
using Meshing
using FileIO # MeshIO should also be installed
using GeometryBasics
using StaticArrays
using Colors
using LinearAlgebra

println("Start [$(splitext(basename(@__FILE__))[1])]")

#---------------------------------------------------------------
# define the application and some basic properties such as title and initial window size
#---------------------------------------------------------------
app = App()
prop!(app, :title, "Glimmer Example - Meshing")
prop!(app, :winInitWidth, 1400)
prop!(app, :winInitHeight,1200)

#---------------------------------------------------------------
# create the MeshCat 3D viewer
#---------------------------------------------------------------
scene = Scene(openWindow = false)
set_Z_up!(scene)
grid!(scene, false)
axes!(scene, false)
cameraTransform!(scene, lookAt(SVector(-1.0, 1.0, 1.0)*10, zero3()))
cameraPlanes!(scene, 0.1, 1000.0)

#---------------------------------------------------------------
# Define Variables
#---------------------------------------------------------------
surface_code = addVariable!(app, Variable(name="surface_code", type="string",
        value="""begin
                    gyroid(v) = cos(v[1])*sin(v[2])+cos(v[2])*sin(v[3])+cos(v[3])*sin(v[1]);
                    gyroid_shell(v) = max(gyroid(v)-0.4,-gyroid(v)-0.4);
                end"""))
status = addVariable!(app, Variable(name="status", type="string", value = ""))

show_mesh =addVariable!(app, Variable(name="show_mesh", type="bool",value=true))
mesh_opacity = addVariable!(app, Variable(name="mesh_opacity", type="flota64",value=0.8))
mesh_wireframe =addVariable!(app, Variable(name="mesh_wireframe", type="bool",value=false))

show_pc =addVariable!(app, Variable(name="show_pc", type="bool",value=false))
pc_opacity = addVariable!(app, Variable(name="pc_opacity", type="flota64",value=0.5))

sample_resolution = addVariable!(app, Variable(name="sample_resolution", type="int",value=40))

#---------------------------------------------------------------
# Define Controls
#---------------------------------------------------------------
ui = VContainer(
    Card(
        title="Marching Cube",
        content=VContainer(
            Field(
                label = "Surface Definition",
                input = "multiline",
                style = "width: 100%; font-family: monospace; font-size: 20px;",
                hint = "defining the sureface to be constructed",
                variable="surface_code"
            ),
            Label(
                text = "\$(status)",
                class = "h2",
                style="color: red;"
            ),

            Card(
                # title="Controls",
                content=HContainer(
                    Card(
                        # title="Controls",
                        content=VContainer(
                            Card(
                                title="Mesh Controls",
                                content=VContainer(
                                    CheckBox(
                                        label="Show Mesh",
                                        variable="show_mesh"
                                    ),
                                    Slider(
                                        text="Mesh Opacity",
                                        trailingText="[\$()]",
                                        min=0,
                                        max=1,
                                        value=0.5,
                                        step=0.1,
                                        variable="mesh_opacity"
                                    ),
                                    CheckBox(
                                        label="Wireframe",
                                        variable="mesh_wireframe"
                                    ),
                                ),
                            ),

                            Card(
                                title="Points Controls",
                                content=VContainer(
                                    CheckBox(
                                        label="Show Points",
                                        variable="show_pc"
                                    ),
                                    Slider(
                                        text="Points Opacity",
                                        trailingText="[\$()]",
                                        min=0,
                                        max=1,
                                        value=0.5,
                                        step=0.1,
                                        variable="pc_opacity"
                                    ),
                                ),
                            ),

                            Card(
                                title="Sampling",
                                content=VContainer(
                                    Slider(
                                        text="Sample Resolution",
                                        trailingText="[\$()]",
                                        min=4,
                                        max=100,
                                        value=20,
                                        variable="sample_resolution"
                                    ),
                                ),
                            ),

                        ),
                    ),

                    Card(
                        title="Surface Visualization",
                        style="width: 50vw;",
                        content=VContainer(
                            MeshCatViewer(
                                url =  url(scene),
                                width = "100%",
                                height = "600px",
                            ),
                        ),
                    ),

                ),
            ),

        ),
    ),

    Glimmer.exampleSourceAsCard(@__FILE__),     # add the source code of the example as the last control
)
# set the controls for the application
controls!(app, ui)

#---------------------------------------------------------------
# the render function - preparing the image
#---------------------------------------------------------------
# gyroid(v) = cos(v[1])*sin(v[2])+cos(v[2])*sin(v[3])+cos(v[3])*sin(v[1])
# gyroid_shell(v) = max(gyroid(v)-0.4,-gyroid(v)-0.4)
function compile()
    try
        d = Meta.parse(surface_code[])
        eval(d)
        Base.invokelatest(gyroid_shell, [0.0, 0.0, 0.0])    # sanity check
        status[] = ""
    catch e
        status[] = "Encountered an error [$e]"
    end
end

function surface_function(v)
    Base.invokelatest(gyroid_shell, v)
end

# prepare a simple texture
img = zeros(RGB, 1, 4)
img[1,1] = RGB(1.0, 0.0, 0.0)
img[1,2] = RGB(0.0, 1.0, 0.0)
img[1,3] = RGB(0.0, 0.0, 1.0)
img[1,4] = RGB(1.0, 1.0, 0.0)

io = IOBuffer()
# save(io, img)
save(Stream(format"PNG", io), img)
res = String(take!(io))
tex = Glimmer.TextureData(Vector{UInt8}(res))

function render()
    # clear the scene
    clear(scene)

    # compile the UI code
    compile()

    # don't process if there is an error
    if (status[] != "")
        return;
    end

    # generate directly using GeometryBasics API
    # Rect specifies the sampling intervals
    gy_mesh = GeometryBasics.Mesh(
        surface_function,
        Rect(Vec(pi*4,pi*4,pi*4)*-0.5, Vec(pi*4,pi*4,pi*4)*1),
        MarchingCubes(),
        samples=(sample_resolution[], sample_resolution[], sample_resolution[])
    )

    # create the texture coordinates for the new mesh
    verts = coordinates(gy_mesh)
    faces = GeometryBasics.faces(gy_mesh)
    N = length(verts)
    point_attributes = Dict{Symbol, Any}()
    point_attributes[:uv] = Vector{SVector{2, Float32}}(undef, N)
    for (i, v) in enumerate(verts)
        x = norm(v)
        y = 0.1
        point_attributes[:uv][i] = SVector(x, y)
    end

    gy_mesh = GeometryBasics.Mesh(meta(verts; point_attributes...), faces)

    # origin axes
    origin = Axes(
        tr=identityTransform(),
        shaft_scale=0.1,
        axes_scale=5.0,
        name="Origin")
    parent!(origin , root(scene))

    # triangulated mesh
    if (show_mesh[])
        # mat = Material(color=RGBA(0.2, 0.8, 0.2, mesh_opacity[]), wireframe=mesh_wireframe[])
        mat = Material(color=RGBA(1.0, 1.0, 1.0, mesh_opacity[]), wireframe=mesh_wireframe[], texture=tex)
        mesh = Glimmer.Mesh(mesh=gy_mesh, material=mat, name="MC Mesh")
        parent!(mesh , root(scene))
    end

    # vertices
    if (show_pc[])
        points = [SVector{3, Float64}(p) for p in coordinates(gy_mesh)]
        pc_mat = Material(color=RGBA(0.7, 0.2, 0.1, pc_opacity[]), size=0.05)
        pc = PointCloud(points=points, tr=identityTransform(), material=pc_mat, name="Mesh Points")
        parent!(pc , root(scene))
    end
end
renderFunction!(app, render)

#---------------------------------------------------------------
# Run the application
#---------------------------------------------------------------
run(app)

println("End [$(splitext(basename(@__FILE__))[1])]")

GeometryBasics.normal_mesh

end # module

# # save("gyroid.ply", gy_mesh)


text:
using Git
using Test

import JLLWrappers

get_env(name) = get(ENV, name, nothing)
const orig_libpath     = deepcopy(get_env(JLLWrappers.LIBPATH_env))
const orig_execpath    = deepcopy(get_env("GIT_EXEC_PATH"))
const orig_cainfo      = deepcopy(get_env("GIT_SSL_CAINFO"))
const orig_templatedir = deepcopy(get_env("GIT_TEMPLATE_DIR"))

function withtempdir(f::Function)
    mktempdir() do tmp_dir
        cd(tmp_dir) do
            f(tmp_dir)
        end
    end
    return nothing
end

@testset "Git.jl" begin
    withtempdir() do tmp_dir
        @test !isdir("Git.jl")
        @test !isfile(joinpath("Git.jl", "Project.toml"))
        run(`$(git()) clone https://github.com/JuliaVersionControl/Git.jl`)
        @test isdir("Git.jl")
        @test isfile(joinpath("Git.jl", "Project.toml"))
    end

    withtempdir() do tmp_dir
        @test !isdir("Git.jl")
        @test !isfile(joinpath("Git.jl", "Project.toml"))
        run(`$(git()) clone https://github.com/JuliaVersionControl/Git.jl`)
        @test isdir("Git.jl")
        @test isfile(joinpath("Git.jl", "Project.toml"))
    end
end

@testset "Safety" begin
    # Make sure `git` commands don't leak environment variables
    @test orig_libpath == get_env(JLLWrappers.LIBPATH_env)
    @test orig_execpath == get_env("GIT_EXEC_PATH")
    @test orig_cainfo == get_env("GIT_SSL_CAINFO")
    @test orig_templatedir == get_env("GIT_TEMPLATE_DIR")
end


text:
__precompile__()

module JACKAudio

using Compat
import Compat.ASCIIString

using SampledSignals
import SampledSignals: nchannels, samplerate, nframes
using Base.Libc: malloc, free

export JACKClient, sources, sinks, seekavailable

# TODO: Logging is segfaulting when used inside the precompiled callback function
# using Logging

# Logging.configure(level=DEBUG)

include("libjack.jl")

# the ringbuffer size will be this times sizeof(float) rounded up to the nearest
# power of two
const RINGBUF_SAMPLES = 131072
# we'll advance the ringbuf read pointer by this many samples whenever we
# detect an overflow.
const OVERFLOW_ADVANCE = 8192

function __init__()
    global const process_cb = cfunction(process, Cint, (NFrames, Ptr{Ptr{Void}}))
    global const shutdown_cb = cfunction(shutdown, Void, (Ptr{JACKClient}, ))
    global const info_handler_cb = cfunction(info_handler, Void, (Cstring, ))
    global const error_handler_cb = cfunction(error_handler, Void, (Cstring, ))

    ccall((:jack_set_info_function, :libjack), Void, (Ptr{Void},),
        info_handler_cb)
    ccall((:jack_set_error_function, :libjack), Void, (Ptr{Void},),
        error_handler_cb)
end

function error_handler(msg)
    println("libjack: ERROR: $(unsafe_string(msg))")

    nothing
end

function info_handler(msg)
    println("libjack: $(unsafe_string(msg))")

    nothing
end

"""Represents a single-channel stream to or from JACK. Each port has a JACK
ringbuffer to get data in and out of the process callback and a julia RingBuffer
object that can be used for additional buffering. One important difference is
that if the jack ringbuffer fills up it can't be written to, but the Julia
RingBuffer maintains the last N samples"""
immutable JACKPort
    name::ASCIIString
    ptr::PortPtr
    jackbuf::RingBufPtr

    function JACKPort(client, name, porttype)
        ptr = jack_port_register(client, name, JACK_DEFAULT_AUDIO_TYPE, porttype, 0)
        if isnullptr(ptr)
            error("Failed to create port for $name")
        end

        bufptr = jack_ringbuffer_create(RINGBUF_SAMPLES * sizeof(JACKSample))
        if isnullptr(bufptr)
            jack_port_unregister(client, ptr)
            error("Failed to create ringbuffer for $name")
        end

        new(name, ptr, bufptr)
    end
end

# JACKSource and JACKSink defs are almost identical, so DRY it out with some
# metaprogramming magic
for (T, Super, porttype) in
        [(:JACKSource, :SampleSource, :PortIsInput),
         (:JACKSink, :SampleSink, :PortIsOutput)]
    """Represents a multi-channel stream, so it contains multiple jack ports
    that are considered synchronized, i.e. you read or write to all of them
    as a group. There can be multiple JACKSources and JACKSinks in a JACKClient,
    and all the sources and sinks in a client get updated by the same `process`
    method."""
    @eval immutable $T <: $Super
        name::ASCIIString
        client::ClientPtr
        clientname::ASCIIString
        ports::Vector{JACKPort}
        ringcondition::Condition # used to synchronize any in-progress transations

        function $T(client::ClientPtr, clientname, name, nchan)
            ports = JACKPort[]
            for ch in 1:nchan
                pname = portname(name, nchan, ch)
                push!(ports, JACKPort(client, pname, $porttype))
            end

            # TODO: switch to a mutable type, add finalizer, null out pointers
            # after closing/freeing them
            new(name, client, clientname, ports, Condition())
        end
    end

    @eval samplerate(stream::$T) = jack_get_sample_rate(stream.client)
    @eval nchannels(stream::$T) = length(stream.ports)
    @eval Base.eltype(stream::$T) = JACKSample

    @eval function Base.close(s::$T)
        while length(s.ports) > 0
            port = pop!(s.ports)
            jack_port_unregister(s.client, port.ptr)
            jack_ringbuffer_free(port.jackbuf)
        end
    end

    @eval function Base.show(io::IO, s::$T)
        print(io, $T, "(\"$(s.name)\", $(length(s.ports)))")
    end

    @eval isopen(s::$T) = length(s.ports) > 0
end

"""Generate the name of an individual port. This is what shows up in a JACK port
list"""
function portname(name, totalchans, chan)
    suffix = totalchans == 1 ? "" : "_$chan"
    string(name, suffix)
end

"""A `JACKClient` represents a connection to the JACK server. It can contain
multiple `JACKSource`s and `JACKSink`s. It is automatically activated when it is
constructed, so the sources and sinks are available for reading and writing,
respectively."""
type JACKClient
    name::ASCIIString
    ptr::ClientPtr
    sources::Vector{JACKSource}
    sinks::Vector{JACKSink}
    # this is memory allocated separately with malloc that is used to give the
    # process callback all the pointers it needs for the source/sink ports and
    # ringbuffers
    portptrs::Ptr{Ptr{Void}}
    callback::Base.SingleAsyncWork

    # this constructor takes a list of name, channelcount pairs
    function JACKClient{T1 <: Tuple, T2 <: Tuple}(
            name::AbstractString="Julia",
            sources::Vector{T1}=[("In", 2)],
            sinks::Vector{T2}=[("Out", 2)];
            connect=true, active=true)
        status = Ref{Cint}(Failure)
        clientptr = jack_client_open(name, 0, status)
        if isnullptr(clientptr)
            error("Failure opening JACK Client: ", status_str(status[]))
        end
        if status[] & ServerStarted
            # info("Started JACK Server")
        end
        if status[] & NameNotUnique
            name = unsafe_string(jack_get_client_name(clientptr))
            info("Given name not unique, renamed to ", name)
        end
        # println("Opened JACK Client with status: ", status_str(status[]))

        # we malloc 2*nsources + 2*nsinks + 3, because for each source and sink
        # we have the port pointer and the ringbuf pointer, the source and sink
        # lists are null-terminated, and we need to include the callback handle
        nsources = sum([p[2] for p in sources])
        nsinks = sum([p[2] for p in sinks])
        nptrs = 2nsources + 2nsinks + 3
        # TODO: we can switch this malloc and unsafe_store business
        # to an array we push! to
        portptrs = Ptr{Ptr{Void}}(malloc(nptrs*sizeof(Ptr{Void})))
        if isnullptr(portptrs)
            jack_client_close(clientptr)
            error("Failure allocating memory for JACK client \"$name\"")
        end

        # for now we leave the callback field uninitialized because we need the
        # client reference to build the callback closure
        client = new(name, clientptr, JACKSource[], JACKSink[], portptrs)

        # TODO: break out the source/sink addition to separate functions
        # initialize the sources and sinks
        ptridx = 1
        try
            for sourceargs in sources
                source = JACKSource(clientptr, name, sourceargs[1], sourceargs[2])
                push!(client.sources, source)
                # copy pointers to our flat pointer list that we'll give to the callback
                for port in source.ports
                    unsafe_store!(portptrs, port.ptr, ptridx)
                    unsafe_store!(portptrs, port.jackbuf, ptridx+1)
                    ptridx += 2
                end
            end
        catch
            close(client)
            rethrow()
        end
        # list of sources is null-terminated
        unsafe_store!(portptrs, C_NULL, ptridx)
        ptridx += 1

        try
            for sinkargs in sinks
                sink = JACKSink(clientptr, name, sinkargs[1], sinkargs[2])
                push!(client.sinks, sink)
                # copy pointers to our flat pointer list that we'll give to the callback
                for port in sink.ports
                    unsafe_store!(portptrs, port.ptr, ptridx)
                    unsafe_store!(portptrs, port.jackbuf, ptridx+1)
                    ptridx += 2
                end
            end
        catch
            close(client)
            rethrow()
        end
        # list of sinks is null-terminated
        unsafe_store!(portptrs, C_NULL, ptridx)
        ptridx += 1

        client.callback = Base.SingleAsyncWork(data -> managebuffers(client))

        # and finally we store the callback handle so the JACK process callback
        # can trigger the managebuffers function to run in the julia context
        unsafe_store!(portptrs, client.callback.handle, ptridx)

        jack_set_process_callback(clientptr, process_cb, portptrs)
        jack_on_shutdown(clientptr, shutdown_cb, pointer_from_objref(client))

        # useful when debugging, because you'll see errors. not sure how safe it
        # is though
        # process(128, portptrs)

        finalizer(client, close)

        if active
            activate(client)
            if connect
                autoconnect(client)
            end
        end

        client
    end
end

SampledSignals.samplerate(client::JACKClient) = jack_get_sample_rate(client.ptr)

# also allow constructing just by giving channel counts
JACKClient(name::AbstractString,
            sourcecount::Integer, sinkcount::Integer; kwargs...) =
    JACKClient(name, [("In", sourcecount)], [("Out", sinkcount)]; kwargs...)

JACKClient(sourcecount::Integer, sinkcount::Integer; kwargs...) =
    JACKClient("Julia", [("In", sourcecount)], [("Out", sinkcount)]; kwargs...)

function Base.show(io::IO, client::JACKClient)
    print(io, "JACKClient(\"$(client.name)\", [")
    sources = ASCIIString["(\"$(source.name)\", $(nchannels(source)))" for source in client.sources]
    sinks = ASCIIString["(\"$(sink.name)\", $(nchannels(sink)))" for sink in client.sinks]
    join(io, sources, ", ")
    print(io, "], [")
    join(io, sinks, ", ")
    print(io, "])")
end

function Base.close(client::JACKClient)
    if !isnullptr(client.ptr)
        deactivate(client)
    end
    while length(client.sources) > 0
        source = pop!(client.sources)
        close(source)
    end
    while length(client.sinks) > 0
        sink = pop!(client.sinks)
        close(sink)
    end
    closestatus = Cint(Success)
    if !isnullptr(client.ptr)
        closestatus = jack_client_close(client.ptr)
        client.ptr = C_NULL
    end
    if !isnullptr(client.portptrs)
        free(client.portptrs)
        client.portptrs = C_NULL
    end
    if closestatus != Cint(Success)
        error("Error closing client $(client.name): $(status_str(status))")
    end

    nothing
end

isopen(client::JACKClient) = !isnullptr(client.ptr)

sources(client::JACKClient) = client.sources
sinks(client::JACKClient) = client.sinks

# define wrapper methods to more conveniently access the sources and sinks
Base.read!(client::JACKClient, args...) = read!(client.sources[1], args...)
Base.read(client::JACKClient, args...) = read(client.sources[1], args...)
Base.write(client::JACKClient, args...) = write(client.sinks[1], args...)
Base.connect(c1::JACKClient, c2::JACKClient) = connect(c1.sinks[1], c2.sources[1])

# TODO: julia PR to extend Base.isnull rather than using isnullptr
isnullptr(ptr::Ptr) = Ptr{Void}(ptr) == C_NULL
isnullptr(ptr::Cstring) = Ptr{Cchar}(ptr) == C_NULL

"""Connect the given client to the physical input/output ports, by just matching
them up sequentially"""
function autoconnect(client::JACKClient)
    # look for physical output ports (the output from the sound card is an input
    # for us)
    ports = jack_get_ports(client.ptr, Ptr{Cchar}(C_NULL), Ptr{Cchar}(C_NULL),
        PortIsPhysical | PortIsOutput)
    if !isnullptr(ports)
        idx = 1
        for stream in client.sources
            for ch in 1:length(stream.ports)
                isnullptr(unsafe_load(ports, idx)) && break
                localportname = string(client.name, ":",
                                       portname(stream.name, length(stream.ports), ch))
                jack_connect(client.ptr, unsafe_load(ports, idx), pointer(localportname))

                idx += 1
            end
            isnullptr(unsafe_load(ports, idx)) && break
        end
        jack_free(ports)
    end
    ports = jack_get_ports(client.ptr, Ptr{Cchar}(C_NULL), Ptr{Cchar}(C_NULL),
        PortIsPhysical | PortIsInput)
    if !isnullptr(ports)
        idx = 1
        for stream in client.sinks
            for ch in 1:length(stream.ports)
                isnullptr(unsafe_load(ports, idx)) && break
                localportname = string(client.name, ":",
                                       portname(stream.name, length(stream.ports), ch))
                jack_connect(client.ptr, pointer(localportname), unsafe_load(ports, idx))

                idx += 1
            end
            isnullptr(unsafe_load(ports, idx)) && break
        end
        jack_free(ports)
    end
end

function Base.connect(sink::JACKSink, source::JACKSource)
    for (sinkport, sourceport) in zip(sink.ports, source.ports)
        sinkportname = string(sink.clientname, ":", sinkport.name)
        sourceportname = string(source.clientname, ":", sourceport.name)
        jack_connect(sink.client, sinkportname, sourceportname)
    end
end

# this should only get called during construction
function activate(client::JACKClient)
    status = ccall((:jack_activate, :libjack), Cint, (ClientPtr, ), client.ptr)
    if status != Int(Success)
        error("Error activating client $(client.name): $(status_str(status))")
    end

    nothing
end

# this should only get called when closing the client
function deactivate(client::JACKClient)
    status = ccall((:jack_deactivate, :libjack), Cint, (ClientPtr, ), client.ptr)
    if status != Int(Success)
        error("Error deactivating client $(client.name): $(status_str(status))")
    end

    nothing
end

# handle writes from a buffer with matching channel count and sample rate. Up/Down
# mixing and resampling should be the responsibility of SampledSignals.jl
function SampledSignals.unsafe_write(sink::JACKSink, buf::Array, frameoffset, framecount)
    isopen(sink) || return 0
    byteswritten = Csize_t(0)
    totalbytes = Csize_t(framecount * sizeof(JACKSample))
    chanptrs = Ptr{JACKSample}[channelptr(buf, c, frameoffset) for c in 1:nchannels(buf)]
    ports = sink.ports

    n = min(bytesavailable(sink), totalbytes)
    for ch in 1:length(ports)
        jack_ringbuffer_write(ports[ch].jackbuf, chanptrs[ch], n)
        chanptrs[ch] += n
    end
    byteswritten += n
    while byteswritten < totalbytes
        # wait to be notified that some space has freed up in the ringbuf
        wait(sink.ringcondition)
        isopen(sink) || return Int(div(byteswritten, sizeof(JACKSample)))
        n = min(bytesavailable(sink), totalbytes - byteswritten)
        for ch in 1:length(ports)
            jack_ringbuffer_write(ports[ch].jackbuf, chanptrs[ch], n)
            chanptrs[ch] += n
        end
        byteswritten += n
    end

    # by now we know we've written the whole length of the buffer
    framecount
end

function overflowed(source::JACKSource)
    for port in source.ports
        if jack_ringbuffer_write_space(port.jackbuf) < sizeof(JACKSample)
            return true
        end
    end

    false
end

"""advances the read pointer of the ring buffer without actually reading the
data"""
function ringbuf_read_advance(source::JACKSource, amount)
    for port in source.ports
        jack_ringbuffer_read_advance(port.jackbuf, amount)
    end
end

# read from the given source into the buffer, assuming that the channel count,
# sample rate and element type match
function SampledSignals.unsafe_read!(source::JACKSource, buf::Array, frameoffset, framecount)
    isopen(source) || return 0
    byteswritten = Csize_t(0)
    totalbytes = Csize_t(framecount * sizeof(JACKSample))
    chanptrs = Ptr{JACKSample}[channelptr(buf, ch, frameoffset) for ch in 1:nchannels(buf)]
    ports = source.ports

    # do the first read immediately
    nbytes = min(bytesavailable(source), totalbytes)
    for ch in 1:length(ports)
        jack_ringbuffer_read(ports[ch].jackbuf, chanptrs[ch], nbytes)
        chanptrs[ch] += nbytes
    end
    byteswritten += nbytes

    # now we wait to be notified that there's new data to read
    while byteswritten < totalbytes
        # wait to be notified that new samples are available in the ringbuf
        wait(source.ringcondition)
        isopen(source) || return Int(div(byteswritten, sizeof(JACKSample)))
        nbytes = min(bytesavailable(source), totalbytes - byteswritten)
        for ch in 1:length(ports)
            jack_ringbuffer_read(ports[ch].jackbuf, chanptrs[ch], nbytes)
            chanptrs[ch] += nbytes
        end
        byteswritten += nbytes
    end

    # by now we know we've read the whole length of the buffer
    framecount
end

function seekavailable(source::JACKSource)
    ringbuf_read_advance(source, bytesavailable(source))
end

"""Gives the least number of bytes available for writing among all the sink's
channels, aligned to a JACKSample boundry"""
function bytesavailable(sink::JACKSink)
    space = typemax(Csize_t)
    for port in sink.ports
        space = min(space, jack_ringbuffer_write_space(port.jackbuf))
    end

    sampalign(space)
end

"""Gives the number of frames available for writing to this sink"""
navailable(sink::JACKSink) = div(bytesavailable(sink), sizeof(JACKSample))

"""Gives the least number of bytes available for reading among all the source's
channels, aligned to a JACKSample boundry"""
function bytesavailable(source::JACKSource)
    space = typemax(Csize_t)
    for port in source.ports
        space = min(space, jack_ringbuffer_read_space(port.jackbuf))
    end

    sampalign(space)
end

"""Gives the number of frames available for reading from this source"""
navailable(source::JACKSource) = div(bytesavailable(source), sizeof(JACKSample))

# This gets called from a separate thread, so it is VERY IMPORTANT that it not
# allocate any memory or JIT compile when it's being run. Here be segfaults.
function process(nframes, portptrs)
    nbytes::Csize_t = nframes * sizeof(JACKSample)

    ptridx = 1
    # handle sources

    # we want to find the minimum amount of space in any ringbuffer, so we keep
    # all the channels synchronized even if a channel is overflowing
    minbytes = nbytes
    while !isnullptr(unsafe_load(portptrs, ptridx))
        ringbuf = unsafe_load(portptrs, ptridx + 1)
        ptridx += 2
        minbytes = min(minbytes, jack_ringbuffer_write_space(ringbuf))
    end
    # make sure we only write whole samples
    minbytes = sampalign(minbytes)

    # rewind back to the beginning for the actual writes
    ptridx = 1

    while !isnullptr(unsafe_load(portptrs, ptridx))
        source = unsafe_load(portptrs, ptridx)
        ringbuf = unsafe_load(portptrs, ptridx + 1)
        ptridx += 2

        buf = jack_port_get_buffer(source, nframes)
        jack_ringbuffer_write(ringbuf, buf, minbytes)
    end
    # skip over the null terminator
    ptridx += 1

    sinkidx = ptridx

    # handle sinks

    minbytes = nbytes
    # we want to find the minimum number of bytes available in any ringbuffer,
    # so we keep all the channels synchronized
    while !isnullptr(unsafe_load(portptrs, ptridx))
        ringbuf = unsafe_load(portptrs, ptridx + 1)
        ptridx += 2
        minbytes = min(minbytes, jack_ringbuffer_read_space(ringbuf))
    end
    # make sure we only read whole samples
    minbytes = sampalign(minbytes)

    # rewind back to the beginning of the sinks
    ptridx = sinkidx
    while !isnullptr(unsafe_load(portptrs, ptridx))
        sink = unsafe_load(portptrs, ptridx)
        ringbuf = unsafe_load(portptrs, ptridx + 1)
        ptridx += 2

        buf = jack_port_get_buffer(sink, nframes)
        # we know we're able to read at least minbytes
        jack_ringbuffer_read(ringbuf, buf, minbytes)
        if minbytes != nbytes
            memset(buf+minbytes, 0, nbytes - minbytes)
        end
    end
    # skip over the null terminator
    ptridx += 1

    handle = unsafe_load(portptrs, ptridx)

    # notify the managebuffers, which will get called with a reference to
    # this client
    ccall(:uv_async_send, Void, (Ptr{Void},), handle)

    Cint(0)
end

"""Returns the largest x <= value s.t. x has the given alignment (in bytes)"""
align{T<:Unsigned}(value::T, alignment::Integer) = value & ~(T(alignment-1))
sampalign(value) = align(value, sizeof(JACKSample))

# this callback gets called from within the Julia event loop, but is triggered
# by every `process` call. It bumps any tasks waiting to read or write
function managebuffers(client)
    # make sure the client is still open
    isopen(client) || return

    for source in client.sources
        # if we've overflowed, advance the read head
        if overflowed(source)
            advance = min(OVERFLOW_ADVANCE, bytesavailable(source))
            ringbuf_read_advance(source, advance)
        end
        notify(source.ringcondition)
    end

    for sink in client.sinks
        notify(sink.ringcondition)
    end
end

function shutdown(arg)
    nothing
end

memset(buf, val, count) = ccall(:memset, Ptr{Void},
    (Ptr{Void}, Cint, Csize_t),
    buf, 0, count)

end # module


text:
# Portions of the generators.jl package are derived from NetworkX, which
# is licensed under the BSD License as follows:
#
# > NetworkX Python library and are BSD Licensed as follows:
# >
# > Copyright (C) 2004-2012, NetworkX Developers
# > Aric Hagberg <hagberg@lanl.gov>
# > Dan Schult <dschult@colgate.edu>
# > Pieter Swart <swart@lanl.gov>
# > All rights reserved.
# >
# > Redistribution and use in source and binary forms, with or without
# > modification, are permitted provided that the following conditions are
# > met:
# >
# >   * Redistributions of source code must retain the above copyright
# >     notice, this list of conditions and the following disclaimer.
# >
# >   * Redistributions in binary form must reproduce the above
# >     copyright notice, this list of conditions and the following
# >     disclaimer in the documentation and/or other materials provided
# >     with the distribution.
# >
# >   * Neither the name of the NetworkX Developers nor the names of its
# >     contributors may be used to endorse or promote products derived
# >     from this software without specific prior written permission.
# >
# >
# > THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
# > "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
# > LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
# > A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
# > OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
# > SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
# > LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
# > DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
# > THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
# > (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
# > OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

function simple_complete_graph(n::Integer; is_directed=true)
    g = simple_graph(n, is_directed=is_directed)
    for i = 1:n, j=1:n
        if is_directed
            if i != j
                add_edge!(g,i,j)
            end
        else
            if i < j
                add_edge!(g,i,j)
            end
        end
    end
    return g
end

function simple_star_graph(n::Integer; is_directed=true)
    g = simple_graph(n; is_directed=is_directed)
    for i = 2:n
        add_edge!(g,1,i)
    end
    return g
end

function simple_path_graph(n::Integer; is_directed=true)
    g = simple_graph(n; is_directed=is_directed)
    for i = 2:n
        add_edge!(g, i-1, i)
    end
    return g
end

function simple_wheel_graph(n::Integer; is_directed=true)
    g = simple_star_graph(n; is_directed=is_directed)
    for i = 3:n
        add_edge!(g, i-1, i)
    end
    if n != 2
        add_edge!(g, n,2)
    end
    return g
end

function _make_simple_undirected_graph{T<:Integer}(n::T, edgelist::Vector{@compat Tuple{T,T}})
    g = simple_graph(n, is_directed=false)
    for (s,d) in edgelist
        add_edge!(g,s,d)
    end
    return g
end

simple_diamond_graph() =
    _make_simple_undirected_graph(4, [(1,2), (1,3), (2,3), (2,4), (3,4)])

simple_bull_graph() =
    _make_simple_undirected_graph(5, [(1,2), (1,3), (2,3), (2,4), (3,5)])

function simple_chvatal_graph()
    e = [
        (1, 2), (1, 5), (1, 7), (1, 10),
        (2, 3), (2, 6), (2, 8),
        (3, 4), (3, 7), (3, 9),
        (4, 5), (4, 8), (4, 10),
        (5, 6), (5, 9),
        (6, 11), (6, 12),
        (7, 11), (7, 12),
        (8, 9), (8, 12),
        (9, 11),
        (10, 11), (10, 12)
    ]
    return _make_simple_undirected_graph(12,e)
end

function simple_cubical_graph()
    e = [
        (1, 2), (1, 4), (1, 5),
        (2, 3), (2, 8),
        (3, 4), (3, 7),
        (4, 6), (5, 6), (5, 8),
        (6, 7),
        (7, 8)
    ]
    return _make_simple_undirected_graph(8,e)
end

function simple_desargues_graph()
    e = [
        (1, 2), (1, 6), (1, 20),
        (2, 3), (2, 17),
        (3, 4), (3, 12),
        (4, 5), (4, 15),
        (5, 6), (5, 10),
        (6, 7),
        (7, 8), (7, 16),
        (8, 9), (8, 19),
        (9, 10), (9, 14),
        (10, 11),
        (11, 12), (11, 20),
        (12, 13),
        (13, 14), (13, 18),
        (14, 15),
        (15, 16),
        (16, 17),
        (17, 18),
        (18, 19),
        (19, 20)
    ]
    return _make_simple_undirected_graph(20,e)
end

function simple_dodecahedral_graph()
    e = [
        (1, 2), (1, 11), (1, 20),
        (2, 3), (2, 9),
        (3, 4), (3, 7),
        (4, 5), (4, 20),
        (5, 6), (5, 18),
        (6, 7), (6, 16),
        (7, 8),
        (8, 9), (8, 15),
        (9, 10),
        (10, 11), (10, 14),
        (11, 12),
        (12, 13), (12, 19),
        (13, 14), (13, 17),
        (14, 15),
        (15, 16),
        (16, 17),
        (17, 18),
        (18, 19),
        (19, 20)
    ]
    return _make_simple_undirected_graph(20,e)
end

function simple_frucht_graph()
    e = [
        (1, 2), (1, 7), (1, 8),
        (2, 3), (2, 8),
        (3, 4), (3, 9),
        (4, 5), (4, 10),
        (5, 6), (5, 10),
        (6, 7), (6, 11),
        (7, 11),
        (8, 12),
        (9, 10), (9, 12),
        (11, 12)
    ]
    return _make_simple_undirected_graph(20,e)
end

function simple_heawood_graph()
    e = [
    (1, 2), (1, 6), (1, 14),
    (2, 3), (2, 11),
    (3, 4), (3, 8),
    (4, 5), (4, 13),
    (5, 6), (5, 10),
    (6, 7),
    (7, 8), (7, 12),
    (8, 9),
    (9, 10), (9, 14),
    (10, 11),
    (11, 12),
    (12, 13),
    (13, 14)
    ]
    return _make_simple_undirected_graph(14,e)
end

function simple_house_graph()
    e = [ (1, 2), (1, 3), (2, 4), (3, 4), (3, 5), (4, 5) ]
    return _make_simple_undirected_graph(5,e)
end

function simple_house_x_graph()
    g = simple_house_graph()
    add_edge!(g, 1, 4)
    add_edge!(g, 2, 3)
    return g
end

function simple_icosahedral_graph()
    e = [
        (1, 2), (1, 6), (1, 8), (1, 9), (1, 12),
        (2, 3), (2, 6), (2, 7), (2, 9),
        (3, 4), (3, 7), (3, 9), (3, 10),
        (4, 5), (4, 7), (4, 10), (4, 11),
        (5, 6), (5, 7), (5, 11), (5, 12),
        (6, 7), (6, 12),
        (8, 9), (8, 10), (8, 11), (8, 12),
        (9, 10),
        (10, 11), (11, 12)
    ]
    return _make_simple_undirected_graph(12, e)
end

function simple_krackhardt_kite_graph()
    e = [
        (1, 2), (1, 3), (1, 4), (1, 6),
        (2, 4), (2, 5), (2, 7),
        (3, 4), (3, 6),
        (4, 5), (4, 6), (4, 7),
        (5, 7),
        (6, 7), (6, 8),
        (7, 8),
        (8, 9),
        (9, 10)
    ]
    return _make_simple_undirected_graph(10,e)
end

function moebius_kantor_graph()
    e = [
        (1, 2), (1, 6), (1, 16),
        (2, 3), (2, 13),
        (3, 4), (3, 8),
        (4, 5), (4, 15),
        (5, 6), (5, 10),
        (6, 7),
        (7, 8), (7, 12),
        (8, 9),
        (9, 10), (9, 14),
        (10, 11),
        (11, 12), (11, 16),
        (12, 13),
        (13, 14),
        (14, 15),
        (15, 16)
    ]
    return _make_simple_undirected_graph(16,e)
end

function simple_octahedral_graph()
    e = [
        (1, 2), (1, 3), (1, 4), (1, 5),
        (2, 3), (2, 4), (2, 6),
        (3, 5), (3, 6),
        (4, 5), (4, 6),
        (5, 6)
    ]
    return _make_simple_undirected_graph(6,e)
end

function simple_pappus_graph()
    e = [
        (1, 2), (1, 6), (1, 18),
        (2, 3), (2, 9),
        (3, 4), (3, 14),
        (4, 5), (4, 11),
        (5, 6), (5, 16),
        (6, 7),
        (7, 8), (7, 12),
        (8, 9), (8, 15),
        (9, 10),
        (10, 11), (10, 17),
        (11, 12),
        (12, 13),
        (13, 14), (13, 18),
        (14, 15),
        (15, 16),
        (16, 17),
        (17, 18)
    ]
    return _make_simple_undirected_graph(18,e)
end

function simple_petersen_graph()
    e = [
        (1, 2), (1, 5), (1, 6),
        (2, 3), (2, 7),
        (3, 4), (3, 8),
        (4, 5), (4, 9),
        (5, 10),
        (6, 8), (6, 9),
        (7, 9), (7, 10),
        (8, 10)
    ]
    return _make_simple_undirected_graph(10,e)
end

function simple_sedgewick_maze_graph()
    e = [
        (1, 3),
        (1, 6), (1, 8),
        (2, 8),
        (3, 7),
        (4, 5), (4, 6),
        (5, 6), (5, 7), (5, 8)
    ]
    return _make_simple_undirected_graph(8,e)
end

simple_tetrahedral_graph() =
    _make_simple_undirected_graph(4, [(1, 2), (1, 3), (1, 4), (2, 3), (2, 4), (3, 4)])

function simple_truncated_cube_graph()
    e = [
        (1, 2), (1, 3), (1, 5),
        (2, 12), (2, 15),
        (3, 4), (3, 5),
        (4, 7), (4, 9),
        (5, 6),
        (6, 17), (6, 19),
        (7, 8), (7, 9),
        (8, 11), (8, 13),
        (9, 10),
        (10, 18), (10, 21),
        (11, 12), (11, 13),
        (12, 15), (13, 14),
        (14, 22), (14, 23),
        (15, 16),
        (16, 20), (16, 24),
        (17, 18), (17, 19),
        (18, 21),
        (19, 20),
        (20, 24),
        (21, 22),
        (22, 23),
        (23, 24)
    ]
    return _make_simple_undirected_graph(24,e)
end

function simple_truncated_tetrahedron_graph()
    e = [
        (1, 2),(1, 3),(1, 10),
        (2, 3),(2, 7),
        (3, 4),
        (4, 5),(4, 12),
        (5, 6),(5, 12),
        (6, 7),(6, 8),
        (7, 8),
        (8, 9),
        (9, 10),(9, 11),
        (10, 11),
        (11, 12)
    ]
    return _make_simple_undirected_graph(12,e)
end

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


text:
using Documenter, MXNet

makedocs(
  modules = MXNet,
  doctest = false
)

deploydocs(
  deps = Deps.pip("pygments", "mkdocs", "mkdocs-material", "python-markdown-math"),
  repo = "github.com/dmlc/MXNet.jl.git",
  julia = "0.5",
)


text:
using SpecialFunctions: digamma

function emstep(comp::Geometric, samples::AbstractVector{Int}, weights::AbstractWeights, zz::AbstractVector; kwargs...)
    Z = sum(i -> zz[i] * weights[i], 1:length(samples))
    mx = sum(i -> zz[i] * weights[i] * samples[i], 1:length(samples)) / Z

    Geometric(1 / (1 + mx))
end

function emstep(::Poisson, samples::AbstractVector{Int}, weights::AbstractWeights, zz::AbstractVector; kwargs...)
    Z = sum(i -> zz[i] * weights[i], 1:length(samples))
    mx = sum(i -> zz[i] * weights[i] * samples[i], 1:length(samples)) / Z

    Poisson(mx)
end

function emstep(::Bernoulli, samples::Union{AbstractVector{Int}, AbstractVector{Bool}}, weights::AbstractWeights, zz::AbstractVector; reg=1e-12, kwargs...)
    Z = sum(i -> zz[i] * weights[i], 1:length(samples))
    mx = sum(i -> zz[i] * weights[i] * samples[i], 1:length(samples)) / Z

    Bernoulli(mx)
end

function emstep(comp::Binomial, samples::AbstractVector{Int}, weights::AbstractWeights, zz::AbstractVector; kwargs...)
    Z = sum(i -> zz[i] * weights[i], 1:length(samples))
    mx = sum(i -> zz[i] * weights[i] * samples[i], 1:length(samples)) / Z

    Binomial(comp.n, mx / comp.n)
end

# This is implemented using the algorithm by Huang et al.
function emstep(comp::NegativeBinomial, samples::AbstractVector{Int}, weights::AbstractWeights, zz::AbstractVector;
                p_eps=1e-3, r_eps=1e-3, reg=1e-12, kwargs...)
    dig = digamma(comp.r + 1e-12)
    Z = sum(i -> zz[i] * weights[i], 1:length(samples))
    mδ = comp.r * sum(i -> zz[i] * weights[i] * (digamma(comp.r + samples[i] + reg) - dig), 1:length(samples)) / Z
    mx = sum(i -> zz[i] * weights[i] * samples[i], 1:length(samples)) / Z

    beta = 1 - 1 / (1 - comp.p + reg) - 1 / log(comp.p + reg)

    theta = beta * mδ / (mx - (1 - beta) * mδ)

    theta = clamp(theta, p_eps, 1 - p_eps)

    new_r = -mδ / log(theta)
    new_r = max(new_r, r_eps)

    new_p = theta

    NegativeBinomial(new_r, new_p)
end


text:
import Pkg
Pkg.instantiate() # downloads all dependencies for the current project

Pkg.add("JSON")
Pkg.add("HTTP")
Pkg.add("DataFrames")
Pkg.add("JSON3")
Pkg.add("JSONTables")

Pkg.add("Plots")
Pkg.add("StatsPlots")

using HTTP, JSON, Plots, StatsPlots, DataFrames, JSON3, JSONTables

# JSON.parsefile("COVID-19 cases worldwide.json")

f = open("COVID-19 cases worldwide.json","r")

dicttxt = read(f,String)

dict2=JSON.parse(dicttxt)

println(dict2)


text:
#####################
# AFFINE OPERATIONS #
#####################

+(x::RevisedAffineForm) = x
-(x::RevisedAffineForm) = RevisedAffineForm(-x.c, -x.ϵ, x.Δ)

for op in (:+, :-)
    @eval function $op(x::RevisedAffineForm{N}, y::RevisedAffineForm{N}) where {N}
        c = $op(Interval(x.c), Interval(y.c))
        ϵ = $op(Interval.(x.ϵ), Interval.(y.ϵ))
        Δ = x.Δ ⊕₊ y.Δ
        return RevisedAffineForm(c, ϵ, Δ)
    end

    @eval function $op(x::RevisedAffineForm, a::Real)
        c = $op(Interval(x.c), Interval(a))
        Δ = x.Δ ⊕₊ radius(c)
        return RevisedAffineForm(mid(c), x.ϵ, Δ)
    end
end

+(a::Real, x::RevisedAffineForm) = x + a
-(a::Real, x::RevisedAffineForm) = -x + a

#########################
# NON-AFFINE OPERATIONS #
#########################
for (op, op_up) in ((:*, :mul_round), (:/, :div_round))
    @eval function $op(x::RevisedAffineForm, a::Real)
        c = $op(Interval(x.c), Interval(a))
        ϵ = $op(Interval.(x.ϵ), Interval(a))
        Δ = $op_up(x.Δ, float(abs(a)), RoundUp)
        return RevisedAffineForm(c, ϵ, Δ)
    end
end

*(a::Real, x::RevisedAffineForm) = *(x, a)

"""
Standard multiplication, see eq. (2.2) of section 2.2 of
*Iwona Skalna, Parametric Interval Algebraic Systems, Springer, 2018*
"""
function *(x::RevisedAffineForm{N}, y::RevisedAffineForm{N}) where {N}
    c = Interval(x.c) * Interval(y.c) + 0.5 * dot(Interval.(x.ϵ), Interval.(y.ϵ))
    ϵ = Interval(x.c) * Interval.(y.ϵ) + Interval(y.c) * Interval.(x.ϵ)
    Δ =
        abs(x.c) ⊗₊ y.Δ ⊕₊ abs(y.c) ⊗₊ x.Δ ⊕₊
        (_sum_up(abs, x.ϵ) ⊕₊ x.Δ) ⊗₊ (_sum_up(abs, y.ϵ) ⊕₊ y.Δ) ⊖₊
        0.5 ⊗₋ _dot_down(abs, x.ϵ, y.ϵ)
    return RevisedAffineForm(c, ϵ, Δ)
end

###########
# HELPERS #
###########

_sum_up(a) = reduce(⊕₊, a)
_sum_up(f, a) = mapreduce(f, ⊕₊, a)
_dot_up(a, b) = reduce(⊕₊, aᵢ ⊗₊ bᵢ for (aᵢ, bᵢ) in zip(a, b))
_dot_up(f, a, b) = reduce(⊕₊, f(aᵢ) ⊗₊ f(bᵢ) for (aᵢ, bᵢ) in zip(a, b))

_sum_down(a) = reduce(⊕₋, a)
_sum_down(f, a) = mapreduce(f, ⊕₋, a)
_dot_down(a, b) = reduce(⊕₋, aᵢ ⊗₋ bᵢ for (aᵢ, bᵢ) in zip(a, b))
_dot_down(f, a, b) = reduce(⊕₋, f(aᵢ) ⊗₋ f(bᵢ) for (aᵢ, bᵢ) in zip(a, b))


text:
using DataTables
using Estuaries
using Feather
using DataUtils
import DataFrames

const nrows = 10

function makefeather(filename::String, data::DataTable)
    Feather.write(filename, convert(DataFrames.DataFrame, data))
end

# data = DataTable(A=rand(nrows), B=[randstring(5) for i ∈ 1:nrows])
# makefeather("testfile.feather", data)

E = Estuary(Feather.Source("testfile.feather"), Feather.Sink("testfile.feather"))


text:
include("lpc.jl")

pcm, fs = wavread("voiced.wav")
x = pcm[:, 1]
pv = plot_lpc_spectrum(x, fs; title="LP spectral envelope (voiced)")

pcm, fs = wavread("unvoiced.wav")
x = pcm[:, 1]
pu = plot_lpc_spectrum(x, fs; title="LP spectral envelope (unvoiced)")

p = [pv; pu]
PlotlyJS.relayout!(p; width=800)
PlotlyJS.savefig(p, "lp.html"; js=:embed)


text:
not_LoadError(x::LoadError) = x.error
not_LoadError(x) = x

macro test_error(testee, tester)
    quote
        ($(esc(tester)))(
            try
                $(esc(testee))
            catch err
               $not_LoadError(err)
            end)
    end
end


text:
module TAMode

using OrdinaryDiffEq
using StaticArrays
using LinearAlgebra
import LabelledArrays
using Turing
import Statistics
import Distances: sqeuclidean

include("types.jl")
include("reactCode.jl")
include("bothLigands.jl")
include("compModel.jl")

const solTol = 1.0e-6

function domainDef(u, p, t)
    return any(x -> x < -solTol, u)
end

function getAutocrine(params::Union{Rates{T}, comprates{T}, Lsrates{T}})::Vector{T} where {T <: Real}
    if params isa Rates
        N = 55
    elseif params isa comprates
        N = 27
    elseif params isa Lsrates
        N = 30
    end

    return vec(runTAMinit([1.0e6], params, zeros(T, N)))
end

function runTAMinit(tps::AbstractVector{Float64}, params::Union{Rates{T}, comprates{T}, Lsrates{T}}, solInit::Vector{T})::Matrix{T} where {T <: Real}
    @assert all(x -> x >= 0.0, tps)

    prob = ODEProblem(TAMreact, solInit, maximum(tps), params)

    sol = AutoTsit5(Rodas4P(autodiff = (T == Float64)), stiffalgfirst = true)
    solut = solve(prob, sol; saveat = tps, reltol = solTol, isoutofdomain = domainDef, maxiters = 1e6)

    if solut.retcode != :Success
        println("Solving failed. Returning NaN.")
        println(solut.retcode)
        return fill(NaN, length(tps), length(solInit))
    else
        solut = solut.u
    end

    if length(tps) > 1
        solut = vcat(transpose.(solut)...)
    else
        solut = reshape(solut[1], (1, length(solInit)))
    end

    return solut
end

function compTAM(tps::AbstractVector{Float64}, params::Union{comprates{T}, Vector{T}}) where {T <: Real}
    if params isa Vector
        params = compParamm(params)
    end

    solInit = getAutocrine(params)
    u0 = repeat(solInit; outer = [compSize])

    sol = runTAMinit(tps, params, u0)

    return reshape(sol, (length(tps), 27, compSize))
end

function runTAM(tps::AbstractVector{Float64}, params::Union{Rates{T}, Lsrates{T}, Vector{T}}, ligStim)::Matrix{T} where {T <: Real}
    params = deepcopy(params) # Make sure we're not mutating the input

    if params isa Vector
        if length(params) == 15
            params = param(params)
        else
            @assert length(params) == 9
            params = Lsparam(params)
        end
    end

    solInit = getAutocrine(params)

    if params isa Rates
        params.gasCur += ligStim
    else
        @assert params isa Lsrates
        params.curL .+= ligStim
    end

    return runTAMinit(tps, params, solInit)
end

include("fitting.jl")
include("plots.jl")

end # module


text:
text = (open(f->read(f, String), "crud.txt"))

println(text)

max_stars_repo_name: codinginbrazil/functional-programming

max_issues_repo_name: codinginbrazil/functional-programming

max_forks_repo_name: codinginbrazil/functional-programming


text:
# retrieve the parameter to change
collision_cost = float(ARGS[1])
@eval @everywhere collision_cost = $collision_cost

@everywhere begin
    using AutoViz, Reel, POMDPs, StatsBase, POMDPToolbox, SARSOP, QMDP, AutomotiveDrivingModels, JLD, Parameters, Distributions, GridInterpolations

    include("../../admenv/crosswalk/occluded_crosswalk_env.jl")
    include("../../admenv/crosswalk/pomdp_types.jl")
    include("../../admenv/crosswalk/spaces.jl")
    include("../../admenv/crosswalk/transition.jl")
    include("../../admenv/crosswalk/observation.jl")
    include("../../admenv/crosswalk/belief.jl")
    include("../../admenv/crosswalk/adm_helpers.jl")
    include("../../admenv/crosswalk/decomposition.jl")

    include("config.jl")
    include("ego_control.jl")
    include("sensor.jl")
    include("policy.jl")
    include("state_estimation.jl")
    include("constant_pedestrian.jl")
    include("simulation.jl")
    include("baseline_policy.jl")
    include("render_helpers.jl")

    function evaluate(env::CrosswalkEnv, ego_model::CrosswalkDriver, config::EvalConfig)
        crash = 0
        avg_step = 0
        success = 0
        time_out = 0
        avg_step = 0
        ego0 = initial_ego(env, config.rng)
        reset_model!(ego_model, ego0)
        models = Dict{Int, DriverModel}()
        models[1] = ego_model
        scene = initial_scene(models, env, config)
        push!(scene, ego0)
        nticks = config.time_out
        rec = SceneRecord(nticks+1, config.sim_dt)

        # execute the simulation
        simulate!(rec, scene, env, models, nticks, config.rng, config.n_ped, config.callbacks)
        avg_step = nframes(rec)

        # analyse output
        if is_crash(rec[0])
            crash = 1
            #if !isfile("crash.gif")
            #    overlay = QMDPOverlay(ego_model)
            #    reset_overlay!(overlay, ego0)
            #    duration, fps, render_rec = animate_record(rec, overlay)
            #    speed_factor = 1
            #    film = roll(render_rec, fps = fps*speed_factor, duration = duration/speed_factor)
            #    write("crash.gif", film)
            #end
            println("Episode CRASH")
        elseif is_terminal(rec[0], env)
            success = 1
            println("Episode SUCCESS")
        elseif nframes(rec) >= config.time_out
            time_out = 1
            println("Episode  TIME OUT")
        else
            println("Episode nothing happened, frames: $(nframes(rec))")
        end

        println("End of the simulation: Crashes $crash , Successes $success, Time out $time_out, Avg Step $(avg_step)")
        return crash, success, time_out, avg_step
    end

    config = EvalConfig(time_out = 200,
                        n_ped=2,
                        n_episodes = 1000)

    params = EnvParams(ped_rate = 0.01)
    env = CrosswalkEnv(params)
    cam = FitToContentCamera(0.0)
    pomdp = OCPOMDP()
    pomdp.collision_cost = collision_cost
    sarsop_policy = load_policy(pomdp, "/scratch/boutonm/crosswalk$(abs(pomdp.collision_cost)).policy")

    sensor = POMDPSensor(pomdp = pomdp, sensor = SimpleSensor(0.5, 0.5))
    policy = QMDPEval(env, pomdp, sarsop_policy);
    up = MixedUpdater(pomdp, config.sim_dt);
    a0 = 0.
    o0 = Dict{Int64, OCObs}()
    ego0 = initial_ego(env, config.rng)
    b0 = Dict{Int64, OCDistribution}()
    b0[-1] = initialstate_distribution(policy.pomdp, ego0.state)
    update_freq = 5
    ego_model = CrosswalkDriver(update_freq, 0, env, a0, b0, o0, ego0, sensor, policy, up)
    reset_model!(ego_model, ego0)

end
envs = fill(env, config.n_episodes)
ego_models = fill(ego_model, config.n_episodes)
configs = Vector{EvalConfig}(config.n_episodes)
for i=1:config.n_episodes
    configs[i] = EvalConfig(time_out = config.time_out, n_ped = config.n_ped, n_episodes = config.n_episodes,
                            rng = MersenneTwister(i))
end
evals = pmap(evaluate, envs, ego_models, configs)

crashes = zeros(config.n_episodes)
successes = zeros(config.n_episodes)
time_outs = zeros(config.n_episodes)
steps = zeros(config.n_episodes)
for (i,t) in enumerate(evals)
    crashes[i], successes[i], time_outs[i], steps[i] = t
end

c_rate = mean(crashes)*100
suc_rate = mean(successes)*100
to_rate = mean(time_outs)*100
avg_time = mean(steps)*config.sim_dt

println("Outcome of $(config.n_episodes) episodes : \n",
         "$c_rate \% crashes \n ",
         "$suc_rate\% successes \n ",
         "$to_rate\% time out \n ",
         "$(avg_time)(s) avg time \n ")

scenario = "crosswalk-min"
goal_reward = pomdp.goal_reward
collision_cost = pomdp.collision_cost
p_birth = pomdp.p_birth
n_ped = config.n_ped

result_file = "../results_sarsop.csv"
if !isfile(result_file)
    f = open(result_file, "w")
    write(f, "scenario, goal_reward, action_cost, collision_cost, p_birth, n_ped, n_episodes, crashes, successes, time_out, avg_time \n")
else
    f = open(result_file, "a")
end
res_str = string(scenario, ", ", pomdp.goal_reward, ", ", pomdp.action_cost, ", ", pomdp.collision_cost, ", ", pomdp.p_birth, ", ",
                 config.n_ped, ", ", config.n_episodes, ", ",
                 c_rate, ", ", suc_rate, ", ", to_rate, ", ", avg_time, "\n")
write(f, res_str)
close(f)

using DataFrames
convergence_file="conv_sarsop_$p_birth$collision_cost.csv"
convergence_data = DataFrame([crashes, successes, time_outs, steps], [:crashes, :successes, :time_outs, :avg_steps])
writetable(convergence_file, convergence_data)

max_stars_repo_name: MaximeBouton/AutomotivePOMDPs

max_issues_repo_name: MaximeBouton/AutomotivePOMDPs

max_forks_repo_name: MaximeBouton/AutomotivePOMDPs