text:
module Triangle
export basic_triangulation
export basic_triangulation_vertices
export constrained_triangulation
export constrained_triangulation_vertices

if !isfile(joinpath(dirname(@__FILE__), "..", "deps", "deps.jl"))
  error("Triangle is not properly installed. Please run\nPkg.build(\"Triangle\")")
end

include("common.jl")

"""
    basic_triangulation_vertices(vertices::Array{Float64,2})

Compute a Delaunay triangulation for a list of `vertices` in the form of `[x1 y1; x2 y2; ... ; xn yn]`

The function will return an array of array of 3-vertices lists (triangles with the correct vertices order) using the vertex coordinates in each triangle definition.

# Example
```jldoctest
julia> using Triangle

julia> points = Array{Float64,2}([0. 0.; 1. 0.; 0. 1.])
3×2 Array{Float64,2}:
 0.0  0.0
 1.0  0.0
 0.0  1.0

julia> Triangle.basic_triangulation_vertices(points)
1-element Array{Array{Float64,2},1}:
 [0.0 0.0; 1.0 0.0; 0.0 1.0]
```
"""
function basic_triangulation_vertices(vertices::Array{Float64,2})
  basic_triangulation_vertices(vertices, Array{Int64,1}(collect(1:1:size(vertices)[1])))
end

"""
    basic_triangulation(vertices::Array{Float64,2},vertices_map::Array{Int64,1})

Compute a Delaunay triangulation for a list of `vertices` in the form of `[x1 y1; x2 y2; ... ; xn yn]`

A list of indexes is provided in `vertices_map` so that each vertex can have a custom integer identifier.

The function will return an array of array of 3-vertices lists (triangles with the correct vertices order) using the vertex identifiers in each triangle definition.

# Example
```jldoctest
julia> using Triangle

julia> points = Array{Float64,2}([0. 0.; 1. 0.; 0. 1.])
3×2 Array{Float64,2}:
 0.0  0.0
 1.0  0.0
 0.0  1.0

julia> points_map = [1, 2, 3]
3-element Array{Int64,1}:
 1
 2
 3

julia> Triangle.basic_triangulation(points,points_map)
1-element Array{Array{Int64,1},1}:
 [1, 2, 3]
```
"""
function basic_triangulation(vertices::Array{Float64,2}, vertices_map::Array{Int64,1})
  Triangulate.basic_triangulation(vertices, vertices_map)
end

"""
    basic_triangulation_vertices(vertices::Array{Float64,2},vertices_map::Array{Int64,1})

Compute a Delaunay triangulation for a list of `vertices` in the form of `[x1 y1; x2 y2; ... ; xn yn]`

A list of indexes is provided in `vertices_map` so that each vertex can have a custom integer identifier.

The function will return an array of array of 3-vertices lists (triangles with the correct vertices order) using the vertex coordinates in each triangle definition.

# Example
```jldoctest
julia> using Triangle

julia> points = Array{Float64,2}([0. 0.; 1. 0.; 0. 1.])
3×2 Array{Float64,2}:
 0.0  0.0
 1.0  0.0
 0.0  1.0

julia> points_map = [1, 2, 3]
3-element Array{Int64,1}:
 1
 2
 3

julia> Triangle.basic_triangulation_vertices(points,points_map)
1-element Array{Array{Float64,2},1}:
 [0.0 0.0; 1.0 0.0; 0.0 1.0]
```
"""
function basic_triangulation_vertices(vertices::Array{Float64,2}, vertices_map::Array{Int64,1})
  map_to_vertices( vertices, vertices_map, basic_triangulation(vertices, vertices_map) )
end

"""
    constrained_triangulation(vertices::Array{Float64,2}, vertices_map::Array{Int64,1}, edges_list::Array{Int64,2})

Compute a Constrained Delaunay triangulation for a list of `vertices` in the form of `[x1 y1; x2 y2; ... ; xn yn]` and a list of edges that will be kept.

A list of indexes is provided in `vertices_map` so that each vertex can have a custom integer identifier.

A list of edges (to be included in the final triangulation) is passed in `edges_list` in the form of `[ vertex-identifier-1 vertex-identifier-2; vertex-identifier-1 vertex-identifier-3; ... ; vertex-identifier-N vertex-identifier-M ]`

The function will return an array of array of 3-vertices lists (triangles with the correct vertices order) using the vertex identifiers in each triangle definition.

# Example
```jldoctest
julia> using Triangle

julia> points = [0. 0.; 0. 3.; 1. 3.; 1. 1.; 2. 1.; 2. 0.]
6×2 Array{Float64,2}:
 0.0  0.0
 0.0  3.0
 1.0  3.0
 1.0  1.0
 2.0  1.0
 2.0  0.0

julia> points_map = Array{Int64,1}(collect(1:1:size(points)[1]))
6-element Array{Int64,1}:
 1
 2
 3
 4
 5
 6

julia> edges_list = Array{Int64,2}([1 2; 2 3; 3 4; 4 5; 5 6; 6 1])
6×2 Array{Int64,2}:
 1  2
 2  3
 3  4
 4  5
 5  6
 6  1

julia> Triangle.constrained_triangulation(points,points_map,edges_list)
4-element Array{Array{Int64,1},1}:
 [1, 4, 2]
 [4, 1, 6]
 [2, 4, 3]
 [5, 4, 6]
```
"""
function constrained_triangulation(vertices::Array{Float64,2}, vertices_map::Array{Int64,1}, edges_list::Array{Int64,2})
  Triangulate.constrained_triangulation(vertices, vertices_map, edges_list, fill(true, size(edges_list)[1])  )
end

"""
    constrained_triangulation_vertices(vertices::Array{Float64,2}, vertices_map::Array{Int64,1}, edges_list::Array{Int64,2})

Compute a Constrained Delaunay triangulation for a list of `vertices` in the form of `[x1 y1; x2 y2; ... ; xn yn]` and a list of edges that will be kept.

A list of indexes is provided in `vertices_map` so that each vertex can have a custom integer identifier.

A list of edges (to be included in the final triangulation) is passed in `edges_list` in the form of `[ vertex-identifier-1 vertex-identifier-2; vertex-identifier-1 vertex-identifier-3; ... ; vertex-identifier-N vertex-identifier-M ]`

The function will return an array of array of 3-vertices lists (triangles with the correct vertices order) using the vertex coordinates in each triangle definition.

# Example
```jldoctest
julia> using Triangle

julia> points = [0. 0.; 0. 3.; 1. 3.; 1. 1.; 2. 1.; 2. 0.]
6×2 Array{Float64,2}:
 0.0  0.0
 0.0  3.0
 1.0  3.0
 1.0  1.0
 2.0  1.0
 2.0  0.0

julia> points_map = Array{Int64,1}(collect(1:1:size(points)[1]))
6-element Array{Int64,1}:
 1
 2
 3
 4
 5
 6

julia> edges_list = Array{Int64,2}([1 2; 2 3; 3 4; 4 5; 5 6; 6 1])
6×2 Array{Int64,2}:
 1  2
 2  3
 3  4
 4  5
 5  6
 6  1

julia> Triangle.constrained_triangulation_vertices(points,points_map,edges_list)
4-element Array{Array{Float64,2},1}:
 [0.0 0.0; 1.0 1.0; 0.0 3.0]
 [1.0 1.0; 0.0 0.0; 2.0 0.0]
 [0.0 3.0; 1.0 1.0; 1.0 3.0]
 [2.0 1.0; 1.0 1.0; 2.0 0.0]
```
"""
function constrained_triangulation_vertices(vertices::Array{Float64,2}, vertices_map::Array{Int64,1}, edges_list::Array{Int64,2})
  map_to_vertices( vertices, vertices_map, constrained_triangulation(vertices, vertices_map, edges_list, fill(true, size(edges_list)[1])) )
end

"""
    constrained_triangulation(vertices::Array{Float64,2}, vertices_map::Array{Int64,1}, edges_list::Array{Int64,2}, edges_boundary::Array{Bool,1})

Compute a Constrained Delaunay triangulation for a list of `vertices` in the form of `[x1 y1; x2 y2; ... ; xn yn]` and a list of edges that will be kept. Some of those edges can be marked as the boundary of the mesh.

A list of indexes is provided in `vertices_map` so that each vertex can have a custom integer identifier.

A list of edge (to be included in the final triangulation) is passed in `edges_list` in the form of `[ vertex-identifier-1 vertex-identifier-2; vertex-identifier-1 vertex-identifier-3; ... ; vertex-identifier-N vertex-identifier-M ]`

A list of boundary markers passed in `edges_boundary` in the form of booleans that tell the triangulator if the edge is on the boundary or not (the indexing is the same of `edges_list`).

The function will return an array of array of 3-vertices lists (triangles with the correct vertices order) using the vertex identifiers in each triangle definition.

# Example
```jldoctest
julia> using Triangle

julia> points = [0. 0.; 0. 3.; 1. 3.; 1. 1.; 2. 1.; 2. 0.]
6×2 Array{Float64,2}:
 0.0  0.0
 0.0  3.0
 1.0  3.0
 1.0  1.0
 2.0  1.0
 2.0  0.0

julia> points_map = Array{Int64,1}(collect(1:1:size(points)[1]))
6-element Array{Int64,1}:
 1
 2
 3
 4
 5
 6

julia> edges_list = Array{Int64,2}([1 2; 2 3; 3 4; 4 5; 5 6; 6 1])
6×2 Array{Int64,2}:
 1  2
 2  3
 3  4
 4  5
 5  6
 6  1

julia> edge_boundary = [false, false, true, true, false, false]
6-element Array{Bool,1}:
 false
 false
  true
  true
 false
 false

julia> Triangle.constrained_triangulation(points,points_map,edges_list,edge_boundary)
4-element Array{Array{Int64,1},1}:
 [1, 4, 2]
 [4, 1, 6]
 [2, 4, 3]
 [5, 4, 6]
```
"""
function constrained_triangulation(vertices::Array{Float64,2}, vertices_map::Array{Int64,1}, edges_list::Array{Int64,2}, edges_boundary::Array{Bool,1})
  Triangulate.constrained_triangulation(vertices, vertices_map, edges_list, edges_boundary)
end

"""
    constrained_triangulation_vertices(vertices::Array{Float64,2}, vertices_map::Array{Int64,1}, edges_list::Array{Int64,2}, edges_boundary::Array{Bool,1})

Compute a Constrained Delaunay triangulation for a list of `vertices` in the form of `[x1 y1; x2 y2; ... ; xn yn]` and a list of edges that will be kept. Some of those edges can be marked as the boundary of the mesh.

A list of indexes is provided in `vertices_map` so that each vertex can have a custom integer identifier.

A list of edge (to be included in the final triangulation) is passed in `edges_list` in the form of `[ vertex-identifier-1 vertex-identifier-2; vertex-identifier-1 vertex-identifier-3; ... ; vertex-identifier-N vertex-identifier-M ]`

A list of boundary markers passed in `edges_boundary` in the form of booleans that tell the triangulator if the edge is on the boundary or not (the indexing is the same of `edges_list`).

The function will return an array of array of 3-vertices lists (triangles with the correct vertices order) using the vertex coordinates in each triangle definition.

# Example
```jldoctest
julia> using Triangle

julia> points = [0. 0.; 0. 3.; 1. 3.; 1. 1.; 2. 1.; 2. 0.]
6×2 Array{Float64,2}:
 0.0  0.0
 0.0  3.0
 1.0  3.0
 1.0  1.0
 2.0  1.0
 2.0  0.0

julia> points_map = Array{Int64,1}(collect(1:1:size(points)[1]))
6-element Array{Int64,1}:
 1
 2
 3
 4
 5
 6

julia> edges_list = Array{Int64,2}([1 2; 2 3; 3 4; 4 5; 5 6; 6 1])
6×2 Array{Int64,2}:
 1  2
 2  3
 3  4
 4  5
 5  6
 6  1

julia> edge_boundary = [false, false, true, true, false, false]
6-element Array{Bool,1}:
 false
 false
  true
  true
 false
 false

julia> Triangle.constrained_triangulation_vertices(points,points_map,edges_list,edge_boundary)
4-element Array{Array{Float64,2},1}:
 [0.0 0.0; 1.0 1.0; 0.0 3.0]
 [1.0 1.0; 0.0 0.0; 2.0 0.0]
 [0.0 3.0; 1.0 1.0; 1.0 3.0]
 [2.0 1.0; 1.0 1.0; 2.0 0.0]
```
"""
function constrained_triangulation_vertices(vertices::Array{Float64,2}, vertices_map::Array{Int64,1}, edges_list::Array{Int64,2}, edges_boundary::Array{Bool,1})
  map_to_vertices( vertices, vertices_map, constrained_triangulation(vertices, vertices_map, edges_list, edges_boundary) )
end

