text:
## RandomShuffle ##

"""
    RandomShuffle

type with no fields representing a completely random shuffle, which is guaranteed to be
equivalent to [`Random.shuffle`](https://docs.julialang.org/en/v1/stdlib/Random/#Random.shuffle)).

This algorithm is set as the default. See [`DEFAULTS`](@ref) and has the singleton instance
[`randshuffle`](@ref).

# Examples
```jldoctest
julia> mt1 = MersenneTwister(1234); mt2 = MersenneTwister(1234);

julia> shuffle(mt1, 1:10, randshuffle) == Random.shuffle(mt2, 1:10)
true

```
"""
struct RandomShuffle <: AbstractRandomShuffle end

"""
    randshuffle

The singleton instance of type [`RandomShuffle`](@ref)
"""
const randshuffle = RandomShuffle()

shuffle!(r::AbstractRNG, c::AbstractArray, s::RandomShuffle) = Random.shuffle!(r, c)

## GilbertShannonReeds ##

"""
    GilbertShannonReeds

type with no fields representing the [Gilbert-Shannon-Reeds](https://en.wikipedia.org/wiki/Gilbert–Shannon–Reeds_model)
model of card shuffling. See singleton instance [`gsrshuffle`](@ref).

An in-place [`shuffle!`](@ref) is not implemented for this algorithm. However, an in-place
[`nshuffle!`](@ref) is implemented.

# Examples
```jldoctest
julia> mt = MersenneTwister(1234);

julia> shuffle(mt, 1:7, gsrshuffle)
7-element Vector{Int64}:
 5
 6
 1
 2
 7
 3
 4

```
"""
struct GilbertShannonReeds <: AbstractRandomShuffle end

"""
    gsrshuffle

The singleton instance of type [`GilbertShannonReeds`](@ref).
"""
const gsrshuffle = GilbertShannonReeds()

function shuffle(
    r::AbstractRNG,
    c::AbstractArray,
    s::GilbertShannonReeds,
    out::AbstractArray = similar(c),
)
    length(out) == length(c) || throw(ArgumentError("array lengths c and out must equal"))

    # flip a coin to determine from which pile each item will come
    flips = rand(r, Bool, length(c))

    # positions in the collection for the first pile and second pile
    pos1 = 0
    pos2 = count(flips)

    # Get items from respective pile
    for (i, frmfirst) in enumerate(flips)
        @inbounds out[i] = c[frmfirst ? pos1 += 1 : pos2 += 1]
    end

    return out
end

function nshuffle!(r::AbstractRNG, c::AbstractArray, n::Integer, s::GilbertShannonReeds)
    tmp = similar(c)

    for i in 1:fld(n, 2)
        shuffle(r, c, s, tmp)
        shuffle(r, tmp, s, c)
    end

    if isodd(n)
        return shuffle(r, c, s, tmp)
    else
        return c
    end
end


text: Dict{String,Any}("thevoid" => Dict{String,Any}("value" => Any[Dict{String,Any}("value" => Any[Dict{String,Any}("value" => Any[Dict{String,Any}("value" => Any[Dict{String,Any}("value" => Any[],"type" => "array")],"type" => "array")],"type" => "array")],"type" => "array")],"type" => "array"))

ext: jl

max_stars_repo_name: vanillajonathan/julia

max_issues_repo_name: vanillajonathan/julia

max_forks_repo_name: vanillajonathan/julia


text:
using HTTP, Sockets, JSON3, Logging
using HTTP.URIs: queryparams, URI

const ROUTER = HTTP.Router()

function serveWelcome(req::HTTP.Request)
    return HTTP.Response(200, "Welcome to Genie! This site is served by Julia $(Base.VERSION) with $(Threads.nthreads()) threads.")
end

function serveRandom(req::HTTP.Request)
    num = parse(Int, get(queryparams(URI(req.target)), "num", "1"))

    return HTTP.Response(200, JSON3.write(rand(Int, num)))
end

HTTP.@register(ROUTER, "GET", "/", serveWelcome)
HTTP.@register(ROUTER, "GET", "/random", serveRandom)

Logging.global_logger(NullLogger())
HTTP.serve(ROUTER, Sockets.localhost, 8000)

max_stars_repo_name: essenciary/performance-tests

max_issues_repo_name: essenciary/performance-tests

max_forks_repo_name: essenciary/performance-tests


text:
mutable struct Object
    vertices::Vector{Point3D}
    faces::Array{Array{Int64,1},1}
    labels::Array{Int64,1}
    name::String
end

function Base.size(o::Object)
    length(o.vertices)
end

function Base.length(o::Object)
    length(o.vertices)
end

Base.lastindex(o::Object) = length(o)

Base.broadcastable(o::Object) = Ref(o)

"""
    make(primitive, name="unnamed")

`primitive` contains two arrays, an array of 3D points, and an array of faces,
where each face is a list of vertex numbers.

Returns a Object.

Example

    make(Cube, "cube")

returns an Object object containing an array of vertices, an array of faces,
and an array of labels.

```
@draw begin
    helloworld()
    tol = 0.01
    a = []
    sethue("black")
    for t in -2pi:tol:2pi
        push!(a, Point3D((50 + cos(5t)) * cos(3t), (50 + cos(5t)) * sin(2t), sin(5t)))
    end
    Knot = make((a, []), "knot")
    pin(Knot, gfunction = (args...) -> poly(args[1], :stroke))
end

The default gfunction expects faces - if there aren't any, use a gfunction that draws vertices.
```
"""
function make(vf, name="unnamed")
    # don't redefine when passed an array
    vertices = deepcopy(vf[1])
    faces    = deepcopy(vf[2])
    labels   = collect(1:length(faces))
    return Object(vertices, faces, labels, name)
end

include("../data/objects.jl")

# don't load all these now!
# include("moreobjects.jl")

"""
    objecttopoly(o::Object)

Return a tuple:

- an array of 2D points representing the vertices of `o`

- an array of 2D polygons representing the faces of `o`

## Example

This example draws the faces of a cube in colors, and marks
the vertices in black.

```
using Luxor, Thebes

@draw begin
    helloworld()
    o = make(Cube)
    scaleby!(o, 100, 100, 100)
    vs, fs = objecttopoly(o)
    setopacity(0.4)
    sethue("black")
    for face in fs
        randomhue()
        poly(face, :fill)
    end
    sethue("black")
    circle.(vs, 3, :fill)
end
```
"""
function objecttopoly(o::Object)
    vertices2D = Point[]
    for v in o.vertices
        r = project(v)
        if r != nothing
            push!(vertices2D, r)
        else
            push!(vertices2D, Point(NaN, NaN))
        end
    end
    facepolys = Array[]
    if length(o.faces) > 0
        for n in 1:length(o.faces)
            push!(facepolys, vertices2D[o.faces[n]])
        end
    end
    filter!(f -> !isnan(f.x) && !isnan(f.y), vertices2D)
    return (vertices2D, facepolys)
end

"""
    sortfaces!(o::Object;
        eyepoint::Point3D=eyepoint())

Find the averages of the z values of the faces in Object, and sort the faces
of o so that the faces are in order of nearest (highest) z relative to eyepoint...

or something like that ? not sure how this works
"""
function sortfaces!(o::Object;
        eyepoint::Point3D=eyepoint())
    avgs = Float64[]
    for f in o.faces
        vs = o.vertices[f]
        s = 0.0
        for v in vs
            s += distance(v, eyepoint)
        end
        avg = s/length(unique(vs))
        push!(avgs, avg)
    end
    neworder = reverse(sortperm(avgs))
    o.faces = o.faces[neworder]
    o.labels = o.labels[neworder]
    return o
end

sortfaces!(o::Array{Object, 1}; kwargs...) =
   map(sortfaces!, o)

function hiddensurface(o::Object)
   if !isempty(o.faces)
       sortfaces!(o)
       @layer begin
           for (n, face) in enumerate(o.faces)
               @layer begin
                   vertices = o.vertices[face]
                   sn = surfacenormal(vertices)
                   ang = anglebetweenvectors(sn, eyepoint())
                   setgrey(rescale(ang, 0, π, 1, 0))
                   pin(vertices, gfunction = (p3, p2) ->
                    begin
                       poly(p2, :fill)
                       sethue("white")
                       poly(p2, :stroke, close=true)
                    end)
               end
           end
       end
   end
end

function wireframe(o::Object)
   if !isempty(o.faces)
       sortfaces!(o)
       @layer begin
           for (n, face) in enumerate(o.faces)
               @layer begin
                   vertices = o.vertices[face]
                   pin(vertices, gfunction = (p3, p2) ->
                    begin
                       poly(p2, :stroke, close=true)
                    end)
               end
           end
       end
   end
end

"""
    pin(o::Object;
        gfunction=(o) -> hiddensurface(o))

Draw a rendering of an object.

The default rendering function is `hiddensurface()`.

You can also use the built-in `wireframe()` rendering function.

## Examples

```
@draw begin
    o = make(Cube)
    axes3D(200)
    scaleby!(o, 200, 200, 200)
    eyepoint(250, 270, 300)
    pin(o) # use an attempted hiddensurface rendering

    o = make(Tetrahedron)
    axes3D(200)
    scaleby!(o, 200, 200, 200)
    eyepoint(250, 270, 300)
    pin(o, gfunction=wireframe) # use a wireframe rendering
end
```
## More help

You could write your own rendering function to draw objects.

```
function a_rendering_function(o::Object)
   if !isempty(o.faces)
       sortfaces!(o)
       @layer begin
           for (n, face) in enumerate(o.faces)
               @layer begin
                   vertices = o.vertices[face]
                   sn = surfacenormal(vertices)
                   ang = anglebetweenvectors(sn, eyepoint())
                   setgrey(rescale(ang, 0, π, 1, 0))
                   pin(vertices, gfunction = (p3, p2) ->
                    begin
                       poly(p2, :fill)
                       sethue("white")
                       poly(p2, :stroke, close=true)
                    end)
               end
           end
       end
   end
end
```
"""
function pin(o::Object;
    gfunction = (o) -> hiddensurface(o))
    gfunction(o)
end

"""
    moveby!(o::Object, x, y, z)
    moveby!(o::Object, pt::Point3D)

Set the position of object to Point3D(x, y, z).
"""
function moveby!(o::Object, x, y, z)
   for n in 1:length(o.vertices)
       nv = o.vertices[n]
       o.vertices[n] = Point3D(nv.x + x, nv.y + y, nv.z + z)
   end
   return o
end

"""
    moveby(o::Object, x, y, z)
    moveby(o::Object, pt::Point3D)

Set the position of a copy of the object to Point3D(x, y, z).
"""
function moveby(o::Object, x, y, z)
   ocopy = deepcopy(o)
   return moveby!(ocopy, x, y, z)
end

moveby(o::Object, pt::Point3D) = moveby(o::Object, pt.x, pt.y, pt.z)
moveby!(o::Object, pt::Point3D) = moveby!(o::Object, pt.x, pt.y, pt.z)

"""
    scaleby!(o::Object, x, y, z)

Scale object by x in x, y in y, and z in z.
"""
function scaleby!(o::Object, x, y, z)
   for n in 1:length(o.vertices)
       nv = o.vertices[n]
       o.vertices[n] = Point3D(nv.x * x, nv.y * y, nv.z * z)
   end
   return o
end

"""
    face(o::Object, n)
"""
function face(o::Object, n)
   facepoints = Point3D[]
   if length(o.faces) > 0
       for i in o.faces[n]
           push!(facepoints, o.vertices[i])
       end
   end
   return facepoints
end

"""
    rotateby!(o::Object, r::Rotation)
    rotateby!(o::Object, angleX, angleY, angleZ)

Rotate an object through rotation `r`, or around the x, y,
and/or z axis by `angleX`, `angleY`, `angleZ`.
"""
function rotateby!(o::Object, angleX, angleY, angleZ)
    for n in 1:length(o.vertices)
        o.vertices[n] = RotXYZ(angleX, angleY, angleZ) * o.vertices[n]
    end
    return o
end

function rotateby!(o::Object, r::Rotation)
    for n in 1:length(o.vertices)
        o.vertices[n] = r * o.vertices[n]
    end
    return o
end

"""
    rotateby(o::Object, angleX, angleY, angleZ)

Rotate a copy of the object by angleX, angleY, angleZ.
"""
function rotateby(o::Object, angleX, angleY, angleZ)
    ocopy = deepcopy(o)
    return rotateby!(ocopy, angleX, angleY, angleZ)
end

function rotateby(o::Object, r::Rotation)
    ocopy = deepcopy(o)
    return rotateby!(ocopy, r)
end

"""
    rotateby!(o::Object, pt::Point3D, angleX, angleY, angleZ)
    rotateby!(o::Object, pt::Point3D, r::Rotation=RotXYZ(0, 0, 0))

Rotate an object around a point by rotation r, or angleX, angleY, angleZ.
"""
function rotateby!(o::Object, pt::Point3D, angleX, angleY, angleZ)
    for n in 1:length(o.vertices)
        v = o.vertices[n] - pt
        v = rotateX(v, angleX)
        v = rotateY(v, angleY)
        v = rotateZ(v, angleZ)
        o.vertices[n] = v + pt
    end
    return o
end

function rotateby!(o::Object, pt::Point3D, r::Rotation=RotXYZ(0, 0, 0))
    for n in 1:length(o.vertices)
        v = o.vertices[n] * r
        o.vertices[n] = v
    end
    return o
end

"""
    rotateby(o::Object, pt::Point3D, angleX, angleY, angleZ)
    rotateby(o::Object, pt::Point3D, r::Rotation=RotXYZ(0, 0, 0))

Rotate a copy of the object around a point by rotation r, or angleX, angleY, angleZ.
"""
function rotateby(o::Object, pt::Point3D, angleX, angleY, angleZ)
    ocopy = deepcopy(o)
    return rotateby!(ocopy, pt, angleX, angleY, angleZ)
end

function rotateby(o::Object, pt::Point3D, r::Rotation=RotXYZ(0, 0, 0))
    ocopy = deepcopy(o)
    return rotateby!(ocopy, pt, r)
end


text:
"""
    function selection_classic(fob::Function, trial::AbstractMatrix, population::AbstractMatrix)

"""
function selection_classic(fob::Function, trial::AbstractMatrix, population::AbstractMatrix)
    #error checks
    if size(trial) != size(population)
        throw(ErrorException("Size of trial matrix $(size(trial)) must be the same as the population matrix $(size(population))"))
    end

    population_size = size(population,1)
    selected = similar(population)

    for individue in 1:population_size
        @inbounds if fob(trial[individue,:]...) < fob(population[individue,:]...)
            @inbounds selected[individue,:] = trial[individue,:]
        else
            @inbounds selected[individue,:] = population[individue,:]
        end
    end

    return selected
end


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct Area <: AggregationTopology
        name::String
        internal::InfrastructureSystemsInternal
    end

A collection of buses.

# Arguments
- `name::String`: area name
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct Area <: AggregationTopology
    "area name"
    name::String
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function Area(name, )
    Area(name, InfrastructureSystemsInternal(), )
end

function Area(; name, )
    Area(name, )
end

# Constructor for demo purposes; non-functional.
function Area(::Nothing)
    Area(;
        name="init",
    )
end

InfrastructureSystems.get_name(value::Area) = value.name
"""Get Area internal."""
get_internal(value::Area) = value.internal


text:
"""
simulation & machine learning of diffusion advection PDE
"""

using IterTools
using Plots
using Random
using Flux
using LinearAlgebra
Random.seed!(1)
using Main.EquivariantOperators # end users should omit Main.

# make grid
n = 2
dx = 0.02
cell = dx * Matrix(I, n, n)
rmax = 1.0
grid = Grid(cell, rmax)
sz = size(grid)
▽ = Op(:▽, cell)

# put small blob as IC
u0 = zeros(sz..., 1)
put!(u0, grid, [0.0, 0.0], [1.0])
g = Op(:Gaussian, cell; σ = dx)
u0 = g(u0)

# diffusion advection
p1 = [0.1]
p2 = [0.5, 0.4]
vf = cat(p2[1] * ones(sz), p2[2] * ones(sz), dims = 3)
f(u, p, t) = p1[1] * (▽ ⋅ ▽(u)) - vf ⊗ ▽(u)

# simulate PDE
using DifferentialEquations
tspan = (0.0, 1.0)
prob = ODEProblem(f, u0, tspan)
sol = solve(prob, Tsit5(), reltol = 1e-3, abstol = 1e-3)

# plot
using Plots
gr()
anim = Animation()

t = 0:0.02:1
for t in t
heatmap(sol(t)[:, :, 1], clim=(0,10))
frame(anim)
end
gif(anim, "f.gif", fps = 10)

##
data = [(sol(t), f(sol(t), 0, 0)) for t in t]
# op = Op(Radfunc(),-1e-6, 2dx, cell)
# ps=Flux.params(op)
p1_ = ones(1)
p2_ = ones(2)
ps = Flux.params(p1_, p2_)

function loss(u, du)
    vf_ = cat(p2_[1] * ones(sz), p2_[2] * ones(sz), dims = 3)
    duhat = p1_[1] .* (▽ ⋅ ▽(u)) - vf_ ⊗ ▽(u)
    @show l = nae(duhat, du)
end

loss(data[1]...)
opt = ADAM(0.1)

Flux.@epochs 20 Flux.train!(loss, ps, data, opt)

@show p1_, p2_
# heatmap(op.kernel[:, :, 1])


text:
module OMEinsum

@doc raw"
   greet()
say hello to the world!
"
greet() = print("Hello World!")

end # module


text:
#=
Created on Saturday 3 August 2019
Last update: -

@author: Michiel Stock
michielfmstock@gmail.com

Vec trick: multiplying vectors with Kronecker systems.
=#

import Base.ReshapedArray

vecmulR!(X,N,V,M) = X .= N * (V * transpose(M))
vecmulL!(X,N,V,M) = X .= (N * V) * transpose(M)
vecmul!(X,N,V,M) = X .= N * V * transpose(M)

vectrick_reshape(v::AbstractVector, d::Int, b::Int) = reshape(v, d, b)
function vectrick_reshape(v::ReshapedArray{<:Any, 1}, d::Int, b::Int)
    return size(v.parent) == (d, b) ? v.parent : reshape(v, d, b)
end

"""
    mul!(x::AbstractVector, K::AbstractKroneckerProduct, v::AbstractVector)

Calculates the vector-matrix multiplication `K * v` and stores the result in
`x`, overwriting its existing value.
"""
function mul!(x::AbstractVector, K::AbstractKroneckerProduct, v::AbstractVector)
    M, N = getmatrices(K)
    a, b = size(M)
    c, d = size(N)
    e = length(v)
    f = length(x)
    f == a * c || throw(DimensionMismatch(
        "Dimension missmatch between kronecker system and result placeholder"))
    e == b * d || throw(DimensionMismatch(
        "Dimension missmatch between kronecker system and vector"))
    V = vectrick_reshape(v, d, b)
    X = vectrick_reshape(x, c, a)
    if b * c * (a + d) < a * d * (b + c)
        mul!(X, N, V * transpose(M))
    else
        mul!(X, N * V, transpose(M))
    end
    return x
end

function Base.:*(K::AbstractKroneckerProduct, v::AbstractVector)
    return mul!(Vector{promote_type(eltype(v), eltype(K))}(undef, first(size(K))), K, v)
end

reshape_cols(x, sizes...) = reshape(x, sizes...)
reshape_rows(x, sizes...) = transpose(reshape(collect(transpose(x)), reverse(sizes)...))

kron_id_a(a, x) = reshape_cols(a * reshape_cols(x, size(a)[2], :), :, size(x)[2])
kron_a_id(a, x) = reshape_rows(a * reshape_rows(x, size(a)[2], :), :, size(x)[2])

function kron_a_b(A, B, x)
    a, b = size(A)
    c, d = size(B)
    if b * c * (a + d) < a * d * (b + c)
        return kron_a_id(A, kron_id_a(B, x))
    else
        return kron_id_a(B, kron_a_id(A, x))
    end
end

function Base.:*(A::AbstractKroneckerProduct, B::StridedMatrix)
    size(A, 2) != size(B, 1) && throw(DimensionMismatch("size(A, 2) != size(B, 1)"))
    return kron_a_b(getmatrices(A)..., B)
end

function Base.:*(A::KroneckerProduct{<:Any, <:Eye, <:AbstractMatrix}, B::StridedMatrix)
    size(A, 2) != size(B, 1) && throw(DimensionMismatch("size(A, 2) != size(B, 1)"))
    return kron_id_a(A.B, B)
end

function Base.:*(A::KroneckerProduct{<:Any, <:AbstractMatrix, <:Eye}, B::StridedMatrix)
    size(A, 2) != size(B, 1) && throw(DimensionMismatch("size(A, 2) != size(B, 1)"))
    return kron_a_id(A.A, B)
end

function Base.sum(K::AbstractKroneckerProduct; dims::Union{Nothing,Int}=nothing)
    A, B = getmatrices(K)
    if dims == nothing
        s = zero(eltype(K))
        sumB = sum(B)
        return sum(sumB * A)
    else
        return kronecker(sum(A, dims=dims), sum(B, dims=dims))
    end
end


text:
include("examples/mlp/MLPDemo.jl")
using AF, MLPDemo

af = ArrayFire{CUDA}()

runDemo(af)

max_stars_repo_name: unbornchikken/julia-ml-proto

max_issues_repo_name: unbornchikken/julia-ml-proto

max_forks_repo_name: unbornchikken/julia-ml-proto


text:
using Luxor
using Colors

function create_logo(text_logo::String, font_family::String; name::Union{String,Symbol} = :svg, width::Int=1250, height::Int=800)

    Drawing(width, height, name)

    background("transparent")
    origin()
    fontsize(300)                             # big fontsize to use for clipping
    fontface(font_family)
    w, h = textextents(text_logo)[3:4]

    translate(-(width / 2) + 50, -(height / 2) + 1.1h)

    textpath(text_logo)                             # make text into a path
    setline(10)                          # fill but keep
    clip()                                    # and use for clipping region

    fontface("Monaco")
    fontsize(8)
    namelist = map(x -> string(x), names(Base)) # get list of function names in Base.

    let
        x = -20
        y = -h
        while y < height
            sethue(rand(7:10) / 10, rand(7:10) / 10, rand(7:10) / 10)
            s = namelist[rand(1:end)]
            text(s, x, y)
            se = textextents(s)
            x += se[5]                            # move to the right
            if x > w
                x = -20                            # next row
                y += 10
            end
        end
    end

    finish()
    preview()
end

file_name = "assets/output.svg"
logo = "TEST"
font_family = "MathJax_Script"
create_logo(logo, font_family)
create_logo(logo, font_family, name = file_name)

create_logo("MB", "MathJax_Script")
create_logo("MB", "MathJax_Script", name = "assets/mb.svg", height = 300, width = 800)

# Create M and B in separate files
files = ["M", "B"]
for letter in files
    create_logo("assets/tmp/$(letter).svg", letter, font_family, width=500, height=300)
end

# Join them with a separation
function join_images(files)
    Drawing(1000, 300, :svg)
    origin()
    image = readsvg("assets/tmp/$(files[1]).svg")
    for (index, file) in enumerate(files[2:end])
        println("Current letter: ", file)
        img = readsvg("assets/tmp/$(file).svg")

        w = img.width
        h = img.height

        placeimage(image, O)
        # setmode("saturate")
        # translate(300, 0)
    # placeimage(img, -w / 2, -h / 2, .5)
    end
    finish()
    preview()
end

join_images(files)

max_stars_repo_name: mattborghi/logo

max_issues_repo_name: mattborghi/logo

max_forks_repo_name: mattborghi/logo


text:
"""

```
load(filename::String)
```

Function to load the json saved results file into the `Result` data structure.

"""
function load(filename::String)
    f = open(filename, "r")
    R = JSON3.read(f, MutableResult)
    # Need to reshape the solute and solvent atom contributions, because
    # the data is read in a single column
    solute_atom = reshape(R.solute_atom, R.nbins, :)
    solvent_atom = reshape(R.solvent_atom, R.nbins, :)
    R.solute_atom = solute_atom
    R.solvent_atom = solvent_atom
    return Result([getfield(R, field) for field in fieldnames(Result)]...)
end


text:
export convex_hull

function default_convex_hull_algorithm(points)
    if length(points) != 0 && length(first(points)) == 2
        return "monotone_chain"
    else
        return nothing
    end
end

"""
    convex_hull(X::LazySet{N}, Y::LazySet{N}; [algorithm]=nothing,
                [backend]=nothing, [solver]=nothing) where {N<:Real}

Compute the convex hull of the given convex sets.

### Input

- `X`         -- convex set
- `Y`         -- convex set
- `algorithm` -- (optional, default: `nothing`) the convex-hull algorithm
- `backend`   -- (optional, default: `nothing`) backend for polyhedral
                 computations (used for higher-dimensional sets)
- `solver`    -- (optional, default: `nothing`) the linear-programming solver
                 used in the backend

### Output

If the input sets are one-dimensional, the result is an `Interval`.
If the input sets are two-dimensional, the result is a `VPolygon`.
Otherwise the result is a `VPolytope`.

### Algorithm

One-dimensional sets are resolved by using `overapproximate` with an `Interval`
(which is exact).
For higher-dimensional sets, we compute the vertices of both `X` and `Y` using
`vertices_list` and then compute the convex hull of the union of those vertices.
"""
function convex_hull(X::LazySet{N}, Y::LazySet{N};
                     algorithm=nothing,
                     backend=nothing,
                     solver=nothing
                    ) where {N<:Real}
    n = dim(X)
    @assert n == dim(Y) "the convex hull requires two sets of the same " *
                        "dimension, but the sets had dimension $n and $(dim(Y))"

    chull = convex_hull!([vertices_list(X); vertices_list(Y)];
                         algorithm=algorithm, backend=backend, solver=solver)
    m = length(chull)
    if m == 0
        return EmptySet{N}(n)
    elseif m == 1
            return Singleton(chull[1])
    elseif n == 1
        @assert m == 2
        low = chull[1][1]
        high = chull[2][1]
        if low > high
            tmp = low
            low = high
            high = tmp
        end
        return Interval(low, high)
    elseif n == 2
        return VPolygon(chull)
    end
    return VPolytope(chull)
end

"""
    convex_hull(points::Vector{VN};
                [algorithm]=nothing,
                [backend]=nothing,
                [solver]=nothing
                ) where {N<:Real, VN<:AbstractVector{N}}

Compute the convex hull of the given points.

### Input

- `points`    -- list of vectors
- `algorithm` -- (optional, default: `nothing`) the convex-hull algorithm; see
                 below for valid options
- `backend`   -- (optional, default: `nothing`) polyhedral computation backend
                 for higher-dimensional point sets
- `solver`    -- (optional, default: `nothing`) the linear-programming solver
                 used in the backend

### Output

The convex hull as a list of vectors with the coordinates of the points.

### Algorithm

A pre-processing step treats the cases with up to two points for one dimension
and up to four points for two dimensions.
For more points in one resp. two dimensions, we use more general algorithms.

For the one-dimensional case we return the minimum and maximum points, in that order.

The two-dimensional case is handled with a planar convex hull algorithm. The
following algorithms are available:

* `"monotone_chain"`        -- compute the convex hull of points in the plane
                               using Andrew's monotone chain method
* `"monotone_chain_sorted"` -- the same as `"monotone_chain"` but assuming that
                               the points are already sorted in counter-clockwise
                               fashion

See the reference docstring of each of those algorithms for details.

The higher dimensional case is treated using the concrete polyhedra library
`Polyhedra`, that gives access to libraries such as `CDDLib` and `ConvexHull.jl`.
These libraries can be chosen from the `backend` argument.

### Notes

For the in-place version use `convex_hull!` instead of `convex_hull`.

### Examples

Compute the convex hull of a random set of points:

```jldoctest ch_label
julia> points = [randn(2) for i in 1:30]; # 30 random points in 2D

julia> hull = convex_hull(points);

julia> typeof(hull)
Vector{Vector{Float64}} (alias for Array{Array{Float64, 1}, 1})
```
"""
function convex_hull(points::Vector{VN};
                     algorithm=nothing,
                     backend=nothing,
                     solver=nothing
                     ) where {N<:Real, VN<:AbstractVector{N}}
    return convex_hull!(copy(points), algorithm=algorithm, backend=backend, solver=solver)