"""
    constrained_triangulation(vertices::Array{Float64,2}, vertices_map::Array{Int64,1}, edges_list::Array{Int64,2}, edges_boundary::Array{Bool,1}, holes::Array{Float64,2})

Compute a Constrained Delaunay triangulation for a list of `vertices` in the form of `[x1 y1; x2 y2; ... ; xn yn]` and a list of edges that will be kept. Some of those edges can be marked as the boundary of the mesh.

A list of `holes` in the form of _[x1 y1; x2 y2; ... ; xn yn]_ can be passed to avoid triangulation in that part of the mesh.

A list of indexes is provided in `vertices_map` so that each vertex can have a custom integer identifier.

A list of edge (to be included in the final triangulation) is passed in `edges_list` in the form of `[ vertex-identifier-1 vertex-identifier-2; vertex-identifier-1 vertex-identifier-3; ... ; vertex-identifier-N vertex-identifier-M ]`

A list of boundary markers passed in `edges_boundary` in the form of booleans that tell the triangulator if the edge is on the boundary or not (the indexing is the same of `edges_list`).

The function will return an array of array of 3-vertices lists (triangles with the correct vertices order) using the vertex identifiers in each triangle definition.

# Example
```jldoctest
julia> using Triangle

julia> points = Array{Float64,2}([0. 0.; 4. 0.; 2. 3.; 8. 0.; 6. 3.; 4. 6.])
6×2 Array{Float64,2}:
 0.0  0.0
 4.0  0.0
 2.0  3.0
 8.0  0.0
 6.0  3.0
 4.0  6.0

julia> points_map = [1, 2, 4, 3, 5, 6]
6-element Array{Int64,1}:
 1
 2
 4
 3
 5
 6

julia> edges_list = Array{Int64,2}([1 2; 2 3; 3 1; 2 4; 4 5; 5 2; 3 5; 5 6; 6 3])
9×2 Array{Int64,2}:
 1  2
 2  3
 3  1
 2  4
 4  5
 5  2
 3  5
 5  6
 6  3

julia> edge_boundary = [false,true,false,false,false,true,true,false,false]
9-element Array{Bool,1}:
 false
  true
 false
 false
 false
  true
  true
 false
 false

julia> holes = [4. 2.]
1×2 Array{Float64,2}:
 4.0  2.0

julia> Triangle.constrained_triangulation(points, points_map, edges_list, edge_boundary, holes)
3-element Array{Array{Int64,1},1}:
 [1,2,3]
 [5,6,3]
 [4,5,2]
```
"""
function constrained_triangulation(vertices::Array{Float64,2}, vertices_map::Array{Int64,1}, edges_list::Array{Int64,2}, edges_boundary::Array{Bool,1}, holes::Array{Float64,2})
  Triangulate.constrained_triangulation(vertices, vertices_map, edges_list, edges_boundary, holes)
end

"""
    constrained_triangulation_vertices(vertices::Array{Float64,2}, vertices_map::Array{Int64,1}, edges_list::Array{Int64,2}, edges_boundary::Array{Bool,1}, holes::Array{Float64,2})

Compute a Constrained Delaunay triangulation for a list of `vertices` in the form of `[x1 y1; x2 y2; ... ; xn yn]` and a list of edges that will be kept. Some of those edges can be marked as the boundary of the mesh.

A list of `holes` in the form of _[x1 y1; x2 y2; ... ; xn yn]_ can be passed to avoid triangulation in that part of the mesh.

A list of indexes is provided in `vertices_map` so that each vertex can have a custom integer identifier.

A list of edge (to be included in the final triangulation) is passed in `edges_list` in the form of `[ vertex-identifier-1 vertex-identifier-2; vertex-identifier-1 vertex-identifier-3; ... ; vertex-identifier-N vertex-identifier-M ]`

A list of boundary markers passed in `edges_boundary` in the form of booleans that tell the triangulator if the edge is on the boundary or not (the indexing is the same of `edges_list`).

The function will return an array of array of 3-vertices lists (triangles with the correct vertices order) using the vertex coordinates in each triangle definition.

# Example
```jldoctest
julia> using Triangle

julia> points = Array{Float64,2}([0. 0.; 4. 0.; 2. 3.; 8. 0.; 6. 3.; 4. 6.])
6×2 Array{Float64,2}:
 0.0  0.0
 4.0  0.0
 2.0  3.0
 8.0  0.0
 6.0  3.0
 4.0  6.0

julia> points_map = [1, 2, 4, 3, 5, 6]
6-element Array{Int64,1}:
 1
 2
 4
 3
 5
 6

julia> edges_list = Array{Int64,2}([1 2; 2 3; 3 1; 2 4; 4 5; 5 2; 3 5; 5 6; 6 3])
9×2 Array{Int64,2}:
 1  2
 2  3
 3  1
 2  4
 4  5
 5  2
 3  5
 5  6
 6  3

julia> edge_boundary = [false,true,false,false,false,true,true,false,false]
9-element Array{Bool,1}:
 false
  true
 false
 false
 false
  true
  true
 false
 false

julia> holes = [4. 2.]
1×2 Array{Float64,2}:
 4.0  2.0

julia> Triangle.constrained_triangulation_vertices(points, points_map, edges_list, edge_boundary, holes)
3-element Array{Array{Float64,2},1}:
 [0.0 0.0; 4.0 0.0; 2.0 3.0]
 [6.0 3.0; 4.0 6.0; 2.0 3.0]
 [8.0 0.0; 6.0 3.0; 4.0 0.0]
```
"""
function constrained_triangulation_vertices(vertices::Array{Float64,2}, vertices_map::Array{Int64,1}, edges_list::Array{Int64,2}, edges_boundary::Array{Bool,1}, holes::Array{Float64,2})
  map_to_vertices( vertices, vertices_map, constrained_triangulation(vertices, vertices_map, edges_list, edges_boundary, holes) )
end

function map_to_vertices(vertices::Array{Float64,2}, vertices_map::Array{Int64,1}, triangle_list::Array{Array{Int64, 1},1})
  vert_size = size(vertices)

  map_positions = Array{Int64,1}(undef, length(vertices_map))
  for i in 1:1:length(vertices_map)
    map_positions[vertices_map[i]] = i
  end

  triangle_list_vertices = Array{Array{Float64,2},1}()
  for i in 1:1:length(triangle_list)
        triangle = Array{Float64}(undef, 3,2)
        triangle[1] = vertices[map_positions[triangle_list[i][1]]]
        triangle[4] = vertices[map_positions[triangle_list[i][1]]+vert_size[1]]
        triangle[2] = vertices[map_positions[triangle_list[i][2]]]
        triangle[5] = vertices[map_positions[triangle_list[i][2]]+vert_size[1]]
        triangle[3] = vertices[map_positions[triangle_list[i][3]]]
        triangle[6] = vertices[map_positions[triangle_list[i][3]]+vert_size[1]]
        push!(triangle_list_vertices, triangle)
  end

  return triangle_list_vertices
end

end


text:
# Davis, Yin. "A Three-Operator Splitting Scheme and its Optimization
# Applications", Set-Valued and Variational Analysis, vol. 25, no. 4,
# pp. 829–858 (2017).

using Base.Iterators
using ProximalAlgorithms.IterationTools
using ProximalOperators: Zero
using LinearAlgebra
using Printf

"""
    DavisYinIteration(; <keyword-arguments>)

Instantiate the Davis-Yin splitting algorithm (see [1]) for solving
convex optimization problems of the form

    minimize f(x) + g(x) + h(x),

where `h` is smooth.

# Arguments
- `x0`: initial point.
- `f=Zero()`: proximable objective term.
- `g=Zero()`: proximable objective term.
- `h=Zero()`: smooth objective term.
- `Lh=nothing`: Lipschitz constant of the gradient of h.
- `gamma=nothing`: stepsize to use, defaults to `1/Lh` if not set (but `Lh` is).

# References
- [1] Davis, Yin. "A Three-Operator Splitting Scheme and its Optimization
Applications", Set-Valued and Variational Analysis, vol. 25, no. 4,
pp. 829–858 (2017).
"""

Base.@kwdef struct DavisYinIteration{R,C<:Union{R,Complex{R}},T<:AbstractArray{C},Tf,Tg,Th}
    f::Tf = Zero()
    g::Tg = Zero()
    h::Th = Zero()
    x0::T
    lambda::R = real(eltype(x0))(1)
    Lh::Maybe{R} = nothing
    gamma::Maybe{R} = Lh !== nothing ? (1 / Lh) : error("You must specify either Lh or gamma")
end

Base.IteratorSize(::Type{<:DavisYinIteration}) = Base.IsInfinite()

struct DavisYinState{T}
    z::T
    xg::T
    grad_h_xg::T
    z_half::T
    xf::T
    res::T
end

function Base.iterate(iter::DavisYinIteration)
    z = copy(iter.x0)
    xg, = prox(iter.g, z, iter.gamma)
    grad_h_xg, = gradient(iter.h, xg)
    z_half = 2 .* xg .- z .- iter.gamma .* grad_h_xg
    xf, = prox(iter.f, z_half, iter.gamma)
    res = xf - xg
    z .+= iter.lambda .* res
    state = DavisYinState(z, xg, grad_h_xg, z_half, xf, res)
    return state, state
end

function Base.iterate(iter::DavisYinIteration, state::DavisYinState)
    prox!(state.xg, iter.g, state.z, iter.gamma)
    gradient!(state.grad_h_xg, iter.h, state.xg)
    state.z_half .= 2 .* state.xg .- state.z .- iter.gamma .* state.grad_h_xg
    prox!(state.xf, iter.f, state.z_half, iter.gamma)
    state.res .= state.xf .- state.xg
    state.z .+= iter.lambda .* state.res
    return state, state
end

# Solver

struct DavisYin{R, K}
    maxit::Int
    tol::R
    verbose::Bool
    freq::Int
    kwargs::K
end

function (solver::DavisYin)(x0; kwargs...)
    stop(state::DavisYinState) = norm(state.res, Inf) <= solver.tol
    disp((it, state)) = @printf("%5d | %.3e\n", it, norm(state.res, Inf))
    iter = DavisYinIteration(; x0=x0, solver.kwargs..., kwargs...)
    iter = take(halt(iter, stop), solver.maxit)
    iter = enumerate(iter)
    if solver.verbose
        iter = tee(sample(iter, solver.freq), disp)
    end
    num_iters, state_final = loop(iter)
    return state_final.xf, state_final.xg, num_iters
end

DavisYin(; maxit=10_000, tol=1e-8, verbose=false, freq=100, kwargs...) =
    DavisYin(maxit, tol, verbose, freq, kwargs)

ext: jl


text:
using AbstractPlotting, GLMakie, MakieLayout, VideoIO, Dates, Observables

file = "a.mp4"
f = VideoIO.openvideo(file)

# some stats about the video
pixelaspectratio = VideoIO.aspect_ratio(f)
h = f.height
w = round(typeof(h), f.width*pixelaspectratio)
duration = VideoIO.get_duration(unsafe_load(f.avin.apFormatContext[1]))
fr = 1/f.framerate

# prepare observables
buff = read(f)
img = Node(rotr90(buff))
ct = Node("00:00:00.000")
cc = Node("-")

@enum MovieState begin
    PAUSE
    PLAY
    STEP
end

const FRAME_CHANNEL = Channel{typeof(rotr90(buff))}(1) # could buffer frames?
const STATE_CHANNEL = Channel{MovieState}(1)

push!(STATE_CHANNEL, PAUSE)

function movie_controller()
    state = PAUSE
    while true
        if isready(STATE_CHANNEL)
            state = take!(STATE_CHANNEL)
            while state == PAUSE
                state = take!(STATE_CHANNEL)
            end
        end
        frame = take!(FRAME_CHANNEL)
        img[] = frame
        if state == PLAY
            sleep(fr)
        elseif state == STEP
            push!(STATE_CHANNEL, PAUSE)
        else
            @warn state
        end
    end
end

function frame_producer()
    while true
        read!(f, buff)
        v = rotr90(buff)
        push!(FRAME_CHANNEL, v)
    end
end

@async movie_controller()
@async frame_producer()

h_control = 0.1 # relative height of the controlers
scene, layout = layoutscene(0, resolution = (Float64(w), h/(1 - h_control)), backgroundcolor = :yellow)
ax = layout[1, 1:5] = LAxis(scene, aspect = DataAspect(), backgroundcolor = :yellow)
layout[2, 1] = LText(scene, ct, textsize = 30, color = :white)
cl = layout[2, 2] = LText(scene, cc, textsize = 30, color = :white)
sl = layout[2, 3] = LSlider(scene, range = range(0, stop = duration, length = min(1000, round(Int, duration/fr))), startvalue = 0)

# slider seeks
sl2 = async_latest(sl.value, 1)
# sl2 = sl.value
on(sl2) do t
    take!(FRAME_CHANNEL)
    seek(f, t)
    push!(STATE_CHANNEL, STEP)
    # stepone!(f, buff, img)
end

# button stepone one frame forward
bt = layout[2,4] = LButton(scene, label = "⊵")
on(bt.clicks) do _
    # stepone!(f, buff, img)
    push!(STATE_CHANNEL, STEP)
end

# the toggle button for the playing
pl = layout[2, 5] = LToggle(scene)
on(pl.active) do p
    push!(STATE_CHANNEL, p ? PLAY : PAUSE)
end

imgplot = image!(ax, 1:w, 1:h, img)