end

function convex_hull!(points::Vector{VN};
                      algorithm=nothing,
                      backend=nothing,
                      solver=nothing) where {N<:Real, VN<:AbstractVector{N}}

    m = length(points)

    # =====================
    # Treat special cases
    # =====================

    # zero or one point
    if m == 0 || m == 1
        return points
    end

    n = length(first(points))

    if n == 1 # dimensional check
        if m == 2
            # two points case in 1d
            return _two_points_1d!(points)
        else
             # general case in 1d
            return _convex_hull_1d!(points)
        end
    elseif n == 2
        return _convex_hull_2d_preprocess!(points, m; algorithm=algorithm)
    else
        # general case in nd
        return _convex_hull_nd!(points, backend=backend, solver=solver)
    end
end

function _convex_hull_2d_preprocess!(points, m=length(points); algorithm=nothing)
    if m == 2
        # two points case in 2d
        return _two_points_2d!(points)
    elseif m == 3
        # three points case in 2d
        return _three_points_2d!(points)
    elseif m == 4
        # four points case in 2d
        return _four_points_2d!(points)
    else
        # general case in 2d
        return _convex_hull_2d!(points, algorithm=algorithm)
    end
end

function _two_points_1d!(points)
    p1, p2 = points[1], points[2]
    if  _isapprox(p1[1], p2[1]) # check for redundancy
        pop!(points)
    end
    return points
end

function _two_points_2d!(points)
    p1, p2 = points[1], points[2]
    if _isapprox(p1[1], p2[1]) && _isapprox(p1[2], p2[2]) # check for redundancy
        pop!(points)
    end
    return points
end

function _three_points_2d!(points::AbstractVector{<:AbstractVector{N}}) where {N<:Real}
    # Algorithm: the function takes three points and uses the formula
    #            from here: https://stackoverflow.com/questions/2122305/convex-hull-of-4-points/2122620#2122620
    #            to decide if the points are ordered in a counter-clockwise fashion or not, the result is saved
    #            in the 'turn' boolean, then returns the points in ccw fashion acting according to 'turn'. For the
    #            cases where the points are collinear we pass the points with the minimum and maximum first
    #            component to the function for two points in 2d(_two_points_2d), if those are equal, we do the same
    #            but with the second component.
    A, B, C = points[1], points[2], points[3]
    turn = right_turn(A, B, C)

    if isapproxzero(turn)
        # ABC are collinear
        if isapprox(A[1], B[1]) && isapprox(B[1], C[1]) && isapprox(C[1], A[1])
            # points are approximately equal in their first component
            if isapprox(A[2], B[2]) && isapprox(B[2], C[2]) && isapprox(C[2], A[2])
                # all points are approximately equal
                pop!(points)
                pop!(points)
                return points
            else
                points[1], points[2] = points[argmin([A[2], B[2], C[2]])], points[argmax([A[2], B[2], C[2]])]
            end
        else
            points[1], points[2] = points[argmin([A[1], B[1], C[1]])], points[argmax([A[1], B[1], C[1]])]
        end
        pop!(points)
        _two_points_2d!(points)
    elseif turn < zero(N)
        # ABC is CW
        points[1], points[2], points[3] = C, B, A
    end
    # else ABC is CCW => nothing to do
    return points
end

# given an index i in {1, 2, 3} return the element among (A, B, C) in the i-th position
@inline function _get_i(i, A, B, C)
    return i == 1 ? A : i == 2 ? B : C
end

function _collinear_case!(points, A, B, C, D)
    # A, B and C collinear, D is the extra point
    if isapprox(A[1], B[1]) && isapprox(B[1], C[1]) && isapprox(C[1], A[1])
        # points are approximately equal in their first component
        if isapprox(A[2], B[2]) && isapprox(B[2], C[2]) && isapprox(C[2], A[2])
            # the three points are approximately equal
            points[1], points[2] = A, D
            pop!(points)
            pop!(points)
            return _two_points_2d!(points)
        else
            # assign the points with max and min value in their second component to the
            # firsts points and the extra point to the third place, then pop the point that was in the middle
            min_y, max_y = arg_minmax(A[2], B[2], C[2])
            points[1] = _get_i(min_y, A, B, C)
            points[2] = _get_i(max_y, A, B, C)
            points[3] = D
            pop!(points)
        end
    else
        # assign the points with max and min value in their first component to the
        # firsts points and the extra point to the third place, then pop the point that was in the middle
        min_x, max_x = arg_minmax(A[1], B[1], C[1])
        points[1] = _get_i(min_x, A, B, C)
        points[2] = _get_i(max_x, A, B, C)
        points[3] = D
        pop!(points)
    end
    return _three_points_2d!(points)
end

function _four_points_2d!(points::AbstractVector{<:AbstractVector{N}}) where {N<:Real}
    A, B, C, D = points[1], points[2], points[3], points[4]
    tri_ABC = right_turn(A, B, C)
    tri_ABD = right_turn(A, B, D)
    tri_BCD = right_turn(B, C, D)
    tri_CAD = right_turn(C, A, D)
    key = 0
    if tri_ABC > zero(N)
        key = key + 1000
    end
    if tri_ABD > zero(N)
        key = key + 100
    end
    if tri_BCD > zero(N)
        key = key + 10
    end
    if tri_CAD > zero(N)
        key = key + 1
    end

    if isapproxzero(tri_ABC)
        return _collinear_case!(points, A, B, C, D)
    elseif isapproxzero(tri_ABD)
        return _collinear_case!(points, A, B, D, C)
    elseif isapproxzero(tri_BCD)
        return _collinear_case!(points, B, C, D, A)
    elseif isapproxzero(tri_CAD)
        return _collinear_case!(points, C, A, D, B)
    end

    # ABC  ABD  BCD  CAD  hull
    # ------------------------
    #  +    +    +    +   ABC
    #  +    +    +    -   ABCD
    #  +    +    -    +   ABDC
    #  +    +    -    -   ABD
    #  +    -    +    +   ADBC
    #  +    -    +    -   BCD
    #  +    -    -    +   CAD
    #  +    -    -    -   [should not happen]
    #  -    +    +    +   [should not happen]
    #  -    +    +    -   ACD
    #  -    +    -    +   DCB
    #  -    +    -    -   DACB
    #  -    -    +    +   ADB
    #  -    -    +    -   ACDB
    #  -    -    -    +   ADCB
    #  -    -    -    -   ACB
    if key == 1111
        points[1], points[2], points[3] = A, B, C # +    +    +    +   ABC
        pop!(points)
    elseif key == 1110
        points[1], points[2], points[3], points[4] = A, B, C, D # +    +    +    -   ABCD
    elseif key == 1101
        points[1], points[2], points[3], points[4] = A, B, D, C #  +    +    -    +   ABDC
    elseif key == 1100
        points[1], points[2], points[3] = A, B, D #  +    +    -    -   ABD
        pop!(points)
    elseif key == 1011
        points[1], points[2], points[3], points[4] = A, D, B, C #  +    -    +    +   ADBC
    elseif key == 1010
        points[1], points[2], points[3] = B, C, D #  +    -    +    -   BCD
        pop!(points)
    elseif key == 1001
        points[1], points[2], points[3] = C, A, D #  +    -    -    +    CAD
        pop!(points)
    elseif key == 0110
        points[1], points[2], points[3] = A, C, D #  -    +    +    -   ACD
        pop!(points)
    elseif key == 0101
        points[1], points[2], points[3] = D, C, B #  -    +    -    +   DCB
        pop!(points)
    elseif key == 0100
        points[1], points[2], points[3], points[4] = D, A, C, B #  -    +    -    -   DACB
    elseif key == 0011
        points[1], points[2], points[3] = A, D, B #  -    -    +    +   ADB
        pop!(points)
    elseif key == 0010
        points[1], points[2], points[3], points[4] = A, C, D, B #  -    -    +    -   ACDB
    elseif key == 0001
        points[1], points[2], points[3], points[4] = A, D, C, B #  -    -    -    +   ADCB
    elseif key == 0000
        points[1], points[2], points[3] = A, C, B #  -    -    -    -   ACB
        pop!(points)
    else
        @assert false "unexpected case in convex_hull"
    end
    return points
end

function _convex_hull_1d!(points::Vector{VN}) where {N<:Real, VN<:AbstractVector{N}}
    points[1:2] = [minimum(points), maximum(points)]
    return resize!(points, 2)
end

function _convex_hull_nd!(points::Vector{VN};
                          backend=nothing,
                          solver=nothing
                          ) where {N<:Real, VN<:AbstractVector{N}}
    V = VPolytope(points)
    Vch = remove_redundant_vertices(V, backend=backend, solver=solver)
    m = length(Vch.vertices)
    points[1:m] = Vch.vertices
    return resize!(points, m)
end

function _convex_hull_2d!(points::Vector{VN};
                          algorithm="monotone_chain"
                         ) where {N<:Real, VN<:AbstractVector{N}}
    if algorithm == nothing
        algorithm = default_convex_hull_algorithm(points)
    end
    if algorithm == "monotone_chain"
        return monotone_chain!(points)
    elseif algorithm == "monotone_chain_sorted"
        return monotone_chain!(points, sort=false)
    else
        error("the convex hull algorithm $algorithm is unknown")
    end
end

"""
    monotone_chain!(points::Vector{VN}; sort::Bool=true
                   ) where {N<:Real, VN<:AbstractVector{N}}

Compute the convex hull of points in the plane using Andrew's monotone chain
method.

### Input

- `points` -- list of 2D vectors; is sorted in-place inside this function
- `sort`   -- (optional, default: `true`) flag for sorting the vertices
              lexicographically; sortedness is required for correctness

### Output

List of vectors containing the 2D coordinates of the corner points of the
convex hull.

### Notes

For large sets of points, it is convenient to use static vectors to get
maximum performance. For information on how to convert usual vectors
into static vectors, see the type `SVector` provided by the
[StaticArrays](http://juliaarrays.github.io/StaticArrays.jl/stable/)
package.

### Algorithm

This function implements Andrew's monotone chain convex hull algorithm to
construct the convex hull of a set of ``n`` points in the plane in
``O(n \\log n)`` time.
For further details see
[Monotone chain](https://en.wikibooks.org/wiki/Algorithm_Implementation/Geometry/Convex_hull/Monotone_chain)
"""
function monotone_chain!(points::Vector{VN}; sort::Bool=true
                        ) where {N<:Real, VN<:AbstractVector{N}}

    @inline function build_hull!(semihull, iterator, points)
        @inbounds for i in iterator
            while length(semihull) >= 2 &&
                    (right_turn(semihull[end-1], semihull[end], points[i])
                         <= zero(N))
                pop!(semihull)
            end
            push!(semihull, points[i])
        end
    end

    if sort
        # sort the points lexicographically
        sort!(points, by=x->(x[1], x[2]))
    end

    # build lower hull
    lower = Vector{VN}()
    iterator = 1:length(points)
    build_hull!(lower, iterator, points)

    # build upper hull
    upper = Vector{VN}()
    iterator = length(points):-1:1
    build_hull!(upper, iterator, points)

    # remove the last point of each segment because they are repeated
    copyto!(points, @view(lower[1:end-1]))
    copyto!(points, length(lower), @view(upper[1:end-1]))
    m = length(lower) + length(upper) - 2
    if m == 2 && _isapprox(points[1], points[2])
        # upper and lower chain consist of a single, identical point
        m = 1
    end
    return resize!(points, m)
end

"""
    convex_hull(U::UnionSetArray{N, PT}; kwargs...) where {N, PT<:AbstractPolytope{N}}

Compute the convex hull of a union of a finite number of polytopes.

### Input

- `U` -- UnionSetArray of polytopes

### Output

A list of the vertices of the convex hull.
"""
function convex_hull(U::UnionSetArray{N, PT}; kwargs...) where {N, PT<:AbstractPolytope{N}}
    vlist = mapreduce(vertices_list, vcat, U.array)
    return convex_hull(vlist; kwargs...)
end

"""
    convex_hull(P::VPolygon, Q::VPolygon; [algorithm]::String="monotone_chain")

Return the convex hull of two polygons in vertex representation.

### Input

- `P`         -- polygon in vertex representation
- `Q`         -- another polygon in vertex representation
- `algorithm` -- (optional, default: "monotone_chain") the algorithm used to
                 compute the convex hull

### Output

A new polygon such that its vertices are the convex hull of the given two polygons.

### Algorithm

A convex hull algorithm is used to compute the convex hull of the vertices of the
given input polygons `P` and `Q`; see `?convex_hull` for details on the available
algorithms. The vertices of the output polygon are sorted in counter-clockwise
fashion.
"""
function convex_hull(P::VPolygon, Q::VPolygon;
                     algorithm::String="monotone_chain")
    vunion = [P.vertices; Q.vertices]
    convex_hull!(vunion; algorithm=algorithm)
    return VPolygon(vunion, apply_convex_hull=false)
end

"""
    convex_hull(P1::VPolytope, P2::VPolytope; [backend]=nothing)

Compute the convex hull of the set union of two polytopes in V-representation.

### Input

- `P1`         -- polytope
- `P2`         -- another polytope
- `backend`    -- (optional, default: `nothing`) the polyhedral
                  computations backend

### Output

The `VPolytope` obtained by the concrete convex hull of `P1` and `P2`.

### Notes

This function takes the union of the vertices of each polytope and then relies
on a concrete convex hull algorithm.
For low dimensions, a specialized implementation for polygons is used.
For higher dimensions, `convex_hull` relies on the polyhedral backend that can
be specified using the `backend` keyword argument.

For performance reasons, it is suggested to use the `CDDLib.Library()` backend.
"""
function convex_hull(P1::VPolytope, P2::VPolytope; backend=nothing)
    vunion = [P1.vertices; P2.vertices]
    convex_hull!(vunion; backend=backend)
    return VPolytope(vunion)
end

"""
    convex_hull(P1::HPoly, P2::HPoly;
               [backend]=default_polyhedra_backend(P1))

Compute the convex hull of the set union of two polyhedra in H-representation.

### Input

- `P1`         -- polyhedron
- `P2`         -- another polyhedron
- `backend`    -- (optional, default: `default_polyhedra_backend(P1)`)
                  the polyhedral computations backend

### Output

The `HPolyhedron` (resp. `HPolytope`) obtained by the concrete convex hull of
`P1` and `P2`.

### Notes

For performance reasons, it is suggested to use the `CDDLib.Library()` backend
for the `convex_hull`.

For further information on the supported backends see
[Polyhedra's documentation](https://juliapolyhedra.github.io/).
"""
function convex_hull(P1::HPoly, P2::HPoly;
                     backend=default_polyhedra_backend(P1))
    require(:Polyhedra; fun_name="convex_hull")
    Pch = Polyhedra.convexhull(polyhedron(P1; backend=backend),
                               polyhedron(P2; backend=backend))
    removehredundancy!(Pch)
    return convert(basetype(P1), Pch)
end


text:
function init_visualization()
    traj = []
    obj_traj = []
    return traj, obj_traj
end

function visualize_environment!(objects)
    fig = Figure()
    ax = Axis(fig[1,1], aspect=DataAspect(), limits=(-10.,10.,-10.,10.))
    obs_dict = Dict()
    for ob in objects
        ox = Observable(SVector{2,Float64}(ob.pose...))
        θ = Observable(0.0)
        scatter!(ax, ox; rotation=θ, marker=:rect, markersize=25, color=ob.color)
        obs_dict[ob.id] = (ox,θ)
    end
    return obs_dict, ax, fig
end

function visualize_trajectory!(bobby, trajectory, obj_trajectory, obs_dict, ax, fig;
                                name="media/test.gif")
    T = length(trajectory)
    θ = Observable(0.0)
    x = Observable(SVector{2,Float64}(0.,0.))
    left_img_bg_path = joinpath(@__DIR__, "..", "models/left_swing_bg.png")
    right_img_bg_path = joinpath(@__DIR__, "..", "models/right_swing_bg.png")
    img = Observable(load(left_img_bg_path))
    scatter!(ax, x; rotation=θ, marker=img, markersize=40)
    # hidedecorations!(ax)

    function gaits!(t)
        stride=10
        x[] = SVector{2,Float64}(trajectory[t][1:2]...)
        θ[] = trajectory[t][3] - π/2.

        if !isnothing(obj_trajectory[t])
            ind = obj_trajectory[t][4]
            ox, oθ = obs_dict[ind]
            ox[] = SVector{2,Float64}(obj_trajectory[t][1:2]...)
            oθ[] = obj_trajectory[t][3] - π/2.
        end

        img[] = t % (stride*2) <= stride ? load(left_img_bg_path) : load(right_img_bg_path)
        sleep(0.01)
    end

    record(fig, name, 1:T, framerate=60) do t
        gaits!(t)
    end
    empty!(ax)

end


text: include("FPGrowth.jl")

max_stars_repo_name: DopedCow/FPGrowth

max_issues_repo_name: DopedCow/FPGrowth

max_forks_repo_name: DopedCow/FPGrowth


text:
"""
    $(SIGNATURES)
Create a {Float32} point cloud from a depth image clampin the z-axis to [clampz]. Note: rotated to Forward Starboard Down
"""
function cloudFromDepthImageClampZ(depths::Array{UInt16,2},
                                   cm::CameraModel,
                                   trans::AffineMap;
                                   depthscale = 0.001f0,
                                   skip::Int = 1,
                                   maxrange::Float32 = 5f0,
                                   clampz = [0f0,1f0],
								   colmap::ColorGradient = cgrad(:bgy,:colorcet))::PointCloud
	  							   # colmap::ColorGradient = ColorGradient([RGBA(0,1,0,1)], [1.0]))::PointCloud # I prever default to be pretty

    #
	cx = Float32(cm.cc[1])
	cy = Float32(cm.cc[2])
	fx = Float32(cm.fc[1])
	fy = Float32(cm.fc[2])
    (row,col) = size(depths)
    cloud = Point3f0[]
    cloudCol = RGB{Float32}[]

    for u = 1:skip:row, v = 1:skip:col
        z = depths[u,v]*depthscale
        if  0 < z < maxrange
            x = (v-cx)/fx * z
            y = (u-cy)/fy * z
            p = trans(Point3f0(z,x,y))
			if clampz[1] <= p[3] <= clampz[2]
            	push!(cloud, p) #NOTE rotated to Forward Starboard Down, TODO: maybe leave in camera frame?
				cmlength = length(colmap.colors)
				push!(cloudCol, colmap.colors[round(Int,((p[3]-clampz[1])/(clampz[2]-clampz[1]))*Float32(cmlength-1))+1])
			end
        end
    end
    pointcloud = PointCloud(cloud, cloudCol)

    return pointcloud
end

"""
    $(SIGNATURES)
Create a {Float32} point cloud from a depth image, repeating every [colmapz = 1f0]. Note: rotated to Forward Starboard Down
""" #TODO fix color map
function cloudFromDepthImage(depths::Array{UInt16,2},
							 cm::CameraModel,
							 trans::AffineMap=Translation(0,0,0);
							 depthscale = 0.001f0,
							 skip::Int = 1,
							 maxrange::Float32 = 5f0,
							 colmapz::Float32 = 1f0,
							 colmap::ColorGradient = cgrad(:bgy,:colorcet))::PointCloud
							 # colmap::ColorGradient = ColorGradient([RGBA(0,1,0,1)], [1.0]))::PointCloud # I prever default to be pretty

	cx = Float32(cm.cc[1])
	cy = Float32(cm.cc[2])
	fx = Float32(cm.fc[1])
	fy = Float32(cm.fc[2])
    (row,col) = size(depths)

    cloud = Point3f0[]
    cloudCol = RGB{Float32}[]

	cmlength = length(colmap.colors)

    for u = 1:skip:row, v = 1:skip:col
        z = depths[u,v]*depthscale
        if  0 < z < maxrange
            x = (v-cx)/fx * z
			y = (u-cy)/fy * z
			p = trans(Point3f0(z,x,y))
            push!(cloud, p) #NOTE rotated to Forward Starboard Down, TODO: maybe leave in camera frame?
			push!(cloudCol, colmap.colors[round(Int,(mod(p[3],colmapz))*Float32(cmlength-1))+1])
        end
    end
    pointcloud = PointCloud(cloud, cloudCol)

    return pointcloud
end

max_stars_repo_name: Affie/Amphitheatre

max_issues_repo_name: Affie/Amphitheatre


text:
using Plots
using Printf
using DifferentialEquations
using DiffEqBiological

seirmodel = @reaction_network SEIR begin
    β, S + I --> E + I
    δ, E --> I
    γ, I --> R
    ϵ, R --> S
end β δ γ ϵ

function solveSEIR(u0, β, γ, δ, ϵ, time_span)

    dprob = DiscreteProblem(seirmodel, u0, time_span, [β, δ, γ, ϵ] )
    jprob = JumpProblem(dprob, Direct(), seirmodel)

    sol = solve(jprob, SSAStepper())

    return sol

end

R0=2.0
N=10000
I=10
γ=0.1
β=γ*R0/N
δ=0.25
ϵ=1/150

sol = solveSEIR([N-I,0,I,0], β, γ, δ, ϵ, (0.0,1080.0))

xlabel = Float64[]

ylabel, yval = Float64[], Float64[]

for i in 1:1080
    rem(i,30)==0 && push!(xlabel, i)
end

# y axis is log scale
for i in 1:10, j in 1:4
    push!(yval, i*10^j)
    i==1 && push!(ylabel, 10^j)
end

title=@sprintf "SEIR (Expire) model/Gillespie R0=%.2f γ=%.3f δ=%.3f ϵ=%.5f" R0 γ δ ϵ
plot_graph = plot(sol, vars=(0,1), title=title, label="Suscepible",
                  color=:blue, seriestype=:line, legend=:bottomright, size=(1200,900),
                  xticks=xlabel, yticks=(yval, ylabel), yaxis=(:log10, (1,Inf)))
plot!(plot_graph, sol, vars=(0,2), label="Infectious", color=:red)
plot!(plot_graph, sol, vars=(0,3), label="Exposed", color=:darkmagenta)
plot!(plot_graph, sol, vars=(0,4), label="Removed", color=:green)

png(plot_graph, "seir_gillespie_exp.png")

max_stars_repo_name: seisyuu-hantatsushi/julia

max_issues_repo_name: seisyuu-hantatsushi/julia

max_forks_repo_name: seisyuu-hantatsushi/julia


text:
module CoverageBase
using Coverage
export testnames, runtests

const need_inlining = ["reflection", "meta"]

function julia_top()
    dir = joinpath(JULIA_HOME, "..", "share", "julia")
    if isdir(joinpath(dir,"base")) && isdir(joinpath(dir,"test"))
        return dir
    end
    dir = JULIA_HOME
    while !isdir(joinpath(dir,"base"))
        dir, _ = splitdir(dir)
        if dir == "/"
            error("Error parsing top dir; JULIA_HOME = $JULIA_HOME")
        end
    end
    dir
end
const topdir = julia_top()
const basedir = joinpath(topdir, "base")
const testdir = joinpath(topdir, "test")
include(joinpath(testdir, "choosetests.jl"))

function testnames()
    Base.JLOptions().can_inline == 1 && return need_inlining

    names, _ = choosetests()
    filter!(x -> !in(x, need_inlining), names)

    # Manually add in `pkg`, which is disabled so that `make testall` passes on machines without internet access
    push!(names, "pkg")
    names
end

function runtests(names)
    cd(testdir) do
        for tst in names
            println(tst)
            include(tst*".jl")
        end
    end
end

end # module


text:
export lossLayer

type lossLayer <:Layer
  bottom   ::Array{Layer}
  top      ::Array{Layer}
  topRange ::Vector

  rₒ      ::Array{Float64}
  rᵢ      ::Array{Float64} #prediction

  Cost     ::Function # Cost function
  ∂Cost    ::Function # derivative of activation function
  loss     ::Float64

  ∂rᵢ_L  ::Array{Float64}

  lock   ::Bool
end

lossLayer(nClass::Int)=
  lossLayer(Layer[],
            Layer[],
            Vector[],

			rand(nClass,batchsize),
			rand(nClass,batchsize),

            leastsq, # crossEntropy
            ∂leastsq,#∂crossEntropy
            0.1,

			rand(nClass,batchsize),

			false)


text:
using JuMP
using LinearAlgebra

# * Free and open-source solvers
using GLPK
# using Cbc
# * Commercial solvers
# using Gurobi

function solve_lp_assignment_problem(costMatrix::AbstractMatrix; optimizer=GLPK.Optimizer, verbose=false,
		binary_var=false, category=Nothing)
    rowNum, colNum = size(costMatrix)
	@assert rowNum <= colNum "Inventory size is smaller than the design set size. Inventory size must be larger."

	# Preparing an optimization model
	m = Model(optimizer)
	# set_optimizer_attribute(m, "OutputFlag", verbose)
	if !verbose
		set_silent(m)
	end

	# https://nbviewer.jupyter.org/github/jump-dev/JuMPTutorials.jl/blob/master/notebook/modelling/network_flows.ipynb
	# Declaring variables
	if !binary_var
		m[:x] = @variable(m, 0 <= x[1:rowNum, 1:colNum] <= 1)
	else
		m[:x] = @variable(m, x[1:rowNum, 1:colNum], Bin)
	end

	# Setting the objective
	@objective(m, Min, dot(costMatrix, x))

	# Adding constraints
	# * each design unit must be matched with a inventory unit
	m[:design] = @constraint(m, [i = 1:rowNum], sum(x[i, :]) == 1)

	# * each inventory unit can be used at most once
	if rowNum == colNum
		# * one-to-one match
		m[:inventory_eq] = @constraint(m, [j = 1:colNum], sum(x[:, j]) == 1)
	else
		m[:inventory_ineq] = @constraint(m, [j = 1:colNum], sum(x[:, j]) <= 1)
	end

	# * additional categorical constraint
	if category isa Vector
		# ! check one-based indexing
		@assert all(map(v->all(v .>= 1), category))
		m[:inventory_category] = @constraint(m, [k = 1:length(category)], sum(sum(x[:, j]) for j in category[k]) <= 1)
	elseif category isa Matrix
		# ! pycall convert non-ragged vector of vectors into a matrix...
		@assert all(category .>= 1)
		m[:inventory_category] = @constraint(m, [k = 1:size(category, 1)], sum(sum(x[:, j]) for j in category[k,:]) <= 1)
	end

	#  the optimization problem
	optimize!(m)
	if termination_status(m) == MOI.OPTIMAL
	    x_opt = value.(x)
	    cost = objective_value(m)
	elseif termination_status(m) == MOI.TIME_LIMIT && has_values(m)
	    x_opt = value.(x)
	    cost = objective_value(m)
		println("The model terminates due to time limits.")
	    error("The model terminates due to time limits.")
	else
	    error("The model was not solved correctly.")
	end

    assignment = zeros(Int, rowNum)
    for i = 1:rowNum
		for j =1:colNum
        	if x_opt[i, j] ≈ 1.0
        	        assignment[i] = j
        	end
		end
    end
    return assignment, cost, m
end