# just layout stuff
rowsize!(layout, 1, Relative(1 - h_control))
rowsize!(layout, 2, Relative(h_control))
colsize!(layout, 1, Relative(0.2))
colsize!(layout, 2, Relative(0.2))
hidedecorations!(ax)
tightlimits!(ax)

# show current time stamp
on(img) do _
    t = VideoIO.gettime(f)
    ct[] = string(Time(0) + Millisecond(round(Int, 1000t)))
end

# some mouse interactions
mousestate = addmousestate!(ax.scene, imgplot)
onmouseover(mousestate) do state
    coordinate = round.(Int, state.pos)
    cc[] = string(coordinate)
end
onmouseout(mousestate) do state
    cc[] = "-"
end
onmouseleftclick(mousestate) do state
    @async begin
        cl.color[] = :red
        sleep(2)
        cl.color[] = :white
    end
end

scene


text:
tests = [
  "utils",
  "solver"
]

for t in tests
	f = "$t.jl"
	println("* running $f ...")
    t = @elapsed include(f)
    println("done (took $t seconds).")
end


text:
function AbstractPlotting.convert_arguments(::AbstractPlotting.PointBased, spec::KmerFrequencySpectra{1})
    tv = _find_plottable_subset(spec)
    return ([Point2f0(i - 1, j) for (i, j) in enumerate(tv)],)
end

function AbstractPlotting.convert_arguments(p::AbstractPlotting.SurfaceLike, spec::KmerFrequencySpectra{2})
    return AbstractPlotting.convert_arguments(p, _find_plottable_subset(spec))
end

AbstractPlotting.plottype(::KmerFrequencySpectra{1}) = AbstractPlotting.BarPlot
AbstractPlotting.plottype(::KmerFrequencySpectra{2}) = AbstractPlotting.Heatmap


text:
# Modified from Graphs.jl.
# Workaround for edge_index not being O(1) on SimpleGraph.
# edge_index was only needed to test for cycles, so
# this implementation skips that check.

# Graphs.jl is licensed under the MIT License:
#
# Copyright (c) 2012: John Myles White and other contributors.
#
# Permission is hereby granted, free of charge, to any person obtaining
# a copy of this software and associated documentation files (the
# "Software"), to deal in the Software without restriction, including
# without limitation the rights to use, copy, modify, merge, publish,
# distribute, sublicense, and/or sell copies of the Software, and to
# permit persons to whom the Software is furnished to do so, subject to
# the following conditions:
#
# The above copyright notice and this permission notice shall be
# included in all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
# EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
# NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
# LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
# OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
# WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

mutable struct TopologicalSortVisitor
    vertices::Vector{Int}
    parents::Vector{Int}

    function TopologicalSortVisitor(n::Int)
        vs = Int[]
        sizehint!(vs, n)
        new(vs, zeros(Int,n))
    end
end

function depth_first_visit_impl!(
    adjlist,
    offsets,
    vertex_stack,                   # an (initialized) stack of vertex
    index_stack,                    # stack with out edge indices
    vertexcolormap::Vector{Int},    # an (initialized) color-map to indicate status of vertices
    visitor::TopologicalSortVisitor)  # the visitor

    while !isempty(vertex_stack)
        u = pop!(vertex_stack)
        out_idx = pop!(index_stack)
        #uegs = out_edges(u,graph)
        len_uegs = offsets[u+1]-offsets[u]
        found_new_vertex = false

        while out_idx <= len_uegs && !found_new_vertex
            v = adjlist[offsets[u]+out_idx-1]
            out_idx += 1
            v_color = vertexcolormap[v]

            if v_color == 0
                found_new_vertex = true
                vertexcolormap[v] = 1
                push!(vertex_stack, u)
                push!(index_stack, out_idx)
                visitor.parents[v] = u

                push!(vertex_stack, v)
                push!(index_stack, 1)
            end
        end

        if !found_new_vertex
            close_vertex!(visitor, u)
            vertexcolormap[u] = 2
        end
    end
end

function traverse_graph(
    adjlist,
    offsets,
    s,
    visitor::TopologicalSortVisitor,
    vertexcolormap,
    vertex_stack,
    index_stack)

    vertexcolormap[s] = 1

    resize!(vertex_stack,1)
    vertex_stack[1] = s
    resize!(index_stack,1)
    index_stack[1] = 1

    depth_first_visit_impl!(adjlist, offsets, vertex_stack, index_stack, vertexcolormap, visitor)
end

function close_vertex!(visitor::TopologicalSortVisitor, v::Int)
    push!(visitor.vertices, v)
end

function reverse_topological_sort_by_dfs(adjlist,offsets,num_vertices, cmap::Vector{Int}, vertex_stack::Vector{Int}=Int[], index_stack::Vector{Int}=Int[])

    @assert length(cmap) == num_vertices
    fill!(cmap,0)
    visitor = TopologicalSortVisitor(num_vertices)

    for s in 1:num_vertices
        if cmap[s] == 0
            traverse_graph(adjlist, offsets, s, visitor, cmap, vertex_stack, index_stack)
        end
    end

    visitor.vertices,visitor.parents
end


text:
# p30.jl - spectral integration, ODE style (compare p12.jl)

# Computation: various values of N, four functions:
Nmax = 50; E = zeros(4,Nmax); clf();
for N = 1:Nmax
    i = 1:N; (D,x) = cheb(N);
    x = x[i]; Di = inv(D[i,i]); w = Di[1,:];
    f = @. abs(x)^3;     E[1,N] = abs(dot(w,f) - .5);
    f = @. exp(-x^(-2)); E[2,N] = abs(dot(w,f) - 2*(exp(-1)+sqrt(pi)*(erf(1)-1)));
    f = @. 1/(1+x^2);    E[3,N] = abs(dot(w,f) - pi/2);
    f = x.^10;           E[4,N] = abs(dot(w,f) - 2/11);
end

# Plot results:
labels = [L"|x|^3",L"\exp(-x^2)",L"1/(1+x^2)",L"x^{10}"];
for iplot = 1:4
    subplot(2,2,iplot)
    semilogy(E[iplot,:]+1e-100,".-",markersize=10);
    axis([0,Nmax,1e-18,1e3]); grid(true);
    xticks(0:10:Nmax); yticks((10.0).^(-15:5:0));
    ylabel("error"); text(32,.004,labels[iplot]);
end


text:
@testset "report no method matching" begin
    # if there is no method matching case, it should be reported
    let
        # NOTE: we can't just wrap them into `let`, closures can't be inferred correctly
        m = gen_virtual_module()
        analyzer, = Core.eval(m, quote
            foo(a::Integer) = :Integer
            $report_call((AbstractString,)) do a
                foo(a)
            end
        end)
        @test length(get_reports(analyzer)) === 1
        report = first(get_reports(analyzer))
        @test report isa NoMethodErrorReport
        @test report.t === Tuple{typeof(m.foo), AbstractString}
    end

    let
        analyzer, = report_call(()->sum([]))
        @test length(get_reports(analyzer)) === 1
        report = first(get_reports(analyzer))
        @test report isa SeriousExceptionReport
        @test report.err isa MethodError
        @test report.err.f === zero
    end

    # if there is no method matching case in union-split, it should be reported
    let
        m = gen_virtual_module()
        analyzer, = Core.eval(m, quote
            foo(a::Integer) = :Integer
            foo(a::AbstractString) = "AbstractString"

            $report_call(a->foo(a), (Union{Nothing,Int},))
        end)

        @test length(get_reports(analyzer)) === 1
        report = first(get_reports(analyzer))
        @test report isa NoMethodErrorReport
        @test report.t == [Tuple{typeof(m.foo), Nothing}]
    end
end

@testset "report local undefined variables" begin
    let
        analyzer, = report_call((Bool,)) do b
            if b
                bar = rand(Int)
                return bar
            end
            return bar # undefined in this pass
        end
        @test length(get_reports(analyzer)) === 1
        r = first(get_reports(analyzer))
        @test r isa LocalUndefVarErrorReport
        @test r.name === :bar
        @test last(r.vst).line == (@__LINE__)-6
    end

    # deeper level
    let
        m = @fixturedef begin
            function foo(b)
                if b
                    bar = rand(Int)
                    return bar
                end
                return bar # undefined in this pass
            end
            baz(a) = foo(a)
        end

        analyzer, = Core.eval(m, :($report_call(baz, (Bool,))))
        @test length(get_reports(analyzer)) === 1
        r = first(get_reports(analyzer))
        @test r isa LocalUndefVarErrorReport
        @test r.name === :bar
        @test last(r.vst).line == (@__LINE__)-10

        # works when cached
        analyzer, = Core.eval(m, :($report_call(baz, (Bool,))))
        @test length(get_reports(analyzer)) === 1
        r = first(get_reports(analyzer))
        @test r isa LocalUndefVarErrorReport
        @test r.name === :bar
        @test last(r.vst).line == (@__LINE__)-18
    end

    # try to exclude false negatives as possible (by collecting reports in after-optimization pass)
    let
        analyzer, = report_call((Bool,)) do b
            if b
                bar = rand()
            end

            return if b
                return bar # this shouldn't be reported
            else
                return nothing
            end
        end
        @test isempty(get_reports(analyzer))
    end

    let
        analyzer, = report_call((Bool,)) do b
            if b
                bar = rand()
            end

            return if b
                return nothing
            else
                # ideally we want to have report for this pass, but tons of work will be
                # needed to report this pass
                return bar
            end
        end
        @test_broken length(get_reports(analyzer)) === 1 &&
            first(get_reports(analyzer)) isa LocalUndefVarErrorReport &&
            first(get_reports(analyzer)).name === :bar
    end

    let
        analyzer, = report_call((Int,)) do a
            function inner(n)
                if n > 0
                   a = n
                end
            end
            inner(rand(Int))
            return a
        end
        @test isempty(get_reports(analyzer))
    end

    # with the current approach, local undefined variables in toplevel frame can't be found
    # since we don't cache toplevel frame and thus it won't be optimized
    let
        res = @analyze_toplevel begin
            foo = let
                if rand(Bool)
                    bar = rand(Int)
                else
                    bar # undefined in this pass
                end
            end
        end
        @test_broken !isempty(res.inference_error_reports)
    end
end

@testset "report undefined (global) variables" begin
    let
        analyzer, = report_call(()->foo)
        @test length(get_reports(analyzer)) === 1
        @test first(get_reports(analyzer)) isa GlobalUndefVarErrorReport
        @test first(get_reports(analyzer)).name === :foo
    end

    # deeper level
    let
        m = @fixturedef begin
            foo(bar) = bar + baz
            qux(a) = foo(a)
        end

        analyzer, = Core.eval(m, :($report_call(qux, (Int,))))
        @test length(get_reports(analyzer)) === 1
        @test first(get_reports(analyzer)) isa GlobalUndefVarErrorReport
        @test first(get_reports(analyzer)).name === :baz

        # works when cached
        analyzer, = Core.eval(m, :($report_call(qux, (Int,))))
        @test length(get_reports(analyzer)) === 1
        @test first(get_reports(analyzer)) isa GlobalUndefVarErrorReport
        @test first(get_reports(analyzer)).name === :baz
    end
end

@testset "report non-boolean condition error" begin
    # simple case
    let
        analyzer, = report_call((Int,)) do a
            a ? a : nothing
        end
        @test length(get_reports(analyzer)) === 1
        er = first(get_reports(analyzer))
        @test er isa NonBooleanCondErrorReport
        @test er.t === Int
    end

    # don't report when a type can be `Bool` (Bool ⊑ type)
    let
        analyzer, = report_call((Integer,)) do a
            a ? a : nothing
        end
        @test isempty(get_reports(analyzer))
    end

    # report union split case
    let
        analyzer, = report_call((Union{Nothing,Bool},)) do a
            a ? a : false
        end
        @test length(get_reports(analyzer)) === 1
        let r = first(get_reports(analyzer))
            @test r isa NonBooleanCondErrorReport
            @test r.t == [Nothing]
            @test occursin("for 1 of union split cases", r.msg)
        end
    end

    let
        analyzer, = report_call() do
            anyary = Any[1,2,3]
            first(anyary) ? first(anyary) : nothing
        end
        @test isempty(get_reports(analyzer)) # very untyped, we can't report on this ...
    end

    # `strict_condition_check` configuration
    let
        # `==(::Missing, ::Any)`
        analyzer, = report_call((Any,Symbol); strict_condition_check = false) do a, b
            a == b ? 0 : 1
        end
        @test isempty(get_reports(analyzer))
        analyzer, = report_call((Any,Symbol); strict_condition_check = true) do a, b
            a == b ? 0 : 1
        end
        @test any(get_reports(analyzer)) do r
            isa(r, NonBooleanCondErrorReport) &&
            r.t == Type[Missing]
        end
    end
end