function solve_upcycling_assignment(design_lengths::AbstractVector, inventory_lengths::AbstractVector,
		match_costs::AbstractMatrix, cutoff_costs::AbstractVector;
		optimizer=Gurobi.Optimizer,
		binary_var=true, allow_cutoff=false, minimize_cutoff=true, verbose=false)
	@assert length(design_lengths) <= length(inventory_lengths) "We must have more inventory elements ($(length(inventory_lengths))) than design elements ($(length(design_lengths)))."
    rowNum, colNum = length(design_lengths), length(inventory_lengths)
	@assert (rowNum, colNum) == size(match_costs) && colNum == length(cutoff_costs)

	# Preparing an optimization model
	m = Model(optimizer)
	if !verbose
		set_silent(m)
	end

	# Declaring variables
	if !binary_var
		@variable(m, 0 <= x[1:rowNum, 1:colNum] <= 1)
		if allow_cutoff
			@variable(m, 0 <= y[1:colNum] <= 1)
		end
	else
		@variable(m, x[1:rowNum, 1:colNum], Bin)
		if allow_cutoff
			@variable(m, y[1:colNum], Bin)
		end
	end

	# Setting the objective, minimize wasted materials
	# @objective(m, Min, sum(x * inventory_lengths - design_lengths) + dot(cutoff_costs, y)) # + dot(match_costs, x)
	if allow_cutoff
		if minimize_cutoff
			@objective(m, Min, sum(y .* inventory_lengths - x' * design_lengths) + dot(cutoff_costs, y))
        else
			@objective(m, Min, dot(match_costs, x) + dot(cutoff_costs, y))
		end
	else
		# minimize cut-off waste
		@objective(m, Min, dot(match_costs, x))
		# equivalent to
		# @objective(m, Min, sum(x * inventory_lengths - design_lengths))
	end

	# * Adding constraints
	# For each member i , only one stock element j is reused
	@constraint(m, [i = 1:rowNum], sum(x[i, :]) == 1)
	if allow_cutoff
		# The use of stock element j \in R for one or more members is constrained by the available length
		@constraint(m, [j = 1:colNum], dot(x[:, j], design_lengths) <= y[j]*inventory_lengths[j])
	else
		# each stock element can only be used once
		@constraint(m, [j = 1:colNum], sum(x[:, j]) <= 1)
		# oversize constraint
		@constraint(m, x*inventory_lengths .>= design_lengths)
	end

	#  the optimization problem
	optimize!(m)
	if verbose
		solution_summary(m, verbose=true)
	end
	if termination_status(m) == MOI.OPTIMAL
		#
	elseif termination_status(m) == MOI.TIME_LIMIT && has_values(m)
		println("The model terminates due to time limits.")
	else
	    error("The model was not solved correctly: $(termination_status(m))")
	end

	x_opt = value.(x)
	y_opt = allow_cutoff ? value.(y) : zeros(colNum)
	# cost = objective_value(m)
	# if allow_cutoff
	# 	cost -= dot(cutoff_costs, y_opt)
	# end
	# @show x_opt

    assignment = zeros(Int, rowNum)
    for i = 1:rowNum
		for j = 1:colNum
        	if abs(x_opt[i, j]-1.0) < 1e-4
        	# if x_opt[i, j] ≈ 1.0
        	    if rowNum ≤ colNum
        	        assignment[i] = j
					if !allow_cutoff
						y_opt[j] = 1
					end
        	    else
        	        assignment[j] = i
					if !allow_cutoff
						y_opt[i] = 1
					end
        	    end
				# break
        	end
		end
    end
	cutoff_waste = sum(y_opt .* inventory_lengths - x_opt' * design_lengths)
    return assignment, (x_opt, y_opt), cutoff_waste
end


text:
export PLASLearner

mutable struct PLASLearner{
    BA1<:NeuralNetworkApproximator,
    BA2<:NeuralNetworkApproximator,
    BC1<:NeuralNetworkApproximator,
    BC2<:NeuralNetworkApproximator,
    V<:NeuralNetworkApproximator,
    R<:AbstractRNG,
} <: AbstractLearner
    policy::BA1
    target_policy::BA2
    qnetwork1::BC1
    qnetwork2::BC2
    target_qnetwork1::BC1
    target_qnetwork2::BC2
    vae::V
    γ::Float32
    τ::Float32
    λ::Float32
    batch_size::Int
    pretrain_step::Int
    update_freq::Int
    update_step::Int
    rng::R
    # Logging
    actor_loss::Float32
    critic_1_loss::Float32
    critic_2_loss::Float32
end

"""
    PLASLearner(;kwargs...)

See [Latent Action Space for Offline Reinforcement Learning](https://arxiv.org/abs/2011.07213)

# Keyword arguments
- `policy`, used to get latent action.
- `target_policy`, similar to `policy`, but used to estimate the target.
- `qnetwork1`, used to get Q-values.
- `qnetwork2`, used to get Q-values.
- `target_qnetwork1`, used to estimate the target Q-values.
- `target_qnetwork2`, used to estimate the target Q-values.
- `vae`, used for mapping hidden actions to actions. This
can be implemented using a `VAE` in a `NeuralNetworkApproximator`.
- `γ::Float32 = 0.99f0`, reward discount rate.
- `τ::Float32 = 0.005f0`, the speed at which the target network is updated.
- `λ::Float32 = 0.75f0`, used for Clipped Double Q-learning.
- `batch_size::Int = 32`
- `pretrain_step::Int = 1000`, the number of pre-training rounds.
- `update_freq::Int = 50`, the frequency of updating the `approximator`.
- `update_step::Int = 0`
- `rng = Random.GLOBAL_RNG`
"""
function PLASLearner(;
    policy,
    target_policy,
    qnetwork1,
    qnetwork2,
    target_qnetwork1,
    target_qnetwork2,
    vae,
    γ = 0.99f0,
    τ = 0.005f0,
    λ = 0.75f0,
    batch_size = 32,
    pretrain_step = 10000,
    update_freq = 50,
    update_step = 0,
    rng = Random.GLOBAL_RNG,
)
    copyto!(policy, target_policy)  # force sync
    copyto!(qnetwork1, target_qnetwork1)  # force sync
    copyto!(qnetwork2, target_qnetwork2)  # force sync
    PLASLearner(
        policy,
        target_policy,
        qnetwork1,
        qnetwork2,
        target_qnetwork1,
        target_qnetwork2,
        vae,
        γ,
        τ,
        λ,
        batch_size,
        pretrain_step,
        update_freq,
        update_step,
        rng,
        0.0f0,
        0.0f0,
        0.0f0,
    )
end

function (l::PLASLearner)(env)
    s = send_to_device(device(l.policy), state(env))
    s = Flux.unsqueeze(s, ndims(s) + 1)
    latent_action = tanh.(l.policy(s))
    action = dropdims(decode(l.vae.model, s, latent_action), dims=2)
end

function RLBase.update!(l::PLASLearner, batch::NamedTuple{SARTS})
    if l.update_step == 0
        update_vae!(l, batch)
    else
        update_learner!(l, batch)
    end
end

function update_vae!(l::PLASLearner, batch::NamedTuple{SARTS})
    s, a, r, t, s′ = send_to_device(device(l.vae), batch)
    a = reshape(a, :, l.batch_size)
    vae_grad = gradient(Flux.params(l.vae)) do
        recon_loss, kl_loss = vae_loss(l.vae.model, s, a)
        0.5f0 * kl_loss + recon_loss
    end
    update!(l.vae, vae_grad)
end

function update_learner!(l::PLASLearner, batch::NamedTuple{SARTS})
    s, a, r, t, s′ = send_to_device(device(l.qnetwork1), batch)

    γ, τ, λ = l.γ, l.τ, l.λ

    latent_action′ = tanh.(l.target_policy(s′))
    action′ = decode(l.vae.model, s′, latent_action′)
    q′_input = vcat(s′, action′)
    q′ = λ .* min.(l.target_qnetwork1(q′_input), l.target_qnetwork2(q′_input)) + (1 - λ) .* max.(l.target_qnetwork1(q′_input), l.target_qnetwork2(q′_input))

    y = r .+ γ .* (1 .- t) .* vec(q′)

    # Train Q Networks
    a = reshape(a, :, l.batch_size)
    q_input = vcat(s, a)

    q_grad_1 = gradient(Flux.params(l.qnetwork1)) do
        q1 = l.qnetwork1(q_input) |> vec
        loss = mse(q1, y)
        ignore() do
            l.critic_1_loss = loss
        end
        loss
    end
    update!(l.qnetwork1, q_grad_1)

    q_grad_2 = gradient(Flux.params(l.qnetwork2)) do
        q2 = l.qnetwork2(q_input) |> vec
        loss = mse(q2, y)
        ignore() do
            l.critic_2_loss = loss
        end
        loss
    end
    update!(l.qnetwork2, q_grad_2)

    # Train Policy
    p_grad = gradient(Flux.params(l.policy)) do
        latent_action = tanh.(l.policy(s))
        action = decode(l.vae.model, s, latent_action)
        actor_loss = -mean(l.qnetwork1(vcat(s, action)))
        ignore() do
            l.actor_loss = actor_loss
        end
        actor_loss
    end
    update!(l.policy, p_grad)

    # polyak averaging
    for (dest, src) in zip(
        Flux.params([l.target_policy, l.target_qnetwork1, l.target_qnetwork2]),
        Flux.params([l.policy, l.qnetwork1, l.qnetwork2]),
    )
        dest .= (1 - τ) .* dest .+ τ .* src
    end
end


text:
nums = []
for line in eachline()
    num1 = parse(Int32, line)
    for num2 in nums
        if num2 + num1 == 2020
            println(num1 * num2)
            return
        end
    end
    push!(nums, num1)
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Functions related to the orbit state vector.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

export orbsv, sv_to_kepler

################################################################################
#                                     API
################################################################################

# Direct getters.
@inline get_epoch(sv::OrbitStateVector) = sv.t
@inline get_r(sv::OrbitStateVector) = sv.r
@inline get_v(sv::OrbitStateVector) = sv.v
@inline get_rv(sv::OrbitStateVector) = sv.r, sv.r

# Getters that require conversions.
@inline get_a(sv::OrbitStateVector) = get_a(sv_to_kepler(sv))
@inline get_e(sv::OrbitStateVector) = get_e(sv_to_kepler(sv))
@inline get_i(sv::OrbitStateVector) = get_i(sv_to_kepler(sv))
@inline get_Ω(sv::OrbitStateVector) = get_Ω(sv_to_kepler(sv))
@inline get_ω(sv::OrbitStateVector) = get_ω(sv_to_kepler(sv))
@inline get_f(sv::OrbitStateVector) = get_f(sv_to_kepler(sv))
@inline get_M(sv::OrbitStateVector) = get_M(sv_to_kepler(sv))

# Conversions.
"""
    sv_to_kepler(sv::OrbitStateVector)

Convert the state vector `sv` to Keplerian elements represented by an instance
of the structure `KeplerianElements`.

"""
@inline sv_to_kepler(sv::OrbitStateVector) = rv_to_kepler(sv.r, sv.v, sv.t)

################################################################################
#                                Initialization
################################################################################

"""
    orbsv(t::Tepoch, r::AbstractVector{T1}, v::AbstractVector{T2}, a::AbstractVector{T3} = [0,0,0]) where {Tepoch<:Number, T1<:Number, T2<:Number, T3<:Number}
    orbsv(t::Number, vec::AbstractVector)

Create a new satellite state vector (see [`OrbitStateVector`](@ref)) using the
position `r`, velocity `v`, and acceleration `a`. It is also possible to pass a
vector `vec` with the information concatenated.

!!! info
    The vectors `r`, `v`, and `a` must have at least 3 elements. In the case
    more elements are available, they will be neglected. On the other hand, the
    vector `v` must have 6 or 9 dimensions, indicating `[r; v]`, or `[r; v; a]`.
"""
function orbsv(
    t::Tepoch,
    r::AbstractVector{T1},
    v::AbstractVector{T2},
    a::AbstractVector{T3} = [0,0,0]
) where {Tepoch<:Number, T1<:Number, T2<:Number, T3<:Number}

    len_r = length(r)
    len_v = length(v)
    len_a = length(a)

    len_r  < 3 && error("The vector `r` must have at least 3 elements.")
    len_v  < 3 && error("The vector `v` must have at least 3 elements.")
    len_a  < 3 && error("The vector `a` must have at least 3 elements.")
    len_r != 3 && @warn("Only the first 3 elements of the vector `r` will be used!")
    len_v != 3 && @warn("Only the first 3 elements of the vector `v` will be used!")
    len_a != 3 && @warn("Only the first 3 elements of the vector `a` will be used!")

    T = promote_type(T1, T2, T3)
    return OrbitStateVector{Tepoch, T}(t = t, r = r[1:3], v = v[1:3], a = a[1:3])
end

function orbsv(t::Number, vec::AbstractVector)
    len_vec = length(vec)

    len_vec ∉ [6,9] && error("The length of input vector must be 6, or 9.")

    if len_vec == 6
        return OrbitStateVector(t = t, r = vec[1:3], v = vec[4:6])
    else
        return OrbitStateVector(t = t, r = vec[1:3], v = vec[4:6], a = vec[7:9])
    end
end

################################################################################
#                                  Overloads
################################################################################

function getindex(sv::OrbitStateVector{Tepoch, T}, ::Colon) where {Tepoch, T}
    return SVector{9, T}(sv.r..., sv.v..., sv.a...)
end

################################################################################
#                                      IO
################################################################################

function show(io::IO, sv::OrbitStateVector{Tepoch, T}) where {Tepoch, T}
    compact = get(io, :compact, true)
    epoch_str = sprint(print, sv.t, context = :compact => compact)
    jd_str = sprint(print, jd_to_date(DateTime, sv.t))
    print(io, "OrbitStateVector{", Tepoch, ", ", T, "}: Epoch = $epoch_str ($jd_str)")
end

function show(
    io::IO,
    mime::MIME"text/plain",
    sv::OrbitStateVector{Tepoch, T}
) where {Tepoch, T}
    # Check if the `io` supports colors.
    color = get(io, :color, false)

    # Compact printing.
    compact = get(io, :compact, true)

    b = (color) ? _b : ""
    d = (color) ? _d : ""

    t_str  = sprint(print, sv.t, context = :compact => compact)
    JD_str = sprint(print, jd_to_date(DateTime, sv.t), context = :compact => compact)
    r_str  = sprint(print, sv.r./1000, context = :compact => compact)
    v_str  = sprint(print, sv.v./1000, context = :compact => compact)

    # Add units.
    max_length = max(length(r_str), length(v_str))
    r_str *= " "^(max_length - length(r_str) + 1) * " km"
    v_str *= " "^(max_length - length(v_str) + 1) * " km/s"

    println(io, "OrbitStateVector{", Tepoch, ", ", T, "}: ")
    println(io, "$b  epoch :$d ", t_str, " (", JD_str, ")")
    println(io, "$b      r :$d ", r_str)
    print(io,   "$b      v :$d ", v_str)

    return nothing
end


text:
include( "../src/bootstrapPosteriorCorEstimation.jl" )

using DataFrames
using Distributions
using Random
using Statistics
using Test
using Optim
using MultivariateStats

function distortVector( vec::Vector{Float64}, pct::Float64, σ::Float64 )
  res          = deepcopy( vec );
  len          = size( vec, 1 );
  indices      = rand( 1:1:len, round( Int64, pct*len ) );
  gaussUniDist = Normal( 0.0, σ )
  for i in indices
    res[i] = vec[i] + rand( gaussUniDist );
  end
  return res;
end

println( "running bootstrap correlation algorithm unit tests." );
predMatTraining, predMatEval, yTraining, yEval, errMatTraining = makeupPredictions();
posterior = bootstrapPosteriorCorEstimation( predMatTraining, yTraining, 64, 10000 );
@test all( posterior .>= 0.0 );
@test all( posterior .<= 1.0 );
@test all( sum( posterior ) ≈ 1.0 );

nSamples    = 10000;
X           = rand( Normal( 0.0, 1.5 ), nSamples );
t_1         = Float64[ ( 0.3 * x^3 + 4.0 ) for x in X ];
t_2         = Float64[ ( 0.2 * x^2 + 0.1 * x ) for x in X ];
t_3         = Float64[ ( 0.5 * x^3 + 3.4 * x^2 + 1.5 * x + 0.95 ) for x in X ];
TMat        = hcat( t_1, t_2, t_3 );
predictions = Vector{ Matrix{Float64} }();
mat1 = hcat( distortVector( t_1, 0.05, 1.85 ),
             distortVector( t_2, 0.30, 0.75 ),
             distortVector( t_3, 0.25, 0.5 ) );
push!( predictions, mat1 );

mat2 = hcat( distortVector( t_1, 0.35, 1.25 ),
             distortVector( t_2, 0.10, 0.95 ),
             distortVector( t_3, 0.45, 0.35 ) );
push!( predictions, mat2 );

mat3 = hcat( distortVector( t_1, 1.35, 2.25 ),
             distortVector( t_2, 0.03, 0.07 ),
             distortVector( t_3, 0.25, 5.0 ) );
push!( predictions, mat3  );
nrRuns         = 5000;
samplingFactor = 0.4;
posterior      = bootstrapPosteriorCorEstimation( predictions, TMat, samplingFactor, nrRuns );
@test all( posterior .>= 0.0 );
@test all( posterior .<= 1.0 );
@test all( sum( posterior ) ≈ 1.0 );


text:
const _apex_data_file   = string(@__DIR__, "/apexsh.dat")
const _apex_so_filename =
            string(
                    @__DIR__, "/apex",
                    isapple ? "_macos" : islinux ? "_linux" : "_win",
                    ".so",
                )

_check_apex() = _check_so(_apex_so_filename)
_check_apex()

# export _loadapexsh
_loadapexsh(year_32::Float32) =
        ccall((:loadapxsh_, _apex_so_filename),
                Nothing,
                (Ptr{UInt8}, Ref{Float32},),
                _apex_data_file, year_32)

# export _apxg2q
_apxg2q(glat::Float32, glon::Float32,
        height::Float32, vecflagin::Int32,
        qlatout::Ref{Float32}, qlonout::Ref{Float32},
        f1::Ref{Float32}, f2::Ref{Float32}, f::Ref{Float32}) =
            ccall((:apxg2q_, _apex_so_filename),
                Nothing,
                (Ref{Float32}, Ref{Float32},
                Ref{Float32}, Ref{Int32},
                Ref{Float32}, Ref{Float32},
                Ref{Float32}, Ref{Float32}, Ref{Float32},),
                    glat, glon,
                    height, vecflagin,
                    qlatout, qlonout,
                    f1, f2, f,)

# export _apxq2g
_apxq2g(qlat::Float32, qlon::Float32,
        height::Float32, precision::Float32,
        glatout::Ref{Float32}, glonout::Ref{Float32},
        error::Ref{Float32}) =
            ccall((:apxq2g_, _apex_so_filename),
                    Nothing,
                    (Ref{Float32}, Ref{Float32},
                    Ref{Float32}, Ref{Float32},
                    Ref{Float32}, Ref{Float32}, Ref{Float32},),
                        qlat, qlon,
                        height, precision,
                        glatout, glonout, error,)

function test_apex()
    dt = Dates.DateTime(2010, 6, 6, 1, 2, 3)
    year = 2010
    yearf = Dates.year(dt) + Dates.dayofyear(dt) / Dates.daysinyear(dt)
    year32f = Float32(yearf)
    glat = Float32(25.0)
    glon = Float32(11.0)
    height = Float32(350.0)
    vecflagin = Int32(0)
    qlatout = Ref{Float32}(0.0)
    qlonout = Ref{Float32}(0.0)
    f1 = Ref{Float32}(0.0)
    f2 = Ref{Float32}(0.0)
    f = Ref{Float32}(0.0)

    _loadapexsh(year32f)
    _apxg2q(glat, glon, height, vecflagin,
            qlatout, qlonout, f1, f2, f)

    qlat = qlatout[]
    qlon = qlonout[]
    precision = 1.0f-10
    glatout = Ref{Float32}(0.0)
    glonout = Ref{Float32}(0.0)
    gerror  = Ref{Float32}(0.0)
    _apxq2g(qlat, qlon, height, precision,
            glatout, glonout,  gerror)

    # println("# ", qlatout[], ", ", qlonout[], ", ", f1[], ", ", f2[], ", ", f[])
    # println("# ", glatout[], ",\t", glonout[], ",\t", gerror[])
    # println("# ", glat, ",\t", glon, )
    # 16.185024, 85.22224, -9999.0, -9999.0, -9999.0
    # 25.0, 11.000001, 0.0

    # @test isapprox(qlatout[], 16.185024)
    # @test isapprox(qlonout[], 85.22224)
    # @test isapprox(glatout[], 25.0)
    # @test isapprox(glonout[], 11.000001)
    qlatout[], qlonout[], glatout[], glonout[]

end


text:
function predict(net::Vector{Layer},input)
  net[1].rₒ=input
  for idx=1:length(net)-1
    forward!(net[idx])
  end
  return net[end-1].rₒ
end


text:
mutable struct MapSimulateState{T,U}
    score::Float64
    noise::Float64
    subtraces::Vector{U}
    retval::Vector{T}
    num_nonempty::Int
end

function process!(gen_fn::Map{T,U}, args::Tuple,
                  key::Int, state::MapSimulateState{T,U}) where {T,U}
    local subtrace::U
    local retval::T
    kernel_args = get_args_for_key(args, key)
    subtrace = simulate(gen_fn.kernel, kernel_args)
    state.noise += project(subtrace, EmptySelection())
    state.num_nonempty += (isempty(get_choices(subtrace)) ? 0 : 1)
    state.score += get_score(subtrace)
    state.subtraces[key] = subtrace
    retval = get_retval(subtrace)
    state.retval[key] = retval
end

function simulate(gen_fn::Map{T,U}, args::Tuple) where {T,U}
    len = length(args[1])
    state = MapSimulateState{T,U}(0., 0., Vector{U}(undef,len), Vector{T}(undef,len), 0)
    for key=1:len
        process!(gen_fn, args, key, state)
    end
    VectorTrace{MapType,T,U}(gen_fn,
        PersistentVector{U}(state.subtraces), PersistentVector{T}(state.retval),
        args, state.score, state.noise, len, state.num_nonempty)
end


text:
using ProjC
using Test

@testset "ProjC.jl" begin
    # Write your own tests here.
end


text:
using StructuralEstimation
srand(122)

# 2 Model
bdisc = 0.95 # discount factor
truepar = [-1.9;1.;2.] # fixed costs, ..., entry costs
# State 1: Exogenous state
nX1 = 2
S1 = EntryState(;dense=false)

nX2 = 5
X2 = 1:nX2
M = [0.8 0.2 0.0 0.0 0.0;
	 0.2 0.6 0.2 0.0 0.0;
	 0.0 0.2 0.6 0.2 0.0;
	 0.0 0.0 0.2 0.6 0.2;
	 0.0 0.0 0.0 0.2 0.8;]

S = States(S1,
           CommonState(X2, M))

Z = [zeros(nX2*2, 3), # don't buy
             [-ones(nX2) log.(X2) -ones(nX2); # buy
			  -ones(nX2) log.(X2)  zeros(nX2)]]

U = LinearUtility(Z, bdisc, copy(truepar))

M, N, T = 2, 100, 1000
D = simulate(U, S, M, N, T)

n_transitions = zeros(nX2, nX2)
for n = 1:T*N-1
	if D.id[n] == D.id[n+1,1]
		n_transitions[D.xs[n,2], D.xs[n+1,2]] += 1
	end
end

estimated_F2 = n_transitions./sum(n_transitions,2)

estimated_S =  States(S1, CommonState(X2, estimated_F2))

U = LinearUtility(Z, bdisc, copy(truepar)*0);
am_nfxp = fit_nfxp(U, estimated_S, D);
estimated_S =  States(S1, CommonState(X2, estimated_F2));
U = LinearUtility(Z, bdisc, copy(truepar)*0);
am_npl = fit_npl(U, estimated_S, D);
@time 1+1
U = LinearUtility(Z, bdisc, copy(truepar)*0)
@time am_nfxp = fit_nfxp(U, estimated_S, D)
estimated_S =  States(S1, CommonState(X2, estimated_F2));
U = LinearUtility(Z, bdisc, copy(truepar)*0);
@time am_npl = fit_npl(U, estimated_S, D)
estimated_S =  States(S1, CommonState(X2, estimated_F2));
U = LinearUtility(Z, bdisc, copy(truepar)*0);
@profile am_npl = fit_npl(U, estimated_S, D)


text:
# functions for Gibbs updates of non-SMC components
using Distributions

function calculatePredictiveParams(s_state::SamplerState,B_sum::Float64,K_sum::Float64,n_obs::Int64)::Tuple{Float64,Float64,Float64,Float64}

    ap_alpha = (s_state.a_α[1] + B_sum)::Float64 # a in BetaBern(a,b)
    bp_alpha = (s_state.b_α[1] + n_obs - B_sum)::Float64 # b in BetaBern(a,b)
    ap_lambda = (s_state.a_λ[1] + K_sum)::Float64 # r in NB(r,p)
    bp_lambda = (1/(s_state.b_λ[1] + n_obs + 1))::Float64 # p in NB(r,p)

    return ap_alpha,bp_alpha,ap_lambda,bp_lambda

end

function updateBandK!(s_state::SamplerState,particle_container::Array{Array{ParticleState,1},1},
                      L::Array{Float64,2},
                      eig_pgf::Array{Float64,1},
                      lp_b::Array{Float64,1},
                      p_b::Array{Float64,1},
                      lp_k::Array{Float64,1},
                      p_k::Array{Float64,1}
                      )

  T = s_state.ne_data[1]
  particle_path = s_state.particle_path
  k_range = s_state.k_range
  b_range = s_state.b_range
  num_k = length(k_range)

  # sufficient statistics
  # assert( s_state.K[1] == s_state.B[1] == zero(Int64) )
  K_sum = convert(Float64,sum(s_state.K))
  B_sum = convert(Float64,sum(s_state.B))

  t_order = randperm(T)

  # pre-allocate some things
  edge = zeros(Int64,1,2)
  rwp = zeros(Float64,1,2)
  frwp = zeros(Float64,1)
  # eig_pgf = zeros(Float64,s_state.nv_data[1])
  root_vtx = zeros(Int64,1)
  rt = zeros(Int64,1,2)
  w = zeros(Float64,1,2)
  nbd = zeros(Int64,size(particle_container[1][1].nbd_list,1))

  # p_k = zeros(Float64,num_k)
  # lp_k = zeros(Float64,num_k)
  # p_b = zeros(Float64,2)
  # lp_b = zeros(Float64,2)

  for t in t_order

    if t>1 # skip t=1, deterministically inserts an edge with a new vertex

      resetArray!(lp_k)
      resetArray!(lp_b)
      resetArray!(p_k)
      resetArray!(p_b)
      resetArray!(nbd)
      resetArray!(w)

      # update sufficient statistics, predictive parameters
      K_sum += -s_state.K[t]
      B_sum += -s_state.B[t]
      ap_alpha,bp_alpha,ap_lambda,bp_lambda = calculatePredictiveParams(s_state,B_sum,K_sum,T-2)

      particle_tm1 = particle_container[particle_path[t-1]][t-1]
      particle_t = particle_container[particle_path[t]][t]
      nv_tm1 = particle_tm1.n_vertices[1]

      if particle_t.new_vertex[1]

        ## update B_t
        edge[:] = particle_t.vertex_unmap[particle_t.edge_list[t,:]]
        root_vtx[:] = find( particle_tm1.vertex_map[edge] .> 0)::Array{Int64,1}
        rt[:] = edge[root_vtx[1]]
        # assert(countnz(root_vtx)==1)
        if !particle_tm1.has_eigensystem[1]
          updateEigenSystem!(L,nv_tm1,particle_tm1)
          L[:] = zero(Float64)
        end
        # calculate negative binomial pgf on eigenvalues
        nbPred!(eig_pgf,nv_tm1,ap_lambda,bp_lambda,particle_tm1.eig_vals)
        # get one ball
        nbd[:] = particle_tm1.nbd_list[:,edge[root_vtx[1]]]
        # calculate fruitless random walk probs
        randomWalkProbs!(w,rt,nv_tm1,particle_tm1.degrees,eig_pgf,particle_tm1.eig_vecs) # r.w. ends where it started
        frwp[1] = (w[1] + randomWalkProbs(edge[root_vtx[1]],nbd,nv_tm1,particle_tm1.degrees,eig_pgf,particle_tm1.eig_vecs))::Float64
        # println("b " * string(t) * " " * string(frwp[1]))

        lp_b[2] = (log(ap_alpha) - log(ap_alpha + bp_alpha))::Float64
        lp_b[1] = (log(bp_alpha) - log(ap_alpha + bp_alpha) + log(frwp[1]))::Float64
        logSumExpWeightsNorm!(p_b,lp_b)
        s_state.B[t] = wsample(b_range,p_b)

        # update K_t
        for k in k_range
          for i = 1:nv_tm1
            k==zero(Int64) ? eig_pgf[i] = (one(Float64) - particle_tm1.eig_vals[i])::Float64 : eig_pgf[i] *= (one(Float64) - particle_tm1.eig_vals[i])::Float64
          end
          randomWalkProbs!(w,rt,nv_tm1,particle_tm1.degrees,eig_pgf,particle_tm1.eig_vecs)
          frwp[1] = (w[1] +  randomWalkProbs(edge[root_vtx[1]],nbd,nv_tm1,particle_tm1.degrees,eig_pgf,particle_tm1.eig_vecs))::Float64
          # println("k " * string(t) * " " * string(k) * " " * string(frwp[1]))  # NEGATIVE???
          # if isnan(frwp[1])
          #   println("nan frwp " * string(t) * " root " * string(edge[root_vtx[1]]) * " nbd " * string(nbd))
          # end
          lp_k[k+1] = (lgamma(k + ap_lambda) - lgamma(k + one(Int64)) - lgamma(ap_lambda) + k*log(bp_lambda) +
                        log( ap_alpha/(ap_alpha + bp_alpha) +
                            (bp_alpha/(ap_alpha + bp_alpha))*frwp[1] ))::Float64
        end
        logSumExpWeightsNorm!(p_k,lp_k)
        s_state.K[t] = wsample(k_range,p_k) # range is 0:k_trunc

        # update sufficient statistics
        B_sum += s_state.B[t]
        K_sum += s_state.K[t]

      else # new edge w/out vertex
        s_state.B[t] = zero(Int64)

        # update K_t
        if !particle_tm1.has_eigensystem[1]
          updateEigenSystem!(L,nv_tm1,particle_tm1)
          L[:] = zero(Float64)
        end

        for k in k_range # 0:k_trunc
          for i = 1:nv_tm1
            k==zero(Int64) ? eig_pgf[i] = (one(Float64) - particle_tm1.eig_vals[i])::Float64 : eig_pgf[i] *= (one(Float64) - particle_tm1.eig_vals[i])::Float64
          end

          edge[:] = particle_t.vertex_unmap[particle_t.edge_list[t,:]]
          randomWalkProbs!(rwp,edge,nv_tm1,particle_tm1.degrees,eig_pgf,particle_tm1.eig_vecs)

          lp_k[k+1] = (lgamma(ap_lambda + k) - lgamma(k + one(Int64)) - lgamma(ap_lambda) + k*log(bp_lambda) +
                      log(bp_alpha) - log(ap_alpha + bp_alpha) +
                      log( particle_tm1.degrees[edge[1]]*rwp[1] + particle_tm1.degrees[edge[2]]*rwp[2] ))::Float64
        end
        logSumExpWeightsNorm!(p_k,lp_k)
        s_state.K[t] = wsample(k_range,p_k) # range is 0:k_trunc

        # update sufficient statistics
        B_sum += s_state.B[t]
        K_sum += s_state.K[t]

      end

    end

  end

end

function updateAlphaAndLambda!(s_state::SamplerState)

  B_sum = sum(s_state.B)
  K_sum = sum(s_state.K)
  s_state.α[:] = rand(Beta(s_state.a_α[1] + B_sum, s_state.b_α[1] + s_state.ne_data[1] - 1 - B_sum))::Float64
  s_state.λ[:] = rand(Gamma(s_state.a_λ[1] + K_sum, 1/(s_state.b_λ[1] + s_state.ne_data[1] - 1)))::Float64

end


text:
# using StatsBase, DataFrames, JLD2, FileIO

# const DIRS = [:left, :up, :right, :down]
# const VALS = vcat(0,[(2).^(1:13)...])

# rand2_1() = rand() < 0.1 ? 2 : 1
# rand4() = Int8((rand(UInt8) % 4) + 1)

# if false
#     # x = [rand4() for i =1:1_000_000]
#     # using StatsBase
#     # countmap(x)

#     # move left
#     xinc=1
#     xstart=1
#     xend=4

#     # move right
#     xinc=-1
#     xstart=4
#     xend=1
# end
# # a function to simulate the move and return a reward
# function move!(x, xinc, xstart, xend)
#     reward = 0
#     @inbounds for i = xstart:xinc:xend # for each row move the left most piece first #if move_row = 1 then i control the row
#         if x[i] != 0 # if the position is occupied by a number move it
#             # firstly look "behind" to see if there is a number that is the same
#             # this is to deal better with situations like 2 2 4 4
#             @inbounds for k = i+xinc:xinc:xend
#                 if x[k] != 0
#                     if x[k] == x[i]
#                         x[i] += 1
#                         reward += 2^x[i]
#                         x[k] = 0
#                     end
#                     break;
#                 end
#             end
#             # now place it in the first empty slot
#             @inbounds for k = xstart:xinc:i
#                 if x[k] == 0
#                     x[k] = x[i]
#                     x[i] = 0
#                 end
#             end
#         end
#     end
#     (x, reward)
# end

# function move_left!(x)
#     move!(x, 1, 1, 4)
# end

# function move_right!(x)
#     move!(x, -1, 4, 1)
# end

# function move!(grid::Array{T,2}, direction) where T <: Integer
#     reward::Int16 = zero(Int16)
#     if direction == :left
#         for j = 1:4
#             #grid[j,:] .= move_left!(grid[j,:])
#             (tmp, new_reward) = move_left!(@view grid[j,:])
#             reward += new_reward
#         end
#     elseif direction == :right
#         for j = 1:4
#             #grid[j,:] .= move_right!(grid[j,:])
#             (tmp, new_reward) = move_right!(@view grid[j,:])
#             reward += new_reward
#         end
#     elseif direction == :up
#         for j = 1:4
#             #grid[:,j] .= move_left!(grid[:,j])
#             (tmp, new_reward) = move_left!(@view grid[:,j])
#             reward += new_reward
#         end
#     else
#         for j = 1:4
#             #grid[:,j] .= move_right!(grid[:,j])
#             (tmp, new_reward) = move_right!(@view grid[:,j])
#             reward += new_reward
#         end
#     end
#     (grid, reward)
# end

# function has_valid_move(x, direction)
#     if direction == :left
#         xstart = 1
#         xinc = 1
#         xend = 4
#     else
#         xstart = 4
#         xinc = -1
#         xend = 1
#     end
#     ((x[xstart] == 0)  & any(x[xstart+xinc:xinc:xend] .> 0)) |
#     ((x[xstart+xinc] == 0)  & any(x[xstart+2xinc:xinc:xend] .> 0)) |
#     ((x[xstart+2xinc] == 0)  & (x[xend] > 0)) |
#     ((x[xstart] > 0) & (x[xstart] == x[xstart + xinc])) |
#     ((x[xstart+xinc] > 0) & (x[xstart+xinc] == x[xstart+2xinc])) |
#     ((x[xstart+2xinc] > 0) & (x[xstart+2xinc] == x[xend]))
# end

# function valid_moves(grid::Array{T,2}) where T
#     left_has = any([has_valid_move(@view(grid[i,:]), :left) for i=1:4])
#     right_has = any([has_valid_move(@view(grid[i,:]), :right) for i=1:4])
#     up_has = any([has_valid_move(@view(grid[:,i]), :left) for i=1:4])
#     down_has = any([has_valid_move(@view(grid[:,i]), :right) for i=1:4])
#     [:left, :right, :up, :down][[left_has, right_has, up_has, down_has]]
# end

# function simulate_move!(grid)
#     directions = sample(DIRS, 4 , replace = false)
#     for i in 1:3
#         d1 = directions[i]
#         (grid, ok, cart, two_or_four, reward) = simulate_move!(grid, d1)
#         if ok
#            return  (grid, true, d1, cart, two_or_four, reward)
#         end
#     end
#     (grid, ok, cart, two_or_four, reward) = simulate_move!(grid, directions[4])
#     (grid, ok, directions[4], cart, two_or_four, reward)
# end

# function simulate_move!(grid, directions::Array{Symbol,1})
#     for i in 1:3
#         d1 = directions[i]
#         (grid, ok, cart, two_or_four, reward) = simulate_move!(grid, d1)
#         if ok
#            return  (grid, true, d1, cart, two_or_four, reward)
#         end
#     end
#     (grid, ok, cart, two_or_four, reward) = simulate_move!(grid, directions[4])
#     (grid, ok, directions[4], cart, two_or_four, reward)
# end

# # assume no need to check for validate moves
# function simulate_move!(grid, direction::Symbol)
#     tmp_grid = copy(grid)
#     (grid, reward) = move!(grid, direction)
#     if all(tmp_grid .== grid)
#         return (grid, false, CartesianIndex{2}(-1,-1), -1, 0)
#     else
#         cart = rand(findall(grid .== 0)) # randomly choose one empty slot
#         one_or_two = rand2_1()
#         grid[cart] = one_or_two
#         return (grid, true, cart, one_or_two, reward)
#     end
# end

# function simulate_game!(grid, lim)
#     init_grid = copy(grid)
#     seq = Symbol[]::Array{Symbol,1}
#     cartarr = CartesianIndex{2}[]
#     one_or_two_arr = Int8[]
#     reward_vec = Int16[]

#     ok = true

#     while ok
#         (grid, ok1, move, cart, one_or_two, new_reward) = simulate_move!(grid)
#         ok = ok1 & all(grid .< lim)
#         push!(seq, move)
#         push!(cartarr, cart)
#         push!(one_or_two_arr, one_or_two)
#         push!(reward_vec, new_reward)
#     end
#     (init_grid, grid, seq, cartarr, one_or_two_arr, reward_vec)
# end

# function simulate_game()
#     init = init_game()
#     simulate_game!(init, Inf)
# end

# function init_game()
#     grid = zeros(Int8,4,4)
#     grid[rand(1:4),rand(1:4)] = rand2_1()
#     grid[rand(1:4),rand(1:4)] = rand2_1()
#     grid
# end

# function see_which_move_is_better(grid)
#     score = sum(grid)
#     grid_try = similar(grid)
#     res_dict = Dict{Symbol, Float64}()
#     for direction in [:left, :right, :up, :down]
#          res_dict[direction] = StatsBase.mean([see_which_move_is_better!(grid, grid_try, direction) for j=1:100])
#     end
#     res_dict
# end

# function see_which_move_is_better!(grid, grid_try, direction)
#     grid_try .= grid
#     simulate_one!(grid_try, direction)
#     simulate_game!(grid_try)
#     sum(grid_try)
# end

# # simulate n moves ahead
# function simulate_moves_ahead(grid)
#     grid_try = similar(grid)

#     res_dict = Dict{NTuple{9, Symbol}, Float64}()
#     for moves in Iterators.product([DIRS for i = 1:9]...)
#         grid_try .= grid
#         simulate_one!.([grid_try], moves)
#         res_dict[moves] = sum(grid_try)
#     end
#     res_dict
# end

# function simulate_n_games_get_score(n)
#     reward_vec = Int64[]
#     maxcell = Int64[]
#     for i=1:n
#         init = init_game()
#         (init, fnl, seq, cart, tf, reward) = simulate_game!(init);
#         push!(reward_vec, sum(reward))
#         push!(maxcell, maximum(fnl))
#     end
#     (reward_vec,maxcell)
# end

# function simulate_game!(grid, nn, ϵ)
#     init_grid = copy(grid)
#     seq = Symbol[]::Array{Symbol,1}
#     cartarr = CartesianIndex{2}[]
#     two_or_four_arr = Int64[]
#     reward_vec = Int64[]

#     vm = valid_moves(grid)
#     ok = length(vm) > 0

#     while ok
#         if rand() < ϵ
#             move = rand(vm)
#         else
#             prob = Tracker.data.(nn.(grid_move_to_feature.([grid], vm)))
#             x = [aa[1] for aa in prob]
#             x .= x .- minimum(x) .+ 1
#             #x .= exp.(x)
#             x .= x.^2
#             move = sample(vm, Weights(x), 1)[1]
#         end
#         (grid, cart, two_or_four, new_reward) = simulate_one!(grid, move)
#         push!(seq, move)
#         push!(cartarr, cart)
#         push!(two_or_four_arr, two_or_four)
#         push!(reward_vec, new_reward)
#         vm = valid_moves(grid)
#         ok = length(vm) > 0
#     end

#     push!(seq, :stopped)
#     push!(cartarr, CartesianIndex{2}())
#     push!(two_or_four_arr, 0)
#     push!(reward_vec, 0.0)

#     (init_grid, grid, seq, cartarr, two_or_four_arr, reward_vec)
# end

# function simulate_one_game_and_train!(model, opt, ϵ)
#     grid = init_game()
#     (init, fnl, seq, cart, tf, reward) = simulate_game!(grid, model, ϵ)
#     x = generate_x_based_on_episode(init, seq, cart, tf);
#     y = discounted_rewards(reward, 0.99)[1:end-1];
#     xy = [(x[:,j], y[j]) for j in 1:size(x)[2]]
#     Flux.train!(loss, xy, opt);
#     nothing
#     #lxy= loss(x,y)
#     #println(lxy)
#     #lxy
# end

# function train_model_n(n, model, opt, ϵ)
#     for i=1:n
#         ϵ = ϵ * 0.999
#         simulate_one_game_and_train!(model, opt, ϵ)
#         if i % 100 == 0
#             grid = init_game()
#             (init, fnl, seq, cart, tf, reward) = simulate_game!(grid, model, ϵ)
#             x = generate_x_based_on_episode(init, seq, cart, tf);
#             y = discounted_rewards(reward, 0.99)[1:end-1];
#             #display(model(x)[[1,end]])
#     #         @show y[[1,end]]
#     #         display((ϵ, loss(x,y)))
#             println(maximum(fnl))
#         end
#     end
#     @save "model.jld2" model
#     ϵ
# end

# function grid_move_to_feature(grid, move)
#     input = vec(Flux.onehotbatch(vec(grid),VALS))
#     # create the data
#     vcat(input, Flux.onehot(move, DIRS))
# end

# function grid_to_feature(grid)
#     vec(Flux.onehotbatch(vec(grid),VALS))
# end

# function generate_x_based_on_episode(init, seq, cart, tf)
#     grid = copy(init)

#     res = Array{Bool, 2}(undef, 14*16+4, length(seq)-1)
#     res[:, 1] .= grid_move_to_feature(grid, seq[1])

#     for i = 2:length(seq)-1
#         #display(grid)
#         move!(grid, seq[i])
#         grid[cart[i]] = tf[i]
#         res[:,i] .= grid_move_to_feature(grid, seq[i])
#     end
#     res
# end

# function sim_one_move_then_end(grid, move)
#     grid_copy = copy(grid)
#     simulate_one!(grid_copy, move)
#     simulate_game!(grid_copy)
#     grid_copy
# end

# function sim_one_move_then_end_n_times(grid, move, n)
#     [sum(sim_one_move_then_end(grid, move)) for i = 1:n]
# end

# function sim_iterate_one_move(grid, vm, n)
#     [sim_one_move_then_end_n_times(grid, move, n) for move in vm]
# end

# # replay the game from beginning to second last move
# function sim_seq!(init, seq, cart, tf)
#     for (s, c, t) in zip(seq, cart, tf)
#         move!(init, s)
#         init[c] = t
#     end
#     init
# end

# function replay_move!(agrid, move, cart, tf)
#     move!(agrid, move)
#     agrid[cart] = tf
#     agrid
# end

# function discounted_policy_rewards(reward_vec, λ)
#     l = length(reward_vec)
#     res = Vector{Float64}(undef, l)
#     res[end] = reward_vec[end] # the last move is always of no value
#     for l1 in l-1:-1:1
#         res[l1] = reward_vec[l1] + res[l1+1]*λ
#     end

#     res
# end

# function discounted_rewards(reward_vec, λ)
#     l = length(reward_vec)
#     res = Vector{Float64}(undef, l)
#     res[l] = 0 # the last move is always of no value
#     for l1 in l-1:-1:1
#         res[l1] = reward_vec[l1] + res[l1+1]*λ
#     end

#     res
# end

# function avg_disc_reward(agrid, n, disc)
#     res = Array{Float64, 1}(undef, n)
#     for i=1:n
#       grid_copy = copy(agrid)
#       (_, _, _, _, _, reward_vec) = simulate_game!(grid_copy)
#       res[i] = discounted_rewards(reward_vec, disc)[1]
#     end
#     StatsBase.mean(res)
# end

# function avg_rewards_next_move(grid, move, n)
#   res = Array{Float64, 1}(undef, n)
#   for i=1:n
#     grid_copy = copy(grid)
#     simulate_one!(grid_copy, move)
#     (_, _, _, _, _, reward_vec) = simulate_game!(grid_copy)
#     res[i] = sum(reward_vec)
#   end
#   StatsBase.mean(res)
# end

# function avg_rewards_next_move(grid, n)
#   vm = valid_moves(grid)
#   res = avg_rewards_next_move.([grid], vm, n)
#   Dict(zip(vm, res))
# end

# function rotatedir(move)
#   if move == :left
#     return (:left, :down, :right, :up, :up, :left, :down, :right)
#   elseif move == :down
#     return (:down, :right, :up, :left, :right, :up, :left, :down)
#   elseif move == :right
#     return (:right, :up, :left, :down, :down, :right, :up, :left)
#   elseif move == :up
#     return (:up, :left, :down, :right, :left, :down, :right, :up)
#   else
#     return (move, move, move, move, move, move, move, move)
#   end
# end

# function rotations(agrid::T, move) where T
#   bgrid::T = transpose(agrid)
#     ((agrid,
#     agrid |> rotl90,
#     agrid |> rot180,
#     agrid |> rotr90,
#     bgrid,
#     bgrid |> rotl90,
#     bgrid |> rot180,
#     bgrid |> rotr90
#     ), rotatedir(move))
#     #((agrid,), (move,))
# end

# struct IdentitySkip
#    inner
#    activation
# end

# (m::IdentitySkip)(x) = m.activation.(m.inner(x) .+ x)


text:
using MCObservables
using Base.Test

include("inform.jl")

include("SimpleObservable.jl")
include("BinningObservable.jl")
include("Jackknife.jl")
include("VectorObservable.jl")


text:
using KLU
using KLU: increment!, KLUITypes, decrement, klu!, KLUFactorization,
klu_analyze!, klu_factor!
using Test
using SparseArrays: SparseMatrixCSC, sparse, nnz
using LinearAlgebra

@testset "KLU Wrappers" begin
    Ap = increment!([0,4,1,1,2,2,0,1,2,3,4,4])
    Ai = increment!([0,4,0,2,1,2,1,4,3,2,1,2])
    Ax = [2.,1.,3.,4.,-1.,-3.,3.,6.,2.,1.,4.,2.]
    A0 = sparse(Ap, Ai, Ax)
    A1 = sparse(increment!([0,4,1,1,2,2,0,1,2,3,4,4]),
                increment!([0,4,0,2,1,2,1,4,3,2,1,2]),
                [2.,1.,3.,4.,-1.,-3.,3.,9.,2.,1.,4.,2.], 5, 5)
    @testset "Core functionality for $Tv elements" for Tv in (Float64, ComplexF64)
        @testset "Core functionality for $Ti indices" for Ti ∈ Base.uniontypes(KLUITypes)
            A = convert(SparseMatrixCSC{Tv, Ti}, A0)
            # test the raw vector construction method.
            klua = klu(size(A, 1), decrement(A.colptr), decrement(A.rowval), A.nzval)
            @test klua.F == klu(A).F
            @test nnz(klua) == 18
            R = Diagonal(Tv == ComplexF64 ? complex.(klua.Rs) : klua.Rs)
            @test R \ A[klua.p, klua.q] ≈ (klua.L * klua.U + klua.F)

            b = [8., 45., -3., 3., 19.]
            x = klua \ b
            @test x ≈ float([1:5;])
            @test A*x ≈ b

            z = complex.(b)
            x = ldiv!(klua, z)
            @test x ≈ float([1:5;])
            @test z === x
            # Can't match UMFPACK's ldiv!(<OUTPUT>, <KLU>, <INPUT>)
            # since klu_solve(<KLU>, <INPUT>) modifies <INPUT>, and has no field for <OUTPUT>.
            @test A*x ≈ b

            b = [8., 20., 13., 6., 17.]
            x = klua'\b
            @test x ≈ float([1:5;])

            @test A'*x ≈ b
            z = complex.(b)
            x = ldiv!(adjoint(klua), z)
            @test x ≈ float([1:5;])
            @test x === z

            @test A'*x ≈ b
            x = transpose(klua) \ b
            @test x ≈ float([1:5;])
            @test transpose(A) * x ≈ b

            x = ldiv!(transpose(klua), complex.(b))
            @test x ≈ float([1:5;])
            @test transpose(A) * x ≈ b

            @inferred klua\fill(1, size(A, 2))

            @testset "Permutation vectors" begin
                #Just to test this works, we'll use the existing permutation vectors.
                @test klu_analyze!(klua, klua.p, klua.q).common.status == 0
                klu_factor!(klua)
                x = klua \ b
                @test A*x ≈ b
            end
            @testset "Utility functions" begin
                K = KLUFactorization(A);
                @test size(K) == (5, 5)
                @test size(K, 3) == 1
                @test_throws ArgumentError K.symbolic
                @test_throws ArgumentError K.numeric
                klu_analyze!(K);
                @test K.symbolic.nz == 12
                @test K.symbolic.nzoff == 4
                klu_factor!(K)
                @test K.numeric.lnz == 7 == K.numeric.unz
                @test K.nblocks == 3
                @test propertynames(K) == (:lnz, :unz, :nzoff, :L, :U, :F, :q, :p, :Rs, :symbolic, :numeric,)
                @test_throws ArgumentError size(K, -1)
                @test size(K, 3) == 1
            end
            @testset "Refactorization" begin
                B = convert(SparseMatrixCSC{Tv, Ti}, A1)
                b = Tv[8., 45., -3., 3., 19.]
                F = klu(A)
                klu!(F, B)
                @test F\b ≈ B\b ≈ Matrix(B)\b
                @test klu!(F, B.nzval)\b ≈ B\b ≈ Matrix(B)\b #test just supplying nzval for recompute
            end

            @testset "Singular matrix" begin
                S = sparse([1 2; 0 0])
                S = convert(SparseMatrixCSC{Tv, Ti}, S)
                @test_throws SingularException klu(S)
            end
        end
    end
    @testset "REPL printing of KLU" begin
        A = sparse([1, 2], [1, 2], Float64[1.0, 1.0])
        F = klu(A)
        facstring = sprint((t, s) -> show(t, "text/plain", s), F)
        lstring = sprint((t, s) -> show(t, "text/plain", s), F.L)
        ustring = sprint((t, s) -> show(t, "text/plain", s), F.U)
        fstring = sprint((t, s) -> show(t, "text/plain", s), F.F)
        @test facstring == "$(summary(F))\nL factor:\n$lstring\nU factor:\n$ustring\nF factor:\n$fstring"
    end
end


text:
export alloc_opt!, barrier_noop!, gc_invariant_verifier!, lower_exc_handlers!,
       combine_mul_add!, multi_versioning!, propagate_julia_addrsp!, lower_ptls!,
       lower_simdloop!, late_lower_gc_frame!, final_lower_gc!, remove_julia_addrspaces!,
       demote_float16!, remove_ni!, julia_licm!, cpu_features!

alloc_opt!(pm::PassManager) = API.LLVMAddAllocOptPass(pm)

barrier_noop!(pm::PassManager) = API.LLVMAddBarrierNoopPass(pm)

gc_invariant_verifier!(pm::PassManager, strong::Bool=false) =
    API.LLVMAddGCInvariantVerifierPass(pm, convert(Bool, strong))

lower_exc_handlers!(pm::PassManager) = API.LLVMAddLowerExcHandlersPass(pm)

combine_mul_add!(pm::PassManager) = API.LLVMAddCombineMulAddPass(pm)

multi_versioning!(pm::PassManager) = API.LLVMAddMultiVersioningPass(pm)

propagate_julia_addrsp!(pm::PassManager) = API.LLVMAddPropagateJuliaAddrspaces(pm)

lower_ptls!(pm::PassManager, imaging_mode::Bool=false) =
    API.LLVMAddLowerPTLSPass(pm, convert(Bool, imaging_mode))

lower_simdloop!(pm::PassManager) = API.LLVMAddLowerSimdLoopPass(pm)

late_lower_gc_frame!(pm::PassManager) = API.LLVMAddLateLowerGCFramePass(pm)

final_lower_gc!(pm::PassManager) = API.LLVMAddFinalLowerGCPass(pm)

remove_julia_addrspaces!(pm::PassManager) = API.LLVMAddRemoveJuliaAddrspacesPass(pm)

demote_float16!(pm::PassManager) = API.LLVMAddDemoteFloat16Pass(pm)

remove_ni!(pm::PassManager) = API.LLVMAddRemoveNIPass(pm)

julia_licm!(pm::PassManager) = API.LLVMAddJuliaLICMPass(pm)

if VERSION >= v"1.8.0-DEV.1120"
    cpu_features!(pm::PassManager) = API.LLVMAddCPUFeaturesPass(pm)
else
    cpu_features!(pm::PassManager) = nothing
end


text:
using Test
using Random
using MPI

MPI.Init()

function bcast_array(A, root)
    comm = MPI.COMM_WORLD

    if MPI.Comm_rank(comm) == root
        B = copy(A)
    else
        B = similar(A)
    end

    MPI.Bcast!(B, root, comm)
    B
end

root = 0

if VERSION >= v"0.7.0-rc1"
    Random.seed!(17)
else
    srand(17)
end

matsize = (17,17)
for typ in Base.uniontypes(MPI.MPIDatatype)
    A = rand(typ, matsize...)
    @test bcast_array(A, root) == A
end

# Char
A = ['s', 't', 'a', 'r', ' ', 'w', 'a', 'r', 's']
@test bcast_array(A, root) == A

comm = MPI.COMM_WORLD

g = x -> x^2 + 2x - 1
if MPI.Comm_rank(comm) == root
    f = g
else
    f = nothing
end
f = MPI.bcast(f, root, comm)
@test f(3) == g(3)
@test f(5) == g(5)
@test f(7) == g(7)

A = Dict("foo" => "bar")
if MPI.Comm_rank(comm) == root
    B = A
else
    B = nothing
end
B = MPI.bcast(B, root, comm)
@test B["foo"] == "bar"

MPI.Finalize()