@testset "inference with abstract global variable" begin
    let
        vmod = gen_virtual_module()
        res = @analyze_toplevel context = vmod virtualize = false begin
            s = "julia"
            sum(s)
        end

        @test is_concrete(vmod, :s)
        test_sum_over_string(res)
    end

    @testset "union assignment" begin
        let
            vmod = gen_virtual_module()
            res = @analyze_toplevel context = vmod virtualize = false begin
                global globalvar
                if rand(Bool)
                    globalvar = "String"
                else
                    globalvar = :Symbol
                end
            end

            @test is_abstract(vmod, :globalvar)
            @test vmod.globalvar.t === Union{String,Symbol}
        end

        let
            vmod = gen_virtual_module()
            res = @analyze_toplevel context = vmod virtualize = false begin
                if rand(Bool)
                    globalvar = "String"
                else
                    globalvar = :Symbol
                end

                foo(s::AbstractString) = length(s)
                foo(globalvar) # union-split no method matching error should be reported
            end

            @test is_abstract(vmod, :globalvar)
            @test vmod.globalvar.t === Union{String,Symbol}
            @test length(res.inference_error_reports) === 1
            er = first(res.inference_error_reports)
            @test er isa NoMethodErrorReport
            @test isa(er.t, Vector) # should be true
        end

        # sequential
        let
            vmod = gen_virtual_module()
            res = @analyze_toplevel context = vmod virtualize = false begin
                if rand(Bool)
                    globalvar = "String"
                else
                    globalvar = :Symbol
                end

                foo(s::AbstractString) = length(s)
                foo(globalvar) # union-split no method matching error should be reported

                globalvar = 10
                foo(globalvar) # no method matching error should be reported
            end

            @test is_concrete(vmod, :globalvar)
            @test is_analyzed(vmod, :globalvar)
            @test isa_analyzed(vmod.globalvar, Int)
            @test length(res.inference_error_reports) === 2
            let er = first(res.inference_error_reports)
                @test er isa NoMethodErrorReport
                @test isa(er.t, Vector)
            end
            let er = last(res.inference_error_reports)
                @test er isa NoMethodErrorReport
                @test !isa(er.t, Vector)
            end
        end
    end
end

@testset "UndefKeywordError" begin
    isa2(t) = x -> isa(x, t)
    let
        m = gen_virtual_module()
        analyzer, = Core.eval(m, quote
            foo(a; #= can be undef =# kw) =  a, kw
            $report_call(foo, (Any,))
        end)
        @test !isempty(get_reports(analyzer))
        @test any(get_reports(analyzer)) do r
            r isa SeriousExceptionReport || return false
            err = r.err
            err isa UndefKeywordError && err.var === :kw
        end
        # there shouldn't be duplicated report for the `throw` call
        @test !any(isa2(UncaughtExceptionReport), get_reports(analyzer))
    end
end

@testset "DivideError" begin
    isa2(t) = x -> isa(x, t)
    let
        apply(f, args...) = f(args...)

        analyzer, = report_call() do
            apply(div, 1, 0)
        end
        @test !isempty(get_reports(analyzer))
        @test any(isa2(DivideErrorReport), get_reports(analyzer))

        analyzer, = report_call() do
            apply(rem, 1, 0)
        end
        @test !isempty(get_reports(analyzer))
        @test any(isa2(DivideErrorReport), get_reports(analyzer))

        # JET analysis isn't sound
        analyzer, = report_call((Int,Int)) do a, b
            apply(div, a, b)
        end
        @test isempty(get_reports(analyzer))
    end
end

@testset "report `throw` calls" begin
    # simplest case
    let
        analyzer, = report_call(()->throw("foo"))
        @test !isempty(get_reports(analyzer))
        @test first(get_reports(analyzer)) isa UncaughtExceptionReport
    end

    # throws in deep level
    let
        foo(a) = throw(a)
        analyzer, = report_call(()->foo("foo"))
        @test !isempty(get_reports(analyzer))
        @test first(get_reports(analyzer)) isa UncaughtExceptionReport
    end

    # don't report possibly false negative `throw`s
    let
        foo(a) = a ≤ 0 ? throw("a is $(a)") : a
        analyzer, = report_call(foo, (Int,))
        @test isempty(get_reports(analyzer))
    end

    # constant prop sometimes helps exclude false negatives
    let
        foo(a) = a ≤ 0 ? throw("a is $(a)") : a
        analyzer, = report_call(()->foo(0))
        @test !isempty(get_reports(analyzer))
        @test first(get_reports(analyzer)) isa UncaughtExceptionReport
    end

    # report even if there're other "critical" error exist
    let
        m = gen_virtual_module()
        analyzer, = Core.eval(m, quote
            foo(a) = sum(a) # should be reported
            bar(a) = throw(a) # shouldn't be reported first
            $report_call((Bool, String)) do b, s
                b && foo(s)
                bar(s)
            end
        end)
        @test length(get_reports(analyzer)) === 3
        test_sum_over_string(get_reports(analyzer))
    end

    # end to end
    let
        # this should report `throw(ArgumentError("Sampler for this object is not defined")`
        analyzer, = report_call(rand, (Char,))
        @test !isempty(get_reports(analyzer))
        @test first(get_reports(analyzer)) isa UncaughtExceptionReport

        # this should not report `throw(DomainError(x, "sin(x) is only defined for finite x."))`
        analyzer, = report_call(sin, (Int,))
        @test isempty(get_reports(analyzer))

        # again, constant prop sometimes can exclude false negatives
        analyzer, = report_call(()->sin(Inf))
        @test !isempty(get_reports(analyzer))
        @test first(get_reports(analyzer)) isa UncaughtExceptionReport
    end
end

@testset "constant analysis" begin
    # constant prop should limit false positive union-split no method reports
    let
        m = @fixturedef begin
            mutable struct P
                i::Int
                s::String
            end
            foo(p, i) = p.i = i
        end

        # "for one of the union split cases, no matching method found for signature: Base.convert(Base.fieldtype(Base.typeof(x::P)::Type{P}, f::Symbol)::Union{Type{Int64}, Type{String}}, v::Int64)" should be threw away
        analyzer, = Core.eval(m, :($report_call(foo, (P, Int))))
        @test isempty(get_reports(analyzer))

        # works for cache
        analyzer, = Core.eval(m, :($report_call(foo, (P, Int))))
        @test isempty(get_reports(analyzer))
    end

    # more cache test, constant prop should re-run in deeper level
    let
        m = @fixturedef begin
            mutable struct P
                i::Int
                s::String
            end
            foo(p, i) = p.i = i
            bar(args...) = foo(args...)
        end

        # "for one of the union split cases, no matching method found for signature: Base.convert(Base.fieldtype(Base.typeof(x::P)::Type{P}, f::Symbol)::Union{Type{Int64}, Type{String}}, v::Int64)" should be threw away
        analyzer, = Core.eval(m, :($report_call(bar, (P, Int))))
        @test isempty(get_reports(analyzer))

        # works for cache
        analyzer, = Core.eval(m, :($report_call(bar, (P, Int))))
        @test isempty(get_reports(analyzer))
    end

    # constant prop should not exclude those are not related
    let
        m = gen_virtual_module()
        analyzer, = Core.eval(m, quote
            mutable struct P
                i::Int
                s::String
            end
            function foo(p, i, s)
                p.i = i
                p.s = s
            end

            $report_call(foo, (P, Int, #= invalid =# Int))
        end)

        # "for one of the union split cases, no matching method found for signature: Base.convert(Base.fieldtype(Base.typeof(x::P)::Type{P}, f::Symbol)::Union{Type{Int64}, Type{String}}, v::Int64)" should be threw away, while
        # "no matching method found for call signature: Base.convert(Base.fieldtype(Base.typeof(x::P)::Type{P}, f::Symbol)::Type{String}, v::Int64)" should be kept
        @test length(get_reports(analyzer)) === 1
        er = first(get_reports(analyzer))
        @test er isa NoMethodErrorReport
        @test er.t === Tuple{typeof(convert), Type{String}, Int}
    end

    # constant prop should narrow down union-split no method error to single no method matching error
    let
        m = gen_virtual_module()
        analyzer, = Core.eval(m, quote
            mutable struct P
                i::Int
                s::String
            end
            function foo(p, i, s)
                p.i = i
                p.s = s
            end

            $report_call(foo, (P, String, Int))
        end)

        # "for one of the union split cases, no matching method found for signature: Base.convert(Base.fieldtype(Base.typeof(x::P)::Type{P}, f::Symbol)::Union{Type{Int64}, Type{String}}, v::String)" should be narrowed down to "no matching method found for call signature: Base.convert(Base.fieldtype(Base.typeof(x::P)::Type{P}, f::Symbol)::Type{Int}, v::String)"
        @test !isempty(get_reports(analyzer))
        @test any(get_reports(analyzer)) do report
            report isa NoMethodErrorReport &&
            report.t === Tuple{typeof(convert), Type{Int}, String}
        end
        # "no matching method found for call signature: Base.convert(Base.fieldtype(Base.typeof(x::P)::Type{P}, f::Symbol)::Type{String}, v::Int)"
        # won't be reported since `typeinf` early escapes on `Bottom`-annotated statement
    end

    # report-throw away with constant analysis shouldn't throw away reports from the same
    # frame but with the other constants
    let
        m = gen_virtual_module()
        analyzer, = Core.eval(m, quote
            foo(a) = a<0 ? a+string(a) : a
            bar() = foo(-1), foo(1) # constant analysis on `foo(1)` shouldn't throw away reports from `foo(-1)`
            $report_call(bar)
        end)
        @test !isempty(get_reports(analyzer))
        @test any(r->isa(r,NoMethodErrorReport), get_reports(analyzer))
    end

    let
        m = gen_virtual_module()
        analyzer, = Core.eval(m, quote
            foo(a) = a<0 ? a+string(a) : a
            function bar(b)
                a = b ? foo(-1) : foo(1)
                b = foo(-1)
                return a, b
            end
            $report_call(bar, (Bool,))
        end)
        @test !isempty(get_reports(analyzer))
        # FIXME our report uniquify logic might be wrong and it wrongly singlifies the different reports here
        @test_broken count(isa(report, NoMethodErrorReport) for report in get_reports(analyzer)) == 2
    end

    @testset "constant analysis throws away false positive reports" begin
        let
            m = @fixturedef begin
                foo(a) = a > 0 ? a : "minus"
                bar(a) = foo(a) + 1
            end

            # constant propagation can reveal the error pass can't happen
            analyzer, = Core.eval(m, :($report_call(()->bar(10))))
            @test isempty(get_reports(analyzer))

            # for this case, no constant prop' doesn't happen, we can't throw away error pass
            analyzer, = Core.eval(m, :($report_call(bar, (Int,))))
            @test length(get_reports(analyzer)) === 1
            er = first(get_reports(analyzer))
            @test er isa NoMethodErrorReport
            @test er.t == [Tuple{typeof(+),String,Int}]

            # if we run constant prop' that leads to the error pass, we should get the reports
            analyzer, = Core.eval(m, :($report_call(()->bar(0))))
            @test length(get_reports(analyzer)) === 1
            er = first(get_reports(analyzer))
            @test er isa NoMethodErrorReport
            @test er.t === Tuple{typeof(+),String,Int}
        end

        # we should throw-away reports collected from frames that are revealed as "unreachable"
        # by constant prop'
        let
            m = @fixturedef begin
                foo(a) = bar(a)
                function bar(a)
                    return if a < 1
                        baz1(a, "0")
                    else
                        baz2(a, a)
                    end
                end
                baz1(a, b) = a ? b : b
                baz2(a, b) = a + b
            end

            # no constant prop, just report everything
            analyzer, = Core.eval(m, :($report_call(foo, (Int,))))
            @test length(get_reports(analyzer)) === 1
            er = first(get_reports(analyzer))
            @test er isa NonBooleanCondErrorReport &&
                er.t === Int

            # constant prop should throw away the non-boolean condition report from `baz1`
            analyzer, = Core.eval(m, quote
                $report_call() do
                    foo(1)
                end
            end)
            @test isempty(get_reports(analyzer))

            # constant prop'ed, still we want to have the non-boolean condition report from `baz1`
            analyzer, = Core.eval(m, quote
                $report_call() do
                    foo(0)
                end
            end)
            @test length(get_reports(analyzer)) === 1
            er = first(get_reports(analyzer))
            @test er isa NonBooleanCondErrorReport &&
                er.t === Int

            # so `Bool` is good for `foo` after all
            analyzer, = Core.eval(m, :($report_call(foo, (Bool,))))
            @test isempty(get_reports(analyzer))
        end

        # end to end
        let
            res = @analyze_toplevel begin
                function foo(n)
                    if n < 10
                        return n
                    else
                        return "over 10"
                    end
                end

                function bar(n)
                    if n < 10
                        return foo(n) + 1
                    else
                        return foo(n) * "+1"
                    end
                end

                bar(1)
                bar(10)
            end

            @test isempty(res.inference_error_reports)
        end
    end
end

@testset "keyword argument methods" begin
    analyzer, = Core.eval(gen_virtual_module(), quote
        f(a; b = nothing, c = nothing) = return
        $report_call((Any,)) do b
            f(1; b)
        end
    end)
    @test isempty(get_reports(analyzer))
end

@testset "don't early escape if type grows up to `Any`" begin
    vmod = gen_virtual_module()
    res = @analyze_toplevel context = vmod virtualize = false begin
        abstract type Foo end
        struct Foo1 <: Foo
            bar
        end
        struct Foo2 <: Foo
            baz # typo
        end
        struct Foo3 <: Foo
            qux # typo
        end
        bar(foo::Foo) = foo.bar

        f = rand(Bool) ? Foo1(1) : rand(Bool) ? Foo2(1) : Foo3(1)
        bar(f)
    end

    # `bar(::Foo1)` is valid but its return type is `Any`, and so we can't collect possible
    # error points for `bar(::Foo2)` and `bar(::Foo3)` if we bail out on `Any`-grew return type
    @test length(res.inference_error_reports) === 2
    @test any(res.inference_error_reports) do er
        return er isa NoFieldErrorReport &&
            er.typ === vmod.Foo2 &&
            er.name === :bar
    end
    @test any(res.inference_error_reports) do er
        return er isa NoFieldErrorReport &&
            er.typ === vmod.Foo3 &&
            er.name === :bar
    end
end

@static isdefined(CC, :abstract_invoke) && @testset "abstract_invoke" begin
    # non-`Type` `argtypes`
    analyzer, = report_call() do
        invoke(sin, :this_should_be_type, 1.0)
    end
    @test length(get_reports(analyzer)) == 1
    r = first(get_reports(analyzer))
    @test isa(r, InvalidInvokeErrorReport)

    # invalid `argtype`
    analyzer, = report_call() do
        Base.@invoke sin(1.0::Int)
    end
    @test length(get_reports(analyzer)) == 1
    r = first(get_reports(analyzer))
    @test isa(r, InvalidInvokeErrorReport)

    # don't report errors collected in `invoke`d functions
    foo(i::Integer) = throw(string(i))
    analyzer, = report_call((Any,)) do a
        Base.@invoke foo(a::Integer)
    end
    @test !isempty(get_reports(analyzer))
    @test !any(r->isa(r, InvalidInvokeErrorReport), get_reports(analyzer))

    bar(a) = return a + undefvar
    function baz(a)
        a += 1
        Base.@invoke bar(a::Int) # `invoke` is valid, but error should happne within `bar`
    end
    analyzer, = report_call(baz, (Any,))
    @test !isempty(get_reports(analyzer))
    @test !any(r->isa(r, InvalidInvokeErrorReport), get_reports(analyzer))
    # virtual stack trace should include frames of both `bar` and `baz`
    # we already `@assert` it within `typeinf` when `JET_DEV_MODE` is enabled, but test it
    # here just to make sure it's working
    # WARNING the code block below is really line sensitive, in a relative way to the definitions of `bar` and `baz`
    @test all(get_reports(analyzer)) do r
        any(r.vst) do vf
            vf.file === Symbol(@__FILE__) &&
            vf.line == (@__LINE__) - 15 && # `bar`
            vf.linfo.def.name === :bar
        end || return false
        any(r.vst) do vf
            vf.file === Symbol(@__FILE__) &&
            vf.line == (@__LINE__) - 17 && # `baz`
            vf.linfo.def.name === :baz
        end || return false
        return true
    end
end

@testset "additional analysis pass for task parallelism code" begin
    # general case with `schedule(::Task)` pattern
    analyzer, = report_call() do
        t = Task() do
            sum("julia")
        end
        schedule(t)
        fetch(t)
    end
    test_sum_over_string(analyzer)

    # handle `Threads.@spawn` (https://github.com/aviatesk/JET.jl/issues/114)
    analyzer, = report_call() do
        fetch(Threads.@spawn 1 + "foo")
    end
    @test length(get_reports(analyzer)) == 1
    let r = first(get_reports(analyzer))
        @test isa(r, NoMethodErrorReport)
        @test r.t === Tuple{typeof(+), Int, String}
    end

    # handle `Threads.@threads`
    analyzer, = report_call((Int,)) do n
        a = String[]
        Threads.@threads for i in 1:n
            push!(a, i)
        end
        return a
    end
    @test !isempty(get_reports(analyzer))
    @test any(get_reports(analyzer)) do r
        isa(r, NoMethodErrorReport) &&
        r.t === Tuple{typeof(convert), Type{String}, Int}
    end

    # multiple tasks in the same frame
    analyzer, = report_call() do
        t1 = Threads.@spawn 1 + "foo"
        t2 = Threads.@spawn "foo" + 1
        fetch(t1), fetch(t2)
    end
    @test length(get_reports(analyzer)) == 2
    let r = get_reports(analyzer)[1]
        @test isa(r, NoMethodErrorReport)
        @test r.t === Tuple{typeof(+), Int, String}
    end
    let r = get_reports(analyzer)[2]
        @test isa(r, NoMethodErrorReport)
        @test r.t === Tuple{typeof(+), String, Int}
    end

    # nested tasks
    analyzer, = report_call() do
        t0 = Task() do
            t = Threads.@spawn sum("julia")
            fetch(t)
        end
        schedule(t0)
        fetch(t0)
    end
    test_sum_over_string(analyzer)

    # when `schedule` call is separated from `Task` definition
    make_task(s) = Task() do
        sum(s)
    end
    function run_task(t)
        schedule(t)
        fetch(t)
    end
    analyzer, = report_call() do
        t = make_task("julia")

        run_task(t)
    end
    test_sum_over_string(analyzer)
    let r = first(get_reports(analyzer))
        # we want report to come from `run_task`, but currently we invoke JET analysis on `Task` construction
        @test_broken any(r.vst) do vf
            vf.linfo.def.name === :run_task
        end
    end

    # report uncaught exception happened in a task
    analyzer, = report_call() do
        fetch(Threads.@spawn throw("foo"))
    end
    @test length(get_reports(analyzer)) == 1
    @test isa(first(get_reports(analyzer)), UncaughtExceptionReport)

    # don't fail into infinite loop (rather, don't spoil inference termination)
    m = @fixturedef begin
        # adapated from https://julialang.org/blog/2019/07/multithreading/
        import Base.Threads.@spawn

        # sort the elements of `v` in place, from indices `lo` to `hi` inclusive
        function psort!(v, lo::Int=1, hi::Int=length(v))
            if lo >= hi                       # 1 or 0 elements; nothing to do
                return v
            end
            if hi - lo < 100000               # below some cutoff, run in serial
                sort!(view(v, lo:hi), alg = MergeSort)
                return v
            end

            mid = (lo+hi)>>>1                 # find the midpoint

            half = @spawn psort!(v, lo, mid)  # task to sort the lower half; will run
            psort!(v, mid+1, hi)              # in parallel with the current call sorting
                                              # the upper half
            wait(half)                        # wait for the lower half to finish

            temp = v[lo:mid]                  # workspace for merging

            i, k, j = 1, lo, mid+1            # merge the two sorted sub-arrays
            @inbounds while k < j <= hi
                if v[j] < temp[i]
                    v[k] = v[j]
                    j += 1
                else
                    v[k] = temp[i]
                    i += 1
                end
                k += 1
            end
            @inbounds while k < j
                v[k] = temp[i]
                k += 1
                i += 1
            end

            return v
        end
    end
    analyzer, = report_call(m.psort!, (Vector{Int},))
    @test true
end

@static if VERSION ≥ v"1.7.0-DEV.705"

@testset "opaque closure" begin
    # can cache const prop' result with varargs
    function oc_varargs_constprop()
        oc = Base.Experimental.@opaque (args...)->args[1]+args[2]+arg[3] # typo on `arg[3]`
        return Val{oc(1,2,3)}()
    end
    analyzer, = @report_call oc_varargs_constprop()
    @test !isempty(get_cache(analyzer))
end

end # @static if VERSION ≥ v"1.7.0-DEV.705"

@testset "https://github.com/aviatesk/JET.jl/issues/133" begin
    res = @analyze_toplevel begin
        @ccall strlen("foo"::Cstring)::Csize_t
    end
    @test isempty(res.inference_error_reports)
end

@testset "staged programming" begin
    m = @fixturedef begin
        @generated function foo(a)
            if a <: Integer
                return :(a)
            elseif a <: Number
                return :(undefvar) # report me this case
            end
            throw("invalid argument")
        end
    end

    # successful code generation, valid code
    let
        analyzer, = report_call(m.foo, (Int,))
        @test isempty(get_reports(analyzer))
    end

    # successful code generation, invalid code
    let
        analyzer, = report_call(m.foo, (Float64,))
        @test length(get_reports(analyzer)) == 1
        r = first(get_reports(analyzer))
        @test isa(r, GlobalUndefVarErrorReport)
        @test r.name === :undefvar
    end

    # unsuccessful code generation
    let
        analyzer, = report_call(m.foo, (String,))
        @test length(get_reports(analyzer)) == 1
        r = first(get_reports(analyzer))
        @test isa(r, GeneratorErrorReport) && r.err == "invalid argument"
    end
end


text:
config = MOI.Test.TestConfig()
optimize!(mock) = MOIU.mock_optimize!(mock, [[128.0, 0.0, -256.0, 0.0, 160.0, 0.0, -32.0, 0.0, 1.0]; zeros(90)])
for mock in mocks(optimize!)
    Tests.chebyshev_test(mock, config)
end


text:
function _randomping()

    ping = true

    function loop()
        count = 1
        t0 = time()
        while ping
            sleep(1)
            msg = "Ping $(count): $(round(time()-t0,digits=2)) seconds"
            e = Pages.Element(id = string(count),innerHTML = msg)
            e.style["color"] = rand(["red","green","blue"])
            rand() < .5 && Pages.append(e)

            count += 1
            tstart = time()
        end
    end

    function start()
        ping = true

        @async loop()
        return
    end

    function stop()
        ping = false
        return
    end

    return (start=start,stop=stop)
end

randomping = _randomping()


text:
using Router, Channels

channel("/todos/subscribe", "todos#TodosChannel.subscribe", named = :todos)

max_stars_repo_name: essenciary/genie-todo-mvc

max_issues_repo_name: essenciary/genie-todo-mvc

max_forks_repo_name: essenciary/genie-todo-mvc


text:
#
# This file is part of the Guards.jl Julia package,
# MIT license, part of https://github.com/JuliaActors
#

using Guards, Test

gd = guard([1,2,3])
@test gd[2] == 2
@test gd[begin] == 1
@test gd[end] == 3
@test (gd[1] = 10) == 10
@test (@grd gd) == [10,2,3]
@test supertype(typeof(gd)) == Array{Int64, 1}
@test size(gd) == (3,)
@test IndexStyle(typeof(gd)) == IndexLinear()
@test length(gd) == 3

dd = guard(Dict("baz" => 17, "bar" => 4711))
@test dd["baz"] == 17
@test (dd["dummy"] = 10) == 10
@test dd["dummy"] == 10


text:
function propagator(f)
  function (cells...)
    output = last(cells)
    inputs = cells[1:end-1]
    function activate!()
      inputvalues = map(getindex, inputs)
      if !any(isunusable, inputvalues)
        addcontent!(output, f(inputvalues...))
      end
    end
    return Propagator(activate!, inputs)
  end
end

function block(tobuild, inputs...)
  isbuilt = false
  function maybebuild()
    if !(isbuilt || all(isunusable, map(getindex, inputs)))
      tobuild()
      isbuilt = true
    end
  end
  return Propagator(maybebuild, inputs)
end

constant(x) = propagator(() -> x)

function conditional(p, t, f, out)
  Propagator([p, t, f]) do
    predicate = p[]
    predicate != nothing || return
    out[] = predicate ? t[] : f[]
  end
end

switch(p, t, out) = conditional(p, t, Cell(), out)

function constraint(::typeof(*))
  function (x, y, xy)
    propagator(*)(x, y, xy)
    propagator(/)(xy, x, y)
    propagator(/)(xy, y, x)
  end
end

function constraint(::typeof(+))
  function (x, y, xy)
    propagator(+)(x, y, xy)
    propagator(-)(xy, x, y)
    propagator(-)(xy, y, x)
  end
end

function quadratic(x, x²)
  propagator(x -> x^2)(x, x²)
  propagator(sqrt)(x², x)
end


text:
module ConvexHull

using JuMP, GenericLinearAlgebra, RowEchelon, Combinatorics

using Polyhedra
using Printf: @printf

export homogeneous_system, read_ine, write_ine, read_ext, write_ext,
       double_description, canonicalize!,
       get_extrema, is_approx_included

const ε = 10eps()

include("double_description.jl")
include("jump_interface.jl")
include("readers_writers.jl")
include("polyhedron.jl")

@static if VERSION < v"1.1"
    @eval isnothing(x) = x === nothing
end

end


text:
module MPIStateArrays
using LinearAlgebra
using DoubleFloats

using MPI

export MPIStateArray, euclidean_distance, weightedsum

"""
    MPIStateArray{S <: Tuple, T, DeviceArray, N,
                  DATN<:AbstractArray{T,N}, Nm1, DAI1} <: AbstractArray{T, N}

`N`-dimensional MPI-aware array with elements of type `T`. The dimension `N` is
`length(S) + 1`. `S` is a tuple of the first `N-1` array dimensions.

!!! todo

    It should be reevaluated whether all this stuff in the type domain is
    really necessary (some of it was optimistically added for functionality that
    never panned out)

"""
struct MPIStateArray{S <: Tuple, T, DeviceArray, N,
                     DATN<:AbstractArray{T,N}, Nm1, DAI1} <: AbstractArray{T, N}
  mpicomm::MPI.Comm
  Q::DATN

  realelems::UnitRange{Int64}
  ghostelems::UnitRange{Int64}
  sendelems::DAI1

  sendreq::Array{MPI.Request, 1}
  recvreq::Array{MPI.Request, 1}

  host_sendQ::Array{T, N}
  host_recvQ::Array{T, N}

  nabrtorank::Array{Int64, 1}
  nabrtorecv::Array{UnitRange{Int64}, 1}
  nabrtosend::Array{UnitRange{Int64}, 1}

  device_sendQ::DATN
  device_recvQ::DATN

  # FIXME: Later we should relax this if we compute on the GPU and probably
  # should let this be a more generic type...
  weights::Array{T, Nm1}

  commtag::Int
  function MPIStateArray{S, T, DA}(mpicomm, numelem, realelems, ghostelems,
                                   sendelems, nabrtorank, nabrtorecv,
                                   nabrtosend, weights, commtag
                                  ) where {S, T, DA}
    N = length(S.parameters)+1
    numsendelem = length(sendelems)
    numrecvelem = length(ghostelems)
    (Q, device_sendQ, device_recvQ) = try
      (DA{T, N}(undef, S.parameters..., numelem),
       DA{T, N}(undef, S.parameters..., numsendelem),
       DA{T, N}(undef, S.parameters..., numrecvelem))
    catch
      try
        # TODO: Remove me after CUDA upgrade...
        (DA{T, N}(S.parameters..., numelem),
         DA{T, N}(S.parameters..., numsendelem),
         DA{T, N}(S.parameters..., numrecvelem))
      catch
        error("MPIStateArray:Cannot construct array")
      end
    end

    host_sendQ = zeros(T, S.parameters..., numsendelem)
    host_recvQ = zeros(T, S.parameters..., numrecvelem)

    # Length check is to work around a CuArrays bug.
    length(Q) > 0 && fill!(Q, 0)
    length(device_sendQ) > 0 && fill!(device_sendQ, 0)
    length(device_recvQ) > 0 && fill!(device_recvQ, 0)

    nnabr = length(nabrtorank)
    sendreq = fill(MPI.REQUEST_NULL, nnabr)
    recvreq = fill(MPI.REQUEST_NULL, nnabr)

    sendelems = typeof(sendelems) <: DA ? sendelems : DA(sendelems)
    DAI1 = typeof(sendelems)
    new{S, T, DA, N, typeof(Q), N-1, DAI1}(mpicomm, Q, realelems, ghostelems,
                                           sendelems, sendreq, recvreq,
                                           host_sendQ, host_recvQ, nabrtorank,
                                           nabrtorecv, nabrtosend,
                                           device_sendQ, device_recvQ, weights,
                                           commtag)
  end

end

"""
   MPIStateArray{S, T, DA}(mpicomm, numelem; realelems=1:numelem,
                           ghostelems=numelem:numelem-1,
                           sendelems=1:0,
                           nabrtorank=Array{Int64}(undef, 0),
                           nabrtorecv=Array{UnitRange{Int64}}(undef, 0),
                           nabrtosend=Array{UnitRange{Int64}}(undef, 0),
                           weights,
                           commtag=888)

Construct an `MPIStateArray` over the communicator `mpicomm` with `numelem`
elements, using array type `DA` with element type `eltype`. The arrays that are
held in this created `MPIStateArray` will be of size `(S..., numelem)`.

The range `realelems` is the number of elements that this mpirank owns, whereas
the range `ghostelems` is the elements that are owned by other mpiranks.
Elements are stored as 'realelems` followed by `ghostelems`.

  * `sendelems` is an ordered array of elements to be sent to neighboring
    mpiranks
  * `nabrtorank` is the list of neighboring mpiranks
  * `nabrtorecv` is an `Array` of `UnitRange` that give the `ghostelems`
    received from neighboring mpiranks (indexes into the ghost elements arrays,
    not the full element array)
  * nabrtosend` is an `Array` of `UnitRange` for which elements to send to
    which neighboring mpiranks indexing into the `sendelems` ordering
  * `weights` is an optional array which gives weight for each degree of freedom
    to be used when computing the 2-norm of the array
"""
function MPIStateArray{S, T, DA}(mpicomm, numelem;
                                 realelems=1:numelem,
                                 ghostelems=numelem:numelem-1,
                                 sendelems=1:0,
                                 nabrtorank=Array{Int64}(undef, 0),
                                 nabrtorecv=Array{UnitRange{Int64}}(undef, 0),
                                 nabrtosend=Array{UnitRange{Int64}}(undef, 0),
                                 weights=nothing,
                                 commtag=888
                                ) where {S<:Tuple, T, DA}

  N = length(S.parameters)+1
  if weights == nothing
    weights = Array{T}(undef, ntuple(j->0, N-1))
  elseif !(typeof(weights) <: Array)
    weights = Array(weights)
  end
  MPIStateArray{S, T, DA}(mpicomm, numelem, realelems, ghostelems,
                          sendelems, nabrtorank, nabrtorecv,
                          nabrtosend, weights, commtag)
end

# FIXME: should general cases should be handled?
function Base.similar(Q::MPIStateArray{S, T, DA}; commtag=Q.commtag
                     ) where {S, T, DA}
  MPIStateArray{S, T, DA}(Q.mpicomm, size(Q.Q)[end], Q.realelems, Q.ghostelems,
                          Q.sendelems, Q.nabrtorank, Q.nabrtorecv,
                          Q.nabrtosend, Q.weights, commtag)
end

# FIXME: Only show real size
Base.size(Q::MPIStateArray, x...;kw...) = size(Q.Q, x...;kw...)

# FIXME: Only let get index access real elements?
Base.getindex(Q::MPIStateArray, x...;kw...) = getindex(Q.Q, x...;kw...)

# FIXME: Only let set index access real elements?
Base.setindex!(Q::MPIStateArray, x...;kw...) = setindex!(Q.Q, x...;kw...)

Base.eltype(Q::MPIStateArray, x...;kw...) = eltype(Q.Q, x...;kw...)

Base.Array(Q::MPIStateArray) = Array(Q.Q)

function Base.copyto!(dst::MPIStateArray, src::Array)
  copyto!(dst.Q, src)
  dst
end
Base.copyto!(dst::Array, src::MPIStateArray) = copyto!(dst, src.Q)
function Base.copyto!(dst::MPIStateArray, src::MPIStateArray)
  copyto!(dst.Q, src.Q)
  dst
end

"""
    post_Irecvs!(Q::MPIStateArray)

posts the `MPI.Irecv!` for `Q`
"""
function post_Irecvs!(Q::MPIStateArray)
  nnabr = length(Q.nabrtorank)
  for n = 1:nnabr
    # If this fails we haven't waited on previous recv!
    @assert Q.recvreq[n].buffer == nothing

    Q.recvreq[n] = MPI.Irecv!((@view Q.host_recvQ[:, :, Q.nabrtorecv[n]]),
                              Q.nabrtorank[n], Q.commtag, Q.mpicomm)
  end
end

"""
    start_ghost_exchange!(Q::MPIStateArray; dorecvs=true)

Start the MPI exchange of the data stored in `Q`. If `dorecvs` is `true` then
`post_Irecvs!(Q)` is called, otherwise the caller is responsible for this.

This function will fill the send buffer (on the device), copies the data from
the device to the host, and then issues the send. Previous sends are waited on
to ensure that they are complete.
"""
function start_ghost_exchange!(Q::MPIStateArray; dorecvs=true)
  dorecvs && post_Irecvs!(Q)

  # wait on (prior) MPI sends
  finish_ghost_send!(Q)

  # pack data in send buffer
  fillsendbuf!(Q.host_sendQ, Q.device_sendQ, Q.Q, Q.sendelems)

  # post MPI sends
  nnabr = length(Q.nabrtorank)
  for n = 1:nnabr
    Q.sendreq[n] = MPI.Isend((@view Q.host_sendQ[:, :, Q.nabrtosend[n]]),
                           Q.nabrtorank[n], Q.commtag, Q.mpicomm)
  end
end

"""
    finish_ghost_exchange!(Q::MPIStateArray)

Complete the exchange of data and fill the data array on the device. Note this
completes both the send and the receive communication. For more fine level
control see [finish_ghost_exchange_recv!](@ref) and
[finish_ghost_exchange_send!](@ref)
"""
function finish_ghost_exchange!(Q::MPIStateArray)
  finish_ghost_recv!(Q::MPIStateArray)
  finish_ghost_send!(Q::MPIStateArray)
end

"""
    finish_ghost_recv!(Q::MPIStateArray)

Complete the receive of data and fill the data array on the device
"""
function finish_ghost_recv!(Q::MPIStateArray)
  # wait on MPI receives
  MPI.Waitall!(Q.recvreq)

  # copy data to state vectors
  transferrecvbuf!(Q.device_recvQ, Q.host_recvQ, Q, length(Q.realelems))
end

"""
    finish_ghost_send!(Q::MPIStateArray)

Waits on the send of data to be complete
"""
finish_ghost_send!(Q::MPIStateArray) = MPI.Waitall!(Q.sendreq)

# {{{ MPI Buffer handling
fillsendbuf!(h, d, b::MPIStateArray, e) = fillsendbuf!(h, d, b.Q, e)
transferrecvbuf!(h, d, b::MPIStateArray, e) = transferrecvbuf!(h, d, b.Q, e)

function fillsendbuf!(host_sendbuf, device_sendbuf::Array, buf::Array, sendelems)
  host_sendbuf[:, :, :] .= buf[:, :, sendelems]
end

function transferrecvbuf!(device_recvbuf::Array, host_recvbuf, buf::Array,
                          nrealelem)
  buf[:, :, nrealelem+1:end] .= host_recvbuf[:, :, :]
end
# }}}

# Integral based metrics
function LinearAlgebra.norm(Q::MPIStateArray; p::Real=2)

  @assert p == 2

  host_array = Array ∈ typeof(Q).parameters
  h_Q = host_array ? Q : Array(Q)
  Np = size(Q, 1)

  if isempty(Q.weights)
    locnorm2 = knl_norm2(Val(Np), h_Q, Q.realelems)
  else
    locnorm2 = knl_L2norm(Val(Np), h_Q, Q.weights, Q.realelems)
  end

  sqrt(MPI.Allreduce([locnorm2], MPI.SUM, Q.mpicomm)[1])
end

function knl_norm2(::Val{Np}, Q, elems) where {Np}
  DFloat = eltype(Q)
  (_, nstate, nelem) = size(Q)

  nrm = zero(DFloat)

  @inbounds for e = elems, q = 1:nstate, i = 1:Np
    nrm += Q[i, q, e]^2
  end

  nrm
end

function knl_L2norm(::Val{Np}, Q, weights, elems) where {Np}
  DFloat = eltype(Q)
  (_, nstate, nelem) = size(Q)

  nrm = zero(DFloat)

  @inbounds for e = elems, q = 1:nstate, i = 1:Np
    nrm += weights[i, e] * Q[i, q, e]^2
  end

  nrm
end

function euclidean_distance(A::MPIStateArray, B::MPIStateArray)

  host_array = Array ∈ typeof(A).parameters
  h_A = host_array ? A : Array(A)
  Np = size(A, 1)

  host_array = Array ∈ typeof(B).parameters
  h_B = host_array ? B : Array(B)
  @assert Np === size(B, 1)

  if isempty(A.weights)
    locdist = knl_dist(Val(Np), h_A, h_B, A.realelems)
  else
    locdist = knl_L2dist(Val(Np), h_A, h_B, A.weights, A.realelems)
  end

  sqrt(MPI.Allreduce([locdist], MPI.SUM, A.mpicomm)[1])
end

function knl_dist(::Val{Np}, A, B, elems) where {Np}
  DFloat = eltype(A)
  (_, nstate, nelem) = size(A)

  dist = zero(DFloat)

  @inbounds for e = elems, q = 1:nstate, i = 1:Np
    dist += (A[i, q, e] - B[i, q, e])^2
  end

  dist
end

function knl_L2dist(::Val{Np}, A, B, weights, elems) where {Np}
  DFloat = eltype(A)
  (_, nstate, nelem) = size(A)

  dist = zero(DFloat)

  @inbounds for e = elems, q = 1:nstate, i = 1:Np
    dist += weights[i, e] * (A[i, q, e] - B[i, q, e])^2
  end

  dist
end

"""
    weightedsum(A[, states])

Compute the weighted sum of the `MPIStateArray` `A`. If `states` is specified on
the listed states are summed, otherwise all the states in `A` are used.

A typical use case for this is when the weights have been initialized with
quadrature weights from a grid, thus this becomes an integral approximation.

!!! note

    This implementation is not optimal and should be revisited when we work out
    the GPUify version!

"""
function weightedsum(A::MPIStateArray, states=1:size(A, 2))

  host_array = Array ∈ typeof(A).parameters
  h_A = host_array ? A : Array(A)
  Np = size(A, 1)

  isempty(A.weights) && error("`weightedsum` requires weights")
  locwsum = knl_weightedsum(Val(Np), h_A, A.weights, A.realelems, states)

  DFloat = eltype(A)
  # Need to use anomous function version of sum else MPI.jl using MPI_SUM
  DFloat(MPI.Allreduce([locwsum], (x,y)->x+y, A.mpicomm)[1])
end

function knl_weightedsum(::Val{Np}, A, weights, elems, states) where {Np}
  DFloat = eltype(A)

  wsum = zero(DoubleFloat{DFloat})

  @inbounds for e = elems
    for q = states, i = 1:Np
      wsum += weights[i, e] * A[i, q, e]
    end
  end

  wsum
end

using Requires

@init @require CuArrays = "3a865a2d-5b23-5a0f-bc46-62713ec82fae" begin
  using .CuArrays
  using .CuArrays.CUDAnative
  using .CuArrays.CUDAnative.CUDAdrv

  include("MPIStateArrays_cuda.jl")
end

end

max_stars_repo_name: christophernhill-dev/CLIMA

max_issues_repo_name: christophernhill-dev/CLIMA

max_forks_repo_name: christophernhill-dev/CLIMA


text:
using Timestamps, MarketData, FactCheck

include("timestamp.jl")
include("operators.jl")
include("conversion.jl")
include("arraymethods.jl")
include("utils.jl")

exitstatus()


text:
# We use getters so the model never accesses struct fields directly,
# except in methods for specific leaf types where declaring all the gettier
# methods would be overkill.
#
# A benefit of this is parameters can be easily moved between shared and specific
# parameter structs.

# So we don't have to depend on all of Lazy.jl
macro forward(ex, fs)
  @capture(ex, T_.field_) || error("Syntax: @forward T.x f, g, h")
  T = esc(T)
  fs = isexpr(fs, :tuple) ? map(esc, fs.args) : [esc(fs)]
  :($([:($f(x::$T, args...) = (Base.@_inline_meta; $f(x.$field, args...)))
       for f in fs]...);
    nothing)
end

@inline dead(o::AbstractOrganism) = o.dead[]
@inline set_dead!(o::AbstractOrganism, val) = o.dead[] = val
@inline environment(o::AbstractOrganism) = o.environment

@inline vars(o::AbstractOrgan) = o.vars
@inline flux(o::AbstractOrgan) = o.J
@inline flux1(o::AbstractOrgan) = o.J1

@inline tstep(v) = v.t[1]
@inline set_tstep!(v, val) = v.t[1] = val

const varfuncs = [:shape, :rate, :temp, :θE, :tempcorrection, :height]
for (v, sv) in zip(varfuncs, Symbol.(:set_, varfuncs, :!))
    @eval @inline ($v)(vars) = vars.$v[tstep(vars)]
    @eval @inline ($sv)(vars, val) = vars.$v[tstep(vars)] = val
end

@inline depth(v) = height(v)

rate_formula(p) = p.rate_formula
assimilation_pars(p) = p.assimilation_pars
shape_pars(p) = p.shape_pars
allometry_pars(p) = p.allometry_pars
maturity_pars(p) = p.maturity_pars
trans_pars(p) = p.trans_pars
production_pars(p) = p.production_pars
rejection_pars(p) = p.rejection_pars
germination_pars(p) = p.germination_pars
turnover_pars(p) = p.turnover_pars
feedback_pars(p) = p.feedback_pars
tempcorr_pars(p) = p.tempcorr_pars
core_pars(p) = p.core_pars
su_pars(p) = p.su_pars
catabolism_pars(p) = p.catabolism_pars
maintenance_pars(p) = p.maintenance_pars

@forward AbstractOrgan.vars θE, temp, set_temp!, tempcorrection, set_tempcorrection!,
         height, set_height!, rate, set_rate!, shape, set_shape!, tstep, set_tstep!

@forward AbstractOrgan.params rate_formula, assimilation_pars, shape_pars, allometry_pars, maturity_pars,
                              trans_pars, production_pars, rejection_pars, germination_pars, turnover_pars

@forward AbstractOrgan.shared maintenance_pars, feedback_pars, su_pars, tempcorr_pars, catabolism_pars, core_pars,
                              y_V_E, y_E_EC, y_E_EN, n_N_P, n_N_V, n_N_E, n_N_EC, n_N_EN, w_V, w_C, w_N, w_E

n_N_P(p) = production_pars(p).n_N_P
y_V_E(p) = core_pars(p).y_V_E
y_E_EC(p) = core_pars(p).y_E_EC
y_E_EN(p) = core_pars(p).y_E_EN
n_N_V(p) = core_pars(p).n_N_V
n_N_E(p) = core_pars(p).n_N_E
n_N_EC(p) = core_pars(p).n_N_EC
n_N_EN(p) = core_pars(p).n_N_EN
w_V(p) = core_pars(p).w_V
w_C(p) = core_pars(p).w_C
w_N(p) = core_pars(p).w_N
w_E(p) = core_pars(p).w_E

j_E_mai(o::AbstractOrgan) = j_E_mai(maintenance_pars(o))

κtra(o::AbstractOrgan) = κtra(trans_pars(o))
κtra(o::Nothing) = 0.0

κmat(o::AbstractOrgan) = κmat(maturity_pars(o))
κmat(::Nothing) = 0.0

κsoma(o::AbstractOrgan) = oneunit(κtra(o)) - κtra(o) - κmat(o)

mass(o, u) = u.V * w_V(o)


text:
"""
    SimplexLatticeDesign( Components::Int, Spaces::Int )

Returns an array of tuples (`Components` in length) which represent design points.
"""
function SimplexLatticeDesign( Components::Int, Spaces::Int )
    @assert( ( Components > 0 ) && ( Spaces > 0 ), "Components and Spaces must be positive integers!" )
    DesignPoints = Inf
    try
        DesignPoints = Int( factorial( Components + Spaces - 1) / ( factorial(Spaces) * factorial(Components - 1) ))
    catch
        @warn("Cannot estimate number of design points. May execute for the lifetime of the universe.")
    end
    PossiblePoints = reverse(0:Spaces) ./ Spaces
    DesignSpace = repeat([PossiblePoints], Components)
    Possible = [ prod for prod in Base.product(DesignSpace...) if (sum(prod) == 1) ]

    if !isinf(DesignPoints)
        @assert( DesignPoints == length(Possible) )
    end
    return Possible
end

"""
    SimplexCentroidDesign( Components::Int, Order::Union{ UnitRange{ Int }, Int} )

Returns an array of tuples (`Components` in length) which represent design points.

Please note: the Order can be an Int or a UnitRange.
"""
function SimplexCentroidDesign( Components::Int, Order::Union{ UnitRange{ Int }, Int} )
    HighestOrder = Order[end]
    @assert( ( Components > 0 ) && ( HighestOrder > 0 ), "Components and the highest orders must be positive integers!" )
    @assert( HighestOrder <= Components, "The highest order of the design must be less than the number of components." )
    DesignPoints = Inf
    try
        DesignPoints = 2 ^ ( Components ) - 2 ^ ( Components - HighestOrder )
        if isa(Order, Int)
            DesignPoints = 2 ^ ( HighestOrder ) - 1
        end
    catch
        @warn("Cannot estimate number of design points. May execute for the lifetime of the universe.")
    end

    PossiblePoints = zeros( Components )
    DesignSpace = []
    for order in Order
        PossiblePoints[ 1:order ] .= 1.0 / order
        push!( DesignSpace, unique( collect( Combinatorics.permutations(PossiblePoints, Components) ) ) )
    end
    DesignSpace = ( isa(Order, Int) ) ? DesignSpace[1] : reduce(vcat, DesignSpace)
    if !isinf(DesignPoints)
        println( DesignPoints )
        println(length(DesignSpace))
        @assert( DesignPoints == length(DesignSpace) )
    end
    return DesignSpace
end

# """
# The Sequential Generation of D-Optimum Experimental Designs. Wynn. The Annals of Mathematical Statistics. 41(5). 1970
# """
# function FedorovWynn(Components::Int, Points::Int)
#     X = rand( Points, Components )
#     while ( 1 / LinearAlgebra.cond(X' * X) ) < 1e-6
#         X = rand( Points, Components )
#     end
#     M = X' * X
#     Minv = Base.inv( M )
#     f(x) = sum(x)
# end


text:
struct GBFP <: MainPathAlgorithm end

struct FBMP <: MainPathAlgorithm end

struct GBFPResult{T<:Integer} <: MainPathResult
    edges::Vector{SimpleEdge{T}}
end

function vmax_outneighbors(g, v, w)
    nb = outneighbors(g, v)
    length(nb) == 0 && return nb
    wnb = w[nb]
    return nb[wnb .== maximum(wnb)]
end

function vmax_inneighbors(g, v, w)
    nb = inneighbors(g, v)
    length(nb) == 0 && return nb
    wnb = w[nb]
    return nb[wnb .== maximum(wnb)]
end

function mainpath(
    g::AbstractGraph{T},
    cutoff::Float64,
    ::GBFP;
    normalize=:global) where T <: Integer

    gkp = genetic_knowper(g, normalize=normalize)
    start = findall(gkp .>= cutoff)

    forwardstate = MainPathState(Vector{SimpleEdge{T}}())
    backwardstate = MainPathState(Vector{SimpleEdge{T}}())
    traverse_graph!(g, start, BreadthFirst(neighborfn = (g, v) -> vmax_outneighbors(g, v, gkp)), forwardstate)
    traverse_graph!(g, start, BreadthFirst(neighborfn = (g, v) -> vmax_inneighbors(g, v, gkp)), backwardstate)

    combined = vcat(forwardstate.edges, reverse.(backwardstate.edges)) |> unique
    return GBFPResult(combined)
end

function mainpath(
    g::AbstractGraph{T},
    start::Vector{T},
    ::GBFP;
    normalize=:global) where T <: Integer

    gkp = genetic_knowper(g, normalize=normalize)

    forwardstate = MainPathState(Vector{SimpleEdge{T}}())
    backwardstate = MainPathState(Vector{SimpleEdge{T}}())
    traverse_graph!(g, start, BreadthFirst(neighborfn = (g, v) -> vmax_outneighbors(g, v, gkp)), forwardstate)
    traverse_graph!(g, start, BreadthFirst(neighborfn = (g, v) -> vmax_inneighbors(g, v, gkp)), backwardstate)

    combined = vcat(forwardstate.edges, reverse.(backwardstate.edges)) |> unique
    return GBFPResult(combined)
end

function mainpath(
    g::AbstractGraph{T},
    start::Vector{T},
    vweights,
    ::FBMP) where T <: Integer

    forwardstate = MainPathState(Vector{SimpleEdge{T}}())
    backwardstate = MainPathState(Vector{SimpleEdge{T}}())
    traverse_graph!(g, start, BreadthFirst(neighborfn = (g, v) -> vmax_outneighbors(g, v, vweights)), forwardstate)
    traverse_graph!(g, start, BreadthFirst(neighborfn = (g, v) -> vmax_inneighbors(g, v, vweights)), backwardstate)

    combined = vcat(forwardstate.edges, reverse.(backwardstate.edges)) |> unique
    return GBFPResult(combined)
end


text:
"""
        mapr(t, a, b, c, d) -> Real

Maps the value `t` from the interval `[a, b]` to `[c, d]`
"""
function mapr(t::Real, a::Real, b::Real, c::Real, d::Real)
    return c + ((d - c) / (b - a)) * (t - a)
end
"""
        mapr(t, a, b) -> Real

Maps the value `t` from the interval `[0, 1]` to `[a, b]`
"""
mapr(t::Real, a::Real, b::Real) = mapr(t, 0, 1, a, b)
mapr(t::Real, (a, b)::Tuple, (c, d)::Tuple) = mapr(t, a, b, c, d)
mapr(t::Real, (a, b)::Tuple) = mapr(t, a, b)
mapr(t::Real, A::AbstractRange, B::AbstractRange) = mapr(t, extrema(A), extrema(B))
mapr(t::Real, A::AbstractRange) = mapr(t, extrema(A))


text:
include("../perfutil.jl")

include("eig.jl")
include("factorizations.jl")

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
using SanitationSystemMassFlow
using Base.Test

SSB = SanitationSystemMassFlow

@testset begin

    # -----------
    # techs

    transC_rel = Dict{String, Float64}("phosphor" => 5.0,
                                       "nitrogen" => 20.0,
                                        "water" => 100.0,
                                       "totalsolids" => 1000.0)

    # --
    transC_A = Dict{String, Dict{Product, Float64}}()
    transC_A["phosphor"] = Dict(Product("a1") => 0.5,
                                Product("a2") => 0.5,
                                Product("airloss") => 0.0,
                                Product("soilloss") => 0.0,
                                Product("waterloss") => 0.0)
    transC_A["nitrogen"] = Dict(Product("a1") => 0.5,
                                Product("a2") => 0.5,
                                Product("airloss") => 0.0,
                                Product("soilloss") => 0.0,
                                Product("waterloss") => 0.0)
    transC_A["water"] = Dict(Product("a1") => 0.5,
                             Product("a2") => 0.5,
                             Product("airloss") => 0.0,
                             Product("soilloss") => 0.0,
                             Product("waterloss") => 0.0)
    transC_A["totalsolids"] = Dict(Product("a1") => 0.5,
                                   Product("a2") => 0.5,
                                   Product("airloss") => 0.0,
                                   Product("soilloss") => 0.0,
                                   Product("waterloss") => 0.0)

    A = Tech(String[], ["a1", "a2"], "A", "group1", 0.5,
             transC_A,
             transC_rel)

    # --
    transC_B = Dict{String, Dict{Product, Float64}}()
    transC_B["phosphor"] = Dict(Product("b1") => 1.0,
                                Product("airloss") => 0.0,
                                Product("soilloss") => 0.0,
                                Product("waterloss") => 0.0)
    transC_B["nitrogen"] = Dict(Product("b1") => 1.0,
                                Product("airloss") => 0.0,
                                Product("soilloss") => 0.0,
                                Product("waterloss") => 0.0)
    transC_B["water"] = Dict(Product("b1") => 1.0,
                             Product("airloss") => 0.0,
                             Product("soilloss") => 0.0,
                             Product("waterloss") => 0.0)
    transC_B["totalsolids"] = Dict(Product("b1") => 1.0,
                                   Product("airloss") => 0.0,
                                   Product("soilloss") => 0.0,
                                   Product("waterloss") => 0.0)

    B = Tech(String[], ["b1"], "B", "group1", 0.5, transC_B, transC_rel)

    # --
    transC_C = Dict{String, Dict{Product, Float64}}()
    transC_C["phosphor"] = Dict(Product("c1") => 0.5,
                                Product("airloss") => 0.2,
                                Product("soilloss") => 0.3,
                                Product("waterloss") => 0.0)
    transC_C["nitrogen"] = Dict(Product("c1") => 1.0,
                                Product("airloss") => 0.0,
                                Product("soilloss") => 0.0,
                                Product("waterloss") => 0.0)
    transC_C["water"] = Dict(Product("c1") => 0.6,
                             Product("airloss") => 0.0,
                             Product("soilloss") => 0.0,
                             Product("waterloss") => 0.4)
    transC_C["totalsolids"] = Dict(Product("c1") => 1.0,
                                   Product("airloss") => 0.0,
                                   Product("soilloss") => 0.0,
                                   Product("waterloss") => 0.0)
    C = Tech(["a1"], ["c1"], "C", "group1", 0.5,
             transC_C, transC_rel)

    # --

    transC_D = Dict{String, Dict{Product, Float64}}()
    transC_D["phosphor"] = Dict(Product("d1") => 0.1,
                                Product("d2") => 0.1,
                                Product("d3") => 0.8,
                                Product("airloss") => 0.0,
                                Product("soilloss") => 0.0,
                                Product("waterloss") => 0.0)
    transC_D["nitrogen"] = Dict(Product("d1") => 0.3,
                                Product("d2") => 0.3,
                                Product("d3") => 0.1,
                                Product("airloss") => 0.0,
                                Product("soilloss") => 0.3,
                                Product("waterloss") => 0.0)
    transC_D["water"] = Dict(Product("d1") => 0.0,
                             Product("d2") => 0.0,
                             Product("d3") => 1.0,
                             Product("airloss") => 0.0,
                             Product("soilloss") => 0.0,
                             Product("waterloss") => 0.0)
    transC_D["totalsolids"] = Dict(Product("d1") => 0.5,
                                   Product("d2") => 0.5,
                                   Product("d3") => 0.0,
                                   Product("airloss") => 0.0,
                                   Product("soilloss") => 0.0,
                                   Product("waterloss") => 0.0)

D = Tech(["a2", "b1"], String["d1", "d2", "d3"], "D", "group1", 0.5,
         transC_D, transC_rel)

# --
transC_E = Dict{String, Dict{Product, Float64}}()
transC_E["phosphor"] = Dict(Product("recovered") => 1.0,
                            Product("airloss") => 0.0,
                            Product("soilloss") => 0.0,
                            Product("waterloss") => 0.0)
transC_E["nitrogen"] = Dict(Product("recovered") => 1.0,
                            Product("airloss") => 0.0,
                            Product("soilloss") => 0.0,
                            Product("waterloss") => 0.0)
transC_E["water"] = Dict(Product("recovered") => 1.0,
                         Product("airloss") => 0.0,
                         Product("soilloss") => 0.0,
                         Product("waterloss") => 0.0)
transC_E["totalsolids"] = Dict(Product("recovered") => 0.8,
                               Product("airloss") => 0.0,
                               Product("soilloss") => 0.0,
                               Product("waterloss") => 0.2)

E = Tech(["c1", "d1"], String[], "E", "group1", 0.5, transC_E, transC_rel)

# --

transC_F = Dict{String, Dict{Product, Float64}}()
transC_F["phosphor"] = Dict(Product("recovered") => 1.0,
                            Product("airloss") => 0.0,
                            Product("soilloss") => 0.0,
                            Product("waterloss") => 0.0)
transC_F["nitrogen"] = Dict(Product("recovered") => 1.0,
                            Product("airloss") => 0.0,
                            Product("soilloss") => 0.0,
                            Product("waterloss") => 0.0)
transC_F["water"] = Dict(Product("recovered") => 0.5,
                         Product("airloss") => 0.5,
                         Product("soilloss") => 0.0,
                         Product("waterloss") => 0.0)
transC_F["totalsolids"] = Dict(Product("recovered") => 1.0,
                               Product("airloss") => 0.0,
                               Product("soilloss") => 0.0,
                               Product("waterloss") => 0.0)

F = Tech(["d2", "d3"], String[], "F", "group1", 0.5, transC_F, transC_rel)

# --

transC_G = Dict{String, Dict{Product, Float64}}()
transC_G["phosphor"] = Dict(Product("recovered") => 0.0,
                            Product("airloss") => 0.0,
                            Product("soilloss") => 0.0,
                            Product("waterloss") => 1.0)
transC_G["nitrogen"] = Dict(Product("recovered") => 0.0,
                            Product("airloss") => 0.0,
                            Product("soilloss") => 0.0,
                            Product("waterloss") => 1.0)
transC_G["water"] = Dict(Product("recovered") => 0.0,
                         Product("airloss") => 0.5,
                         Product("soilloss") => 0.5,
                         Product("waterloss") => 0.0)
transC_G["totalsolids"] = Dict(Product("recovered") => 0.0,
                               Product("airloss") => 1.0,
                               Product("soilloss") => 0.0,
                               Product("waterloss") => 0.0)

G = Tech(["a1", "a2", "b1"], String[], "G", "group1", 0.5, transC_G, transC_rel)

# --
# this shoud fail
transC_error = Dict{String, Dict{Product, Float64}}()
transC_error["phosphor"] = Dict(Product("recovered") => 0.0,
                                Product("airloss") => 0.0,
                                Product("soilloss") => 1.0,
                                Product("waterloss") => 1.0)
transC_error["nitrogen"] = Dict(Product("recovered") => 0.0,
                                Product("airloss") => 0.0,
                                Product("soilloss") => 0.0,
                                Product("waterloss") => 1.0)
transC_error["water"] = Dict(Product("recovered") => 0.0,
                             Product("airloss") => 0.5,
                             Product("soilloss") => 0.5,
                             Product("waterloss") => 0.0)
transC_error["totalsolids"] = Dict(Product("recovered") => 0.0,
                                   Product("airloss") => 1.0,
                                   Product("soilloss") => 0.0,
                                   Product("waterloss") => 0.0)

@test_throws ErrorException Tech(["a1", "a2", "b1"], String[], "G", "group1", 0.5, transC_error, transC_rel)

# --

@test length(SSB.get_inputs([A])) == 0
@test length(SSB.get_inputs([A,D])) == 2
@test length(SSB.get_inputs([A,D,E])) == 4

@test length(SSB.get_outputs([A])) == 2
@test length(SSB.get_outputs([A,D])) == 5
@test length(SSB.get_outputs([A,D,E])) == 5

# -----------
# System

# system without  triangle
allSys = build_all_systems([A, B], [C, D, E, F, G])
@test length(allSys) == 2

# -----------
# Mass flow

# define input masses for each source

M_in = Dict("A" => Dict("phosphor" => 600,
                        "nitrogen" => 400,
                        "water" => 260,
                        "totalsolids" => 90),
            "B" => Dict("phosphor" => 60,
                        "nitrogen" => 40,
                        "water" => 26,
                        "totalsolids" => 9))

# test mass balances
for sys in allSys
    m_outs = massflow(sys, M_in)
    @test isapprox(entered(M_in, sys) - recovered(m_outs) - sum(lost(m_outs),2),
                   [0.0, 0.0, 0.0, 0.0],
                   atol=1e-12)

    m_outs = massflow(sys, M_in, MC=true)
    @test isapprox(entered(M_in, sys) - recovered(m_outs) - sum(lost(m_outs),2),
                   [0.0, 0.0, 0.0, 0.0], atol=1e-12)

    m_outs = massflow(sys, M_in, MC=true, scale_reliability=100)
    @test isapprox(entered(M_in, sys) - recovered(m_outs) - sum(lost(m_outs),2),
                   [0.0, 0.0, 0.0, 0.0], atol=1e-12)

end

end


text:
#=
Copyright 2020 Ali Mert Ceylan

Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
=#

import Base: broadcasted, abs2, sum, inv
import Statistics: mean
import LinearAlgebra: adjoint, transpose, det

function flat(ra)
    if(typeof(ra) <: Pair)
        return [ra]
    end

    ta = []
    for a in ra
        if(!(typeof(a) <: Number))
            append!(ta, flat(a))
        end
    end
    return ta
end

params(xs) = Dict(flat(xs))

macro differentiate(e)
    args = e.args

    local varname = string(args[1])
    local var = args[1]
    local val = args[2]

    esc(:($var = BD{typeof($val)}($val, Symbol("d$(string($varname))"), )))
end

function sum(B::BD{T}; kwargs...) where T <: AbstractArray
    s = sum([b for b in B.f[1]], kwargs...)
    return BD{typeof(s)}(s, (dy)->(B.f[2](dy .* ones(size(B)))))
end

# @primitive mean(x;d...),dy  (dy.*one.(x)./(length(x)÷length(dy)))
function mean(bd::BD{T}; kwargs...) where T <: AbstractArray
    s = mean(bd.f[1])
    return BD{typeof(s)}(s, (dy)->(bd.f[2](dy .* ones(size(bd.f[1])) ./ length(bd.f[1]))))# dy.*ones(size(bd.f[1]))./length(bd.f[1])
end

abs2(b::BD{T}) where T <: Number = BD{T}(abs2(b.f[1]), (dy)->(b.f[2](dy * 2b.f[1])))
broadcasted(::typeof(abs2), b::BD{T}) where T <: AbstractArray = BD{T}(abs2.(b.f[1]), (dy)->(b.f[2](dy .* 2b.f[1])))

transpose(b::BD{T}) where T <: AbstractArray = BD{T}(transpose(b.f[1]), (dy)->(transpose(dy)))
adjoint(b::BD{T}) where T <: AbstractArray = BD{T}(adjoint(b.f[1]), (dy)->(adjoint(dy)))

# inv(x),dy,y  (yt=y'; -yt*dy*yt) # addtest(:inv, rand(2,2)) # ref https://people.maths.ox.ac.uk/gilesm/files/NA-08-01.pdf
# inv(b::BD{T}) where T <: AbstractArray = BD{T}(inv(b.f[1]), (dy)->(yt = inv(b.f[1])'; -yt * dy * yt))

# det(x),dy,y  dy*y*inv(x)'
det(b::BD{T}) where T <: AbstractArray = BD{T}(det(b.f[1]), (dy)->(dy * det(b.f[1]) * inv(b.f[1])'))

# TODO diag(x),dy,y   diagm(0=>dy)  # alternative: Diagonal(dy)
# TODO diag(x,i),dy,y   diagm(i=>dy)  # warning: these only works for square matrices
# TODO dot(x1, x2),dy,y  dy*x2  dy*x1  # addtestN(:dot, rand(3,2), rand(3,2))

# TODO @zerograd axes(x,i...)
# TODO @primitive maximum(x;d...),dy,y  (dy.*(y.==x))
# TODO @primitive maximum(f::typeof(abs),x;d...),dy,y  nothing  (dy.*(y.==abs.(x)).*sign.(x))
# TODO @primitive minimum(x;d...),dy,y  (dy.*(y.==x))
# TODO @primitive minimum(f::typeof(abs),x;d...),dy,y  nothing  (dy.*(y.==abs.(x)).*sign.(x))


text:
using Test
using Example2018

# Polynomials form a ring
# 1. addition is associatice
# 2. additive identity
# 3. additice inverse
# 4. addition commutes
# 5. multiplication is ...

@testset "Polynomial form a ring" begin

	z = Polynomial([0])

	for n in 1:100

	deg = rand(0:100)
		p = Polynomial(rand(-1000000:1000000, deg+1))

		@test p + z == p
		@test z + p == p
		@test p + (-p) == z
	end

end


text:
n = parse(Int, readline())
for i in 1:n
    global x = i
    if x % 3 == 0
        print(" ", i)
    else
        truth = true
        while truth
            if x % 10 == 3
                print(" ", i)
                break
            end
            x = div(x, 10)
            truth = (x != 0)
        end
    end
end
println()