text:
# Flip image pixels vertically
@static if use_VideoIO
function vflip!(A)
    nrows, ncols = size(A)
    nrp1 = nrows + 1
    for col = 1:ncols
        for row = 1:div(nrows,  2)
            t = A[nrp1-row, col]
            A[nrp1-row, col] = A[row, col]
            A[row, col] = t
        end
    end

    return A
end
end

# Flip image pixels horizontally
function hflip!(A)
    nrows, ncols = size(A)
    ncp1 = ncols + 1
    for row = 1:nrows
        for col = 1:div(ncols, 2)
            t = A[row, ncp1-col]
            A[row, ncp1-col] = A[row, col]
            A[row, col] = t
        end
    end

    return A
end

function alignscale(s::mjSim)
    s.cam[].lookat = s.m.stat.center
    s.cam[].distance = 1.5*s.m.stat.extent

    # set to free camera
    s.cam[].type = MJCore.mjCAMERA_FREE
end

function str2vec(s::String, len::Int)
    str = zeros(UInt8, len)
    str[1:length(s)] = codeunits(s)
    return str
end

#=
function reloadmodel(s::mjSim)
    width, height = GLFW.GetFramebufferSize(s.window)

    loadmodel(s.modelfile, width, height)
end
=#

function print_state(s::mjSim)
    #println(s.d.qpos)
    c = s.robot.command
    println("== Robot state ==")
    println("velocity: $(round(c.horizontal_velocity[1], digits=2))")
    println("height:   $(round(c.height,    digits=2))")
    println("yaw:      $(round(c.yaw_rate,  digits=2))")
    println("pitch:    $(round(c.pitch,     digits=2))")
    println("roll:     $(round(c.roll,      digits=2))")
    #println("=================")
end


text:
using DSGE
using HDF5
path = dirname(@__FILE__)

# Test in model optimization
custom_settings = Dict{Symbol, Setting}(
    :date_forecast_start  => Setting(:date_forecast_start, quartertodate("2015-Q4")))
m = AnSchorfheide(custom_settings = custom_settings, testing = true)

file = "$path/../reference/optimize.h5"
x0 = h5read(file, "params")
data = h5read(file, "data")'
minimizer = h5read(file, "minimizer")
minimum = h5read(file, "minimum")
H_expected = h5read(file, "H")

# See src/estimate/estimate.jl
update!(m, x0)
n_iterations = 3

x0 = Float64[p.value for p in m.parameters]

out, H = optimize!(m, data; iterations=n_iterations)

@testset "Check optimize minimizers are the same" begin
    @test @test_matrix_approx_eq minimizer out.minimizer
    @test minimum ≈ out.minimum atol=5e-7
    @test @test_matrix_approx_eq H_expected H
end

nothing

ext: jl


text:
using Tokamak: @tk, cpu, dagger
using Base.Test

@testset "Tokamak" begin

@tk diag(A)[i] = A[i,i]
@tk outer(x,y)[i,j] = x[i]*y[j]
@tk trans(A)[i,j] = A[j,i]
@tk add(A,B)[i] = A[i] + B[i]

@testset "CPU" begin

A = rand(5,5)
B = rand(5,5)
x = rand(5)
y = rand(5)

cpu_diag = cpu(diag, Matrix{Float64})
@test cpu_diag(A) == Base.diag(A)

cpu_outer = cpu(outer, Vector{Float64}, Vector{Float64})
@test cpu_outer(x,y) == x*y'

cpu_trans = cpu(trans, Matrix{Float64})
@test cpu_trans(A) == A'

cpu_add = cpu(add, Vector{Float64}, Vector{Float64})
@test cpu_add(x,y) == x+y

end

@testset "Dagger" begin

using Dagger

A = distribute(rand(100,100), Blocks(10,10))
B = distribute(rand(100,100), Blocks(10,10))
x = distribute(1:100, Blocks(10))
y = distribute(1:100, Blocks(10))

dagger_diag = dagger(diag, Matrix{Float64})
@test dagger_diag(A) |> collect == Base.diag(collect(A))

dagger_outer = dagger(outer, Vector{Float64}, Vector{Float64})
@test dagger_outer(x,y) == collect(x)*collect(y)'

dagger_trans = dagger(trans, Matrix{Float64})
@test dagger_trans(A) == A'

dagger_add = dagger(add, Vector{Float64}, Vector{Float64})
@test dagger_add(x,y) == x+y

end

end

max_stars_repo_name: JuliaTagBot/Tokamak

max_issues_repo_name: JuliaTagBot/Tokamak

max_forks_repo_name: FluxML/Tokamak


text:
# Chapter 15. Structs and Objects

# export Point

# Composite Types
struct Point
    x
    y
end

# Exercise 17-2
function Point(x::Float64=0, y::Float64=0)
    Point(x,y)
end

# Mutable Structs
# https://benlauwens.github.io/ThinkJulia.jl/latest/book.html#_mutable_structs
mutable struct MPoint
    x
    y
end

# Rectangles
# https://benlauwens.github.io/ThinkJulia.jl/latest/book.html#_rectangles
"""
Represents a rectangle.

width: Width
height:
corner: Lower left corner
"""
struct Rectangle
    width
    height
    corner::MPoint
end

# Instances as Arguments
# https://benlauwens.github.io/ThinkJulia.jl/latest/book.html#_instances_as_arguments
function printpoint(p)
    println("($(p.x), $(p.y))")
end

# Exercise 15-1
function distancebetweenpoints(p, q)
    x2 = p.x^2 + q.x^2
    y2 = p.y^2 + q.y^2
    ret = sqrt(x2 + y2)

    return ret
end

function movepoint!(p, dx, dy)
    p.x += dx
    p.y += dy

    nothing
end

function moverectangle!(rect, dx, dy)
    movepoint!(rect.corner, dx, dy)
end

# Instances as Return Values
function findcenter(rect)
    Point(
    rect.corner.x + rect.width/2.,
    rect.corner.y + rect.height/2.)
end


text:
using Base.Test
using Cubature
include("../model/AccelerationControlledDifferentialDrive.jl")

let state = initial_state(0, 0, 0, 0, 0)
    control = [10, 10]
    elapsed_time = 2
    (x0, y0, theta0, wl0, wr0) = state
    (al, ar) = control
    r = 0.033
    L = 0.13
    ε_integration = 0.000001
    expected_position_x = x0 + hquadrature(t -> r/2 * (wl0 + al*t + wr0 + ar*t) * cos(theta0 + r/L*(wr0-wl0)*t + r/L*(ar-al)*t^2/2), 0, elapsed_time; abstol=ε_integration)[1] # TODO cubature doesn't seem to care about abstol
    expected_position_y = y0 + hquadrature(t -> r/2 * (wl0 + al*t + wr0 + ar*t) * sin(theta0 + r/L*(wr0-wl0)*t + r/L*(ar-al)*t^2/2), 0, elapsed_time; abstol=ε_integration)[1]
    expected_orientation = theta0 + r/L*(wr0-wl0)*elapsed_time + r/L*(ar-al)*elapsed_time^2/2
    ε_position = 0.001
	ε_orientation = 0.01

    state = new_state(state, control, elapsed_time)

    println()
	println("Constant acceleration linear move results")
    println("State:      ", state)
	println("Expected:   ", [expected_position_x expected_position_y expected_orientation])
	println("Error[mm]:  ", sqrt((expected_position_x - state[1])^2 + (expected_position_y - state[2])^2) * 1000)
	@test state[1] ≈ expected_position_x atol=ε_position
	@test state[2] ≈ expected_position_y atol=ε_position
	@test state[3] ≈ expected_orientation atol=ε_orientation
    println("Constant acceleration linear move pass")
end

let state = initial_state(1, 1.5, π/2, 5, 5)
    control = [-5, 0]
    elapsed_time = 0.1
    (x0, y0, theta0, wl0, wr0) = state
    (al, ar) = control
    r = 0.033
    L = 0.13
    ε_integration = 0.000001
    expected_position_x = x0 + hquadrature(t -> r/2 * (wl0 + al*t + wr0 + ar*t) * cos(theta0 + r/L*(wr0-wl0)*t + r/L*(ar-al)*t^2/2), 0, elapsed_time; abstol=ε_integration)[1]
    expected_position_y = y0 + hquadrature(t -> r/2 * (wl0 + al*t + wr0 + ar*t) * sin(theta0 + r/L*(wr0-wl0)*t + r/L*(ar-al)*t^2/2), 0, elapsed_time; abstol=ε_integration)[1]
    expected_orientation = theta0 + r/L*(wr0-wl0)*elapsed_time + r/L*(ar-al)*elapsed_time^2/2
    ε_position = 0.001
	ε_orientation = 0.01

    state = new_state(state, control, elapsed_time)

    println()
	println("Turn left from constant linear move results")
    println("State:      ", state)
	println("Expected:   ", [expected_position_x expected_position_y expected_orientation])
	println("Error[mm]:  ", sqrt((expected_position_x - state[1])^2 + (expected_position_y - state[2])^2) * 1000)
	@test state[1] ≈ expected_position_x atol=ε_position
	@test state[2] ≈ expected_position_y atol=ε_position
	@test state[3] ≈ expected_orientation atol=ε_orientation
    println("Turn left from constant linear move pass")
end

let state = initial_state(1, 1.5, π/2, 5, 5)
    control = [-5, 0]
    elapsed_time = 1.0
    cycle_time = 0.1
    (x0, y0, theta0, wl0, wr0) = state
    (al, ar) = control
    r = 0.033
    L = 0.13
    ε_integration = 0.000001
    expected_position_x = x0 + hquadrature(t -> r/2 * (wl0 + al*t + wr0 + ar*t) * cos(theta0 + r/L*(wr0-wl0)*t + r/L*(ar-al)*t^2/2), 0, elapsed_time; abstol=ε_integration)[1]
    expected_position_y = y0 + hquadrature(t -> r/2 * (wl0 + al*t + wr0 + ar*t) * sin(theta0 + r/L*(wr0-wl0)*t + r/L*(ar-al)*t^2/2), 0, elapsed_time; abstol=ε_integration)[1]
    expected_orientation = theta0 + r/L*(wr0-wl0)*elapsed_time + r/L*(ar-al)*elapsed_time^2/2
    ε_position = 0.001
	ε_orientation = 0.01

    for i in 1:10
        state = new_state(state, control, cycle_time)
    end

    println()
	println("Turn left from constant linear move multiple update results")
    println("State:      ", state)
	println("Expected:   ", [expected_position_x expected_position_y expected_orientation])
	println("Error[mm]:  ", sqrt((expected_position_x - state[1])^2 + (expected_position_y - state[2])^2) * 1000)
	@test state[1] ≈ expected_position_x atol=ε_position
	@test state[2] ≈ expected_position_y atol=ε_position
	@test state[3] ≈ expected_orientation atol=ε_orientation
    println("Turn left from constant linear move multiple update pass")
end

let (x0, y0, theta0, wl0, wr0) = initial_state(0, 0, 0, 0, 0)
    (al, ar) = [10, 10]
    elapsed_time = 2
    r = 0.033
    L = 0.13
    expected_x = 0.66
    expected_y = 0
    expected_theta = 0
    ε_integration = 0.000001
    ε_position = 0.001
	ε_orientation = 0.01

    calculated_theta = theta0 + r/L*(wr0-wl0)*elapsed_time + r/L*(ar-al)*elapsed_time^2/2
    calculated_x = x0 + hquadrature(t -> r/2 * (wl0 + al*t + wr0 + ar*t) * cos(theta0 + r/L*(wr0-wl0)*t + r/L*(ar-al)*t^2/2), 0, elapsed_time; abstol=ε_integration)[1]
    calculated_y = y0 + hquadrature(t -> r/2 * (wl0 + al*t + wr0 + ar*t) * sin(theta0 + r/L*(wr0-wl0)*t + r/L*(ar-al)*t^2/2), 0, elapsed_time; abstol=ε_integration)[1]

    println()
    println("Cubature test suit test constant acceleration linear move result")
    println("Calculated: ", [calculated_x calculated_y calculated_theta])
	println("Expected:   ", [expected_x expected_y expected_theta])
	println("Error[mm]:  ", sqrt((expected_x - calculated_x)^2 + (expected_y - calculated_y)^2) * 1000)
    @test calculated_x ≈ expected_x atol = ε_position
	@test calculated_y ≈ expected_y atol = ε_position
	@test calculated_theta ≈ expected_theta atol = ε_orientation
    println("Cubature test suit test constant acceleration linear move pass")
end

max_stars_repo_name: grgomrton/robot-model

max_issues_repo_name: grgomrton/robot-model

max_forks_repo_name: grgomrton/robot-model


text:
using Coverage
cov_res = process_folder()
Codecov.submit(cov_res)
Coveralls.submit(cov_res)


text:
import LinearAlgebra: checksquare, BlasInt
import Printf

## Todo
# 1. Add method for partial sweep cut
# 2. Add a "SymmetricMatrix" wrapper to avoid
#    constant symmetry checks.
# 3. Avoid creating the Laplacian matrix directly

"""
`_symeigs_smallest_arpack`
---
Compute eigenvalues and vectors using direct calls to the ARPACK wrappers
to get type-stability. This function works for symmetric matrices.

This function works on Float32 and Float64-valued sparse matrices.
It returns the smallest set of eigenvalues and vectors.

It only works on matrices with more than 21 rows and columns.
(Use a dense eigensolver for smaller problems.)

Functions
---------
- `(evals,evecs) = _symeigs_smallest_arpack(A::SparseMatrixCSC{V,Int},
                        nev::Int,tol::V,maxiter::Int, v0::Vector{V})`
Inputs
------
- `A`: the sparse matrix, must be symmetric
- `nev`: the number of eigenvectors requested
- `tol`: the relative tolerance of the eigenvalue computation
- `maxiter`: the maximum number of restarts
- `v0`: the initial vector for the Lanczos process

Example
-------
This is an internal function.
"""
function _symeigs_smallest_arpack(
            A::SparseMatrixCSC{V,Int},nev::Int,tol::V,maxiter::Int,
            v0::Vector{V}) where V

    n::Int = checksquare(A) # get the size
    @assert n >= 21

    # setup options
    mode = 1
    sym = true
    iscmplx = false
    bmat = String("I")
    ncv = min(max(2*nev,20),n-1)

    whichstr = String("SA")
    ritzvec = true
    sigma = 0.

    #TOL = Array{V}(undef,1)
    #TOL[1] = tol
    TOL = Ref(tol)
    lworkl = ncv*(ncv + 8)
    v = Array{V}(undef, n, ncv)
    workd = Array{V}(undef, 3*n)
    workl = Array{V}(undef, lworkl)
    resid = Array{V}(undef, n)

    resid[:] = v0[:]

    info = zeros(BlasInt, 1)
    info[1] = 1

    iparam = zeros(BlasInt, 11)
    ipntr = zeros(BlasInt, 11)
    ido = zeros(BlasInt, 1)

    iparam[1] = BlasInt(1)
    iparam[3] = BlasInt(maxiter)
    iparam[7] = BlasInt(mode)

    # this is a helpful indexing vector
    zernm1 = 0:(n-1)

    while true
        # This is the reverse communication step that ARPACK does
        # we need to extract the desired vector and multiply it by A
        # unless the code says to stop
        Arpack.saupd(
            ido, bmat, n, whichstr, nev, TOL, resid, ncv, v, n,
            iparam, ipntr, workd, workl, lworkl, info)

        load_idx = ipntr[1] .+ zernm1
        store_idx = ipntr[2] .+ zernm1

        x = workd[load_idx]

        if ido[1] == 1
            workd[store_idx] = A*x
        elseif ido[1] == 99
            break
        else
            error("unexpected ARPACK behavior")
        end
    end

    # Once the process terminates, we need to extract the
    # eigenvectors.

    # calls to eupd
    howmny = String("A")
    select = Array{BlasInt}(undef, ncv)

    d = Array{V}(undef, nev)
    sigmar = ones(V,1)*sigma
    ldv = n
    Arpack.seupd(ritzvec, howmny, select, d, v, ldv, sigmar,
        bmat, n, whichstr, nev, TOL, resid, ncv, v, ldv,
        iparam, ipntr, workd, workl, lworkl, info)
    if info[1] != 0
        error("unexpected ARPACK exception")
    end

    # Now we want to return them in sorted order (smallest first)
    p = sortperm(d)

    d = d[p]
    vectors = v[1:n,p]

    return (d,vectors,iparam[5])
end

"""
`fiedler_vector`
----------------

Compute the Fiedler vector associated with the normalized Laplacian
of the graph with adjacency matrix A.

This function works with Float32 and Float64-valued sparse
matrices and vectors.

It requires a symmetric input matrix representing
the adjacency matrix of an undirected graph. If the input
is a disconnected network then the result may not

The return vector is signed so that the number of positive
entries is at least the number of negative entries. This will
always give a unique, deterministic output
in the case of a non-repeated second eigenvalue.

Functions
---------
- `(v,lam2) = fiedler_vector(A::SparseMatrixCSC{V,Int})`
- `(v,lam2) = fiedler_vector(A::MatrixNetwork{V,Int})`

Inputs
------
- `A`: the adjacency matrix
- Optional Keywoard Inputs
  - `checksym`: ensure the matrix is symmetric
  - `tol`: the residual tolerance in the eigenvalue, eigenvector estimate.
    (This is an absolute error)
  - `maxiter`: the maximum iteration for ARPACK
  - `nev`: the number of eigenvectors estimated by ARPACK
  - `dense`: the threshold for a dense (LAPACK) computation
    instead of a sparse (ARPACK) computation. If the size of
    the matrix is less than `dense` then maxiter and nev are
    ignored and LAPACK computes all the eigenvalues.

Returns
-------
- `v`: The Fiedler vector as an array
- `lam2`: the eigenvalue itself

Example
-------
~~~
# construct an n-node path graph
n = 25
A = sparse(1:n-1,2:n,1.,n,n)
A = A + A'
(v,lam2) = fiedler_vector(A) # returns a cosine-like fiedler vector
# using UnicodePlots; lineplot(v)
~~~
"""
function fiedler_vector(A::SparseMatrixCSC{V,Int};
    tol=1e-12,maxiter=300,dense=96,nev=2,checksym=true) where V
    n = checksquare(A)
    if checksym
        if !issymmetric(A)
            throw(ArgumentError("The input matrix must be symmetric."))
        end
    end

    d = vec(sum(A,dims=1))
    d = sqrt.(d)

    if n == 1
        X = zeros(V,1,2)
        lam2 = 0.
    elseif n <= dense
        ai,aj,av = findnz(A)
        L = sparse(ai,aj,-av./((d[ai].*d[aj])),n,n) # applied sqrt above
        L = Matrix(L) + 2I
        F = eigen!(Symmetric(L))
        lam2 = F.values[2]-1.
        X = F.vectors
    else
        ai,aj,av = findnz(A)
        L = sparse(ai,aj,-av./((d[ai].*d[aj])),n,n) # applied sqrt above
        L = L + sparse(2.0I,n,n)

        (lams,X,nconv) = _symeigs_smallest_arpack(L,nev,tol,maxiter,d)
        lam2 = lams[2]-1.
    end
    x1err = norm(X[:,1]*sign(X[1,1]) - d/norm(d))
    if x1err >= sqrt(tol)
        s = @sprintf("""
        the null-space vector associated with the normalized Laplacian
        was computed inaccurately (diff=%.3e); the Fiedler vector is
        probably wrong or the graph is disconnected""",x1err)
        @warn s
    end

    x = vec(X[:,2])
    if n > 1
        x = x./d # applied sqrt above
    end

    # flip the sign if the number of pos. entries is less than the num of neg. entries
    nneg = sum(x .< 0.)
    if n-nneg < nneg
      x = -x;
    end

    return (x,lam2)
end

fiedler_vector(A::MatrixNetwork{T};kwargs...) where {T} =
    fiedler_vector(sparse_transpose(A);kwargs...) # sparse_transpose converts directly

"""
RankedArray
-----------

This is a data-type that functions as a union between
a sparse and dense array. It allows us to check if
an element is not in the array (i.e. it's too long)
by adding a getindex function to the array type.

It's called a RankedArray because the idea is that
the elements of the array are the ranks of the
nodes in a sorted vector.

This should not be used externally.

Example
-------
v = rand(10)
r = RankedArray(sortperm(v))
haskey(r,11) # returns false
haskey(r,1) # returns true
haskey(r,0) # returns false
"""
mutable struct RankedArray{S}
    data::S
end

import Base: getindex, haskey

haskey(A::RankedArray{S}, x::Int) where {S} = x >= 1 && x <= length(A.data)
getindex(A::RankedArray{S}, i) where {S} = A.data[i]

"""
SweepcutProfile
---------------

This type is the result of a sweepcut operation. It stores
a number of vectors associated with a sweep cut including
the cut and volume associated with each.

Methods
-------
- `bestset(P::SweepcutProfile)` return the best set identified
  in the sweepcut. This is usually what you want.

Example
-------
See `sweepcut`

See also
--------
* `spectral_cut`

"""
struct SweepcutProfile{V,F}
    p::Vector{Int}
    conductance::Vector{F}
    cut::Vector{V}
    volume::Vector{V}
    total_volume::V
    total_nodes::Int

    function SweepcutProfile{V,F}(p::Vector{Int},nnodes::Int,totalvol::V) where {V,F}
        n = length(p)
        new{V,F}(p,Array{F}(undef,n-1),Array{V}(undef,n-1),Array{V}(undef,n-1),totalvol,nnodes)
    end
end

"""
`sweepcut`
--------

A sweepcut is an operation that takes an order to
the vertices of a graph and evaluates the cut and
volume of every partition induced by a prefix of
that ordering. That is, if the order of vertices
is

    v1, v2, v3, ...

the the sweep cut evaluates

    cut({v1}), vol({v1})
    cut({v1,v2}), vol({v1,v2})
    cut({v1,v2,v3}), vol({v1,v2,v3})
    ...

where cut is the total edge weight that connects
vertices in S to the graph. And vol is the
total edge weight connecting

Functions
---------
- `profile = sweepcut(A,x::Vector{T})`
  A standard sweepcut call that will get a sweepcut
  for a vector where x is sorted to get the order
  of the vertices (decreasing order).
  This is useful if you have a vector that should
  indicate good cuts in the graph, such as
  the Fiedler vector.

- `profile = sweepcut(A,p,r,totalvol,maxvol)`
  the in-depth call that all others are converted into.
  We strongly recommend against calling this yourself
  unless you understand the sweepcut code.

Inputs
------
- `A`: the sparse matrix representing the symmetric graph
- `x`: A vector scoring each vertex. This will be sorted and
       turned into one of the other inputs.
- `p`: a partial permutation vector over the vertices of the graph
        This vector needs to list every vertex at most once.
        It could be shorter and need not list every vertex.
- `r`: A general data structure that gives
        the rank of an item in the permutation vector
        p should be sorted in decreasing order so that
        i < j means that r[p[i]] < r[p[j]]
- `totalvol`: the total volume of the graph
- `maxvol`: the maximum volume of any set considered

Returns
-------
A `SweepcutProfile` type with all the information computed
in the sweepcut. Most likely, you want the result `bestset`
as indicated below.

Example
-------
~~~~
A = load_matrix_network("minnesota")
v = fiedler_vector(A)[1] # get the
p = sweepcut(A,v)
S = bestset(p) # get the bestset from the profile
T = spectral_cut(A).set # should give you the same set
# using UnicodePlots; lineplot(p.conductance) # show the conductance
~~~~
"""
function sweepcut(A::SparseMatrixCSC{V,Int}, p::Vector{Int}, r, totalvol::V, maxvol::T) where {V,T}

    n = checksquare(A)
    nlast = length(p)

    if n < nlast
        throw(DimensionMismatch(
            "the permutation vector is too long and should have fewer entries than the graph"))
    end

    F = typeof(one(V)/one(V)) # find the conductance type
    output = SweepcutProfile{V,F}(p,n,totalvol)

    cut = zero(V)
    vol = zero(V)
    colptr = A.colptr
    rowval = rowvals(A)
    nzval = A.nzval

    for (i,v) in enumerate(p)
        deltain = zero(V) # V might be pos. only... so delay subtraction
        deg = zero(V)
        rankv = getindex(r,v)
        for nzi in colptr[v]:(colptr[v+1] - 1)
            nbr = rowval[nzi]
            deg += nzval[nzi]
            if haskey(r,nbr) # our neighbor is ranked
                if getindex(r,nbr) <= rankv # nbr is ranked lower, decrease cut
                    deltain += v == nbr ? nzval[nzi] : 2*nzval[nzi]
                end
            end
        end
        cut += deg
        cut -= deltain
        vol += deg

        # don't assign final values because they are unhelpful
        if i==nlast
            break
        end

        cond = cut/min(vol,totalvol-vol)
        output.conductance[i] = cond
        output.cut[i] = cut
        output.volume[i] = vol
    end

    if nlast == size(A,1)
        @assert abs(cut) <= 1e-12*totalvol
    end
    return output
end

function sweepcut(A::SparseMatrixCSC{V,Int}, x::Vector{T}, vol::V) where {V,T}
    p = sortperm(x,rev=true)
    ranks = Array{Int}(undef, length(x))
    for (i,v) in enumerate(p)
        ranks[v] = i
    end
    r = RankedArray(ranks)
    return sweepcut(A, p, r, vol, Inf)
end

sweepcut(A::SparseMatrixCSC{V,Int}, x::Vector{T}) where {V,T} =
    sweepcut(A, x, sum(A))

function bestset(prof::SweepcutProfile{V,F}) where {V,F}
    nnodes = 0
    if !isempty(prof.conductance)
        bsetind = argmin(prof.conductance)
        bsetvol = prof.volume[bsetind]
        nnodes = bsetind
        if bsetvol > prof.total_volume - bsetvol
            # we want the complement set
            nnodes = prof.total_nodes - bsetind
        end
    end

    bset = Vector{Int}(undef,nnodes)
    if isempty(prof.conductance)
    elseif bsetvol > prof.total_volume - bsetvol
        # ick, we need the complement
        bset[:] = collect(setdiff(BitSet(Int(1):Int(prof.total_nodes)),prof.p[1:bsetind]))
    else
        # easy
        bset[:] = prof.p[1:bsetind]
    end
    return bset
end

"""
`SpectralCut`
-------------
The return type from the `spectral_cut`

Fields
------
- `set`: the small side of the spectral cut
- `A`: the network of the largest strong component of the network
- `lam2`: the eigenvalue of the normalized Laplacian
- `x`: the Fiedler vector for spectral partitioning
- `sweepcut_profile`: the sweepcut profile output
- `comps`: the output from the strong_components function
    check comps.number to get the number of components
- `largest_component`: the index of the largest strong component

The most useful outputs are `set` and `lam2`; the others are provided
for experts who wish to use some of the diagonstics provided.
"""
struct SpectralCut{V,F}
    set::Vector{Int}
    A::SparseMatrixCSC{V,Int}
    lam2::F
    x::Vector{Float64}
    sweepcut_profile::SweepcutProfile{V,F}
    comps::MatrixNetworks.Strong_components_output
    largest_component::Int
end

function show(io::IO,obj::SpectralCut{V,F}) where {V,F}
    println(io,"Spectral cut on adjacency matrix with $(size(obj.A,1)) nodes and $(nnz(obj.A)) non-zeros")
    if obj.comps.number > 1
        println(io, "    formed from the largest connected component of the input matrix")
    end
    println(io,"  conductance = $(minimum(obj.sweepcut_profile.conductance))")
    println(io,"  size = $(length(obj.set))")
end

"""
`spectral_cut`
--------------
`spectral_cut` will produce a spectral cut partition of the graph into
two pieces.

Special cases
* if your graph is disconnected, then we will partition it into the
largest connected component (chosen arbitrary if there are multiple)
and produce a cut of just the largest connected component.
* if your graph is a single node, we will partition it into the empty
cut.

Output
------
The output has type SpectralCut
We always return the smaller side of the cut in terms of total volume.

Inputs
------
- `A`: The sparse matrix or martrix network that you want to
  find the spectral cut for.
- `checksym`: Should we check symmetry?
  Don't set this to false unless you have checked symmetry
  in another call. *This may go away in future interfaces*
- `ccwarn`: Turn off the warning for disconnected graphs.
  This useful in larger subroutines where this is handled
  through another mechanism.
"""
function spectral_cut(A::SparseMatrixCSC{V,Int},checksym::Bool,ccwarn::Bool) where V
    n = checksquare(A)
    if checksym
        if !issymmetric(A)
            throw(ArgumentError("The input matrix must be symmetric."))
        end
    end

    # test for non-negativity
    if !all(nonzeros(A) .>= 0)
        throw(ArgumentError("The input matrix must be non-negative"))
    end

    # need to test for components
    G = MatrixNetwork(n,A.colptr,A.rowval,A.nzval)
    cc = scomponents(G)
    lccind = argmax(cc.sizes)
    B = A
    if cc.number > 1
        if ccwarn
            @warn "The input matrix had $(cc.number) components, running on largest..."
        end
        f = cc.map .== lccind
        B = A[f,f] # using logical indexing is faster than find for large sets (2015/11/14)
        # n=10^5;A=sprand(n,n,25/n);f=rand(n).>0.01;s=find(f);@time sum(A);@time A[f,f]; @time A[s,s];
        subset = findall(f)
    end

    # run the partition
    x,lam2 = fiedler_vector(B;checksym=false)

    totalvol = sum(B)
    sweep = sweepcut(B,x,totalvol)
    bset = bestset(sweep)

    if cc.number > 1
        # map the set back to the original index set
        bset = subset[bset]
    end

    return SpectralCut(bset,B,lam2,x,sweep,cc,lccind)
end

function spectral_cut(A::SparseMatrixCSC{V,Int}) where V
    return spectral_cut(A,true,true)
end

function spectral_cut(A::MatrixNetwork{V}) where V
    return spectral_cut(sparse_transpose(A),true,true)
end


text:
using PGFPlots
using Reel
using Colors

# Setup plotting for animinations
# Set default axis and label styles
Reel.extension(m::MIME"image/svg+xml") = "svg"
Reel.set_output_type("gif") # may be necessary for use in IJulia

function plotReachable(sets,extent,time;zmax=360,timeOffset=0,cmCustom=true,desiredKey=nothing,psi=nothing,nbin=150,width="18cm",height="18cm",title=-1,xlabel="Downrange (kft)",ylabel="Crossrange (kft)",filename="Pics/HCAS_ReachablePlot.png",savePlot=false,saveName=saveName="Pics/HCAS_ReachablePlot.tex",isRow=false,isLast=false,keepReg=nothing)
    #= Plot reachable set
    Inputs:
        sets (Dictionary): Reachable sets at different times (taus)
        extent (list of float): Boundaries of plot [xmin, xmax, ymin, ymax]
        time (int): Time (tau) value to plot (Should be a key in sets)
    Optional Inputs:
        zmax (float): Maximum z-value for colorbar
        timeOffset (int): Offset for time index (tau)
        cmCustom (bool): True if using custom Jet colormap
        desiredKey (string): Previous advisory sequence to plot (should be a key of sets[time]).
            If nothing, the reachable sets of all previous advisory sequences are combined
        psi (float): Show reachable cells that contain a given psi value.
            If nothing, all reachable cells are used
        nbin (int): Number of bins for plotting resolution
        width (string): Width of plot
        height (string): height of plot
        title (string): Title of plot. If title==-1, use a default title
        xlabel (string): X-axis label
        ylabel (string): Y-axis label
        filename (string): Filename for saving png file (if saving)
        savePlot (bool): True if saving the plot as .tex and .png files
        saveName (string): Name for the .tex file
        isRow (bool): True if this plot is a part of a larger row of plots
        isLast (bool): True if this is the last plot in an animation
        keepReg (BitArray): Region to always plot as reachable
    Outputs:
        g (PGFPlots.Axis): Plot of reachable regions
    =#
    set = BitArray(undef,NUMREGIONS)
    if desiredKey != nothing
        set .= sets[time-timeOffset][desiredKey]
    else
        set .= false
        for (k,v) = sets[time-timeOffset]
            set .|= v
        end
    end

    if !cmCustom
        zmax=1
    end

    function getXY(x,y)
        x*=1000
        y*=1000

        if y>YS[end-1]
            y=YS[end-1]
        elseif y<YS[1]
            y=YS[1]
        end
        if x > XS[end-1]
            x=XS[end-1]
        elseif x<XS[1]
            x=XS[1]
        end

        xInd = findall(XS.>x)[1]-1
        yInd = findall(YS.>y)[1]-1
        sumFound=0
        if psi !=nothing
            if psi>PSIS[end-1]
                psi=PSIS[end-1]
            elseif psi<PSIS[1]
                psi = PSIS[1]
            end
            pInd = findall(PSIS.>psi)[1]-1
            ind = indicesToIndex(pInd,yInd,xInd)
            if (keepReg != nothing) && (!keepReg[ind])
                return zmax
            end
            sumFound = set[ind]
        else
            inds = indicesToIndex_Vector(1:NUMP,yInd*ones(Int32,NUMP),xInd*ones(Int32,NUMP))
            if (keepReg != nothing) && (!keepReg[inds[1]])
                return zmax
            end
            sumFound = sum([set[i] for i in inds])
        end
        return sumFound
    end

    if !savePlot && !isRow
        filename=nothing
    end

    if !isRow
        PGFPlots.resetPGFPlotsPreamble()
        PGFPlots.pushPGFPlotsPreamble("\\pgfplotsset{every axis/.append style={title style={font=\\huge},label style={font=\\huge},tick label style={font=\\Large}}}")
    end

    xmin, xmax, ymin, ymax = extent
    if title==-1
        title=@sprintf("Time to Coaltitude: %d s",max(0,time))
        if time<=0
            title=@sprintf("Time at Coaltitude: %d s",-time)
        end
        if isLast
            title="Converged to Steady State"
        end
    end
    style=""
    if ylabel==""
        style = "yticklabels={,,}, scaled y ticks=false"
    end
    cm = ColorMaps.Named("Blues")
    if cmCustom
        cmColors = copy(ColorMaps.Named("Jet").colors)
        cmColors[1] = RGB(1.0,1.0,1.0)
        cm=ColorMaps.RGBArrayMap(cmColors)
    end
    g = Axis([
            Plots.Image(getXY,(xmin,xmax),(ymin,ymax),zmin=0,zmax=zmax,
                xbins=nbin,ybins=nbin,colormap=cm,colorbar=false,filename=filename),
            Plots.Circle(0,0,0.5,style="black,fill=red"),
        ],
        width=width,height=height,xlabel=xlabel,ylabel=ylabel,title=title,style=style)
    if savePlot
        save(saveName,g,include_preamble=false)
    else
        return g
    end
end

function plotReachableRow(sets,extent,times;nbin=200,zmax=360,timeOffset=0,cmCustom=true,savePlot=false,saveName="Pics/HCAS_ReachableRow.tex",keepReg=nothing)
    #= Plot a row of reachable plots, which are the reachable regions at different times
    Inputs:
        sets (Dictionary): Reachable sets at different times (taus)
        extent (list of float): Boundaries of plot [xmin, xmax, ymin, ymax]
        times (array of int): Times (taus) at which we want to plot (Should be keys in sets)
    Optional Inputs:
        zmax (float): Maximum z-value for colorbar
        timeOffset (int): Offset for time index (tau)
        cmCustom (bool): True if using custom Jet colormap
        nbin (int): Number of bins for plotting resolutiones
        saveName (string): Name for the .tex file
        keepReg (BitArray): Region to always plot as reachable
    Outputs:
        g (PGFPlots.GroupPlot): Row of reachable region plots
    =#
    PGFPlots.resetPGFPlotsPreamble()
    g = GroupPlot(length(times),1, groupStyle = "horizontal sep=0.35cm", style="height=4.4cm, width=4.4cm")
    width=nothing; height=nothing;ylabel="Crossrange (kft)";xlabel="Downrange (kft)"
    for time in times
        fn = nothing
        if savePlot
            fn = @sprintf("Pics/HCAS_ReachRow_%d.png",time)
        end
        push!(g,plotReachable(sets,extent,time,nbin=nbin,width=width,height=height,title=-1,xlabel=xlabel,ylabel=ylabel,filename=fn,zmax=zmax,timeOffset=timeOffset,cmCustom=cmCustom,isRow=true,keepReg=keepReg))
        ylabel=""
    end
    if savePlot
        save(saveName,g,include_preamble=false)
        return
    else
        return g
    end
end

function animateReach(sets,extent,timeMax,timeMin,save_name;timeOffset=0,cmCustom=true,zmax=360,nbin=150,keepReg=nothing)
    #= Create an animation of reachable set plots. Begins with timeMax and counts down to timeMin
    Inputs:
        sets (Dictionary): Reachable sets at different times (taus)
        extent (list of float): Boundaries of plot [xmin, xmax, ymin, ymax]
        timeMin (int): Minimum time (tau) value to plot (Should be a key in  sets)
        timeMax (int): Maximum time (tau) value to plot (Should be a key in  sets)
        save_name (string): Name of gif file (without the .gif)
    Optional Inputs:
        zmax (float): Maximum z-value for colorbar
        timeOffset (int): Offset for time index (tau)
        cmCustom (bool): True if using custom Jet colormap
        nbin (int): Number of bins for plotting resolution
        keepReg (BitArray): Region to always plot as reachable
    Outputs:
        frames (Reel.Frames): Animation of reachable regions changing over time
    =#

    currentFrac = 0
    deltaFrac = 0.1
    frames = Frames(MIME("image/svg+xml"), fps=10)
    for frame in timeMax:-1:timeMin
        if (timeMax-frame)/(timeMax-timeMin)>=currentFrac
            println("Completed "*string(round(currentFrac*100))*"%")
            currentFrac+=deltaFrac
        end
        push!(frames,plotReachable(sets,extent,frame,timeOffset=timeOffset,cmCustom=cmCustom,zmax=zmax,nbin=nbin,isLast=frame==timeMin-1,keepReg=keepReg))
    end
    write(save_name*".gif",frames)
    frames
end

max_stars_repo_name: cherrywoods/HorizontalCAS

max_issues_repo_name: ahmed-irfan/HorizontalCAS

max_forks_repo_name: ahmed-irfan/HorizontalCAS


text:
@testset "Testing RedBlack Tree" begin
    @test reverse(ReuseDistance.Left) == ReuseDistance.Right
    @test reverse(ReuseDistance.Right) == ReuseDistance.Left

    @testset "Testing Accessors" begin
        # Initialize the tree with length zero.
        # This will test the corner case handling of setting the minimum node vector
        # length to 1 AND test the length extension logic.
        tree = ReuseDistance.RedBlackTree{Int}(0)

        # Get a new index from the tree
        n = tree[]
        @test n == 1
        tree[n] = 500
        @test tree[n] == 500

        # Null all other pointer fields
        ReuseDistance.setcolor!(tree, n, ReuseDistance.Black)
        ReuseDistance.setparent!(tree, n, 0)
        ReuseDistance.setleft!(tree, n, 0)
        ReuseDistance.setright!(tree, n, 0)

        @test ReuseDistance.getcolor(tree, n) == ReuseDistance.Black
        @test ReuseDistance.getparent(tree, n) == 0
        @test ReuseDistance.leftchild(tree, n) == 0
        @test ReuseDistance.rightchild(tree, n) == 0

        # Setting the color to `Red` should not affect the parent
        ReuseDistance.setcolor!(tree, n, ReuseDistance.Red)
        @test ReuseDistance.getcolor(tree, n) == ReuseDistance.Red
        @test ReuseDistance.getparent(tree, n) == 0
        @test ReuseDistance.leftchild(tree, n) == 0
        @test ReuseDistance.rightchild(tree, n) == 0

        # Similarly, setting the parent should not affect the color.
        ReuseDistance.setparent!(tree, n, typemax(Int64))
        @test ReuseDistance.getcolor(tree, n) == ReuseDistance.Red
        @test ReuseDistance.getparent(tree, n) == typemax(Int64)
        @test ReuseDistance.leftchild(tree, n) == 0
        @test ReuseDistance.rightchild(tree, n) == 0

        ReuseDistance.setleft!(tree, n, 10)
        ReuseDistance.setright!(tree, n, 33)
        @test ReuseDistance.getcolor(tree, n) == ReuseDistance.Red
        @test ReuseDistance.getparent(tree, n) == typemax(Int64)
        @test ReuseDistance.leftchild(tree, n) == 10
        @test ReuseDistance.rightchild(tree, n) == 33

        # Test misc accessors
        @test ReuseDistance.child_direction(tree, n, 10) == ReuseDistance.Left
        @test ReuseDistance.child_direction(tree, n, 33) == ReuseDistance.Right
        @test ReuseDistance.getparent(tree, n) == typemax(Int64)
        m = tree[]
        @test m == 2
        ReuseDistance.setparent!(tree, n, 0)
        ReuseDistance.setparent!(tree, m, n)
        ReuseDistance.setleft!(tree, n, m)
        @test ReuseDistance.child_direction(tree, m) == ReuseDistance.Left

        ReuseDistance.setleft!(tree, n, 0)
        ReuseDistance.setright!(tree, n, m)
        @test ReuseDistance.child_direction(tree, m) == ReuseDistance.Right
        @test ReuseDistance.getparent(tree, m) == n
    end

    @testset "Testing Tree Operations" begin
        function clear!(tree, n)
            ReuseDistance.setcolor!(tree, n, ReuseDistance.Black)
            ReuseDistance.setparent!(tree, n, 0)
            ReuseDistance.setleft!(tree, n, 0)
            ReuseDistance.setright!(tree, n, 0)
            tree[n] = 0
            return n
        end

        # Now, we get a bit more elaborate with our tree setup, testing out tree rotations
        # and the accessors for uncles, nephews etc.
        #
        # Out tree will look like this.
        #
        #         A
        #        / \
        #       B   C
        #      / \   \
        #     D   E   F
        function setup_tree()
            tree = ReuseDistance.RedBlackTree{Int}(10)
            a = clear!(tree, tree[])
            b = clear!(tree, tree[])
            c = clear!(tree, tree[])
            d = clear!(tree, tree[])
            e = clear!(tree, tree[])
            f = clear!(tree, tree[])

            tree.root = 1
            ReuseDistance.setleft!(tree, a, b)
            ReuseDistance.setparent!(tree, b, a)

            ReuseDistance.setright!(tree, a, c)
            ReuseDistance.setparent!(tree, c, a)

            ReuseDistance.setright!(tree, c, f)
            ReuseDistance.setparent!(tree, f, c)

            ReuseDistance.setleft!(tree, b, d)
            ReuseDistance.setparent!(tree, d, b)
            ReuseDistance.setright!(tree, b, e)
            ReuseDistance.setparent!(tree, e, b)
            return tree, (; a, b, c, d, e, f)
        end

        #####
        ##### Tree Rotations
        #####

        # After a right rotation on A
        #
        #         B
        #        / \
        #       D   A
        #          / \
        #         E   C
        #              \
        #               F

        tree, nodes = setup_tree()
        ReuseDistance.rotateright!(tree, nodes.a)
        @test tree.root == nodes.b
        @test ReuseDistance.getparent(tree, nodes.b) == 0

        @test ReuseDistance.leftchild(tree, nodes.b) == nodes.d
        @test ReuseDistance.getparent(tree, nodes.d) == nodes.b
        @test ReuseDistance.rightchild(tree, nodes.b) == nodes.a
        @test ReuseDistance.getparent(tree, nodes.a) == nodes.b

        @test ReuseDistance.leftchild(tree, nodes.a) == nodes.e
        @test ReuseDistance.getparent(tree, nodes.e) == nodes.a

        @test ReuseDistance.rightchild(tree, nodes.a) == nodes.c
        @test ReuseDistance.getparent(tree, nodes.c) == nodes.a

        # After a left rotation on A
        #           C
        #          / \
        #         A   F
        #        /
        #       B
        #      / \
        #     D   E
        tree, nodes = setup_tree()
        ReuseDistance.rotateleft!(tree, nodes.a)
        @test tree.root == nodes.c
        @test ReuseDistance.getparent(tree, nodes.c) == 0
        @test ReuseDistance.rightchild(tree, nodes.c) == nodes.f
        @test ReuseDistance.getparent(tree, nodes.f) == nodes.c

        @test ReuseDistance.leftchild(tree, nodes.c) == nodes.a
        @test ReuseDistance.getparent(tree, nodes.a) == nodes.c

        @test ReuseDistance.rightchild(tree, nodes.a) == 0
        @test ReuseDistance.leftchild(tree, nodes.a) == nodes.b
        @test ReuseDistance.getparent(tree, nodes.b) == nodes.a

        # Finally, perform a rotation that doesn't change the root
        # Right rotation on `B`.
        #
        #
        #         A
        #        / \
        #       D   C
        #        \   \
        #         B   F
        #          \
        #           E
        tree, nodes = setup_tree()
        ReuseDistance.rotateright!(tree, nodes.b)
        @test tree.root == nodes.a
        @test ReuseDistance.leftchild(tree, nodes.a) == nodes.d
        @test ReuseDistance.getparent(tree, nodes.d) == nodes.a
        @test ReuseDistance.leftchild(tree, nodes.d) == 0
        @test ReuseDistance.rightchild(tree, nodes.d) == nodes.b
        @test ReuseDistance.getparent(tree, nodes.b) == nodes.d

        #####
        ##### Search
        #####

        tree, nodes = setup_tree()
        tree[nodes.a] = 20
        tree[nodes.b] = 10
        tree[nodes.c] = 30
        tree[nodes.d] = 5
        tree[nodes.e] = 15
        tree[nodes.f] = 35

        @test ReuseDistance.search(tree, 20) == (nodes.a, ReuseDistance.Found)
        @test ReuseDistance.search(tree, 10) == (nodes.b, ReuseDistance.Found)
        @test ReuseDistance.search(tree, 30) == (nodes.c, ReuseDistance.Found)
        @test ReuseDistance.search(tree, 5) == (nodes.d, ReuseDistance.Found)
        @test ReuseDistance.search(tree, 15) == (nodes.e, ReuseDistance.Found)
        @test ReuseDistance.search(tree, 35) == (nodes.f, ReuseDistance.Found)

        # Look for findings
        @test ReuseDistance.search(tree, 25) == (nodes.c, ReuseDistance.NotFoundLeft)
        @test ReuseDistance.search(tree, 16) == (nodes.e, ReuseDistance.NotFoundRight)
    end

    #####
    ##### Insertion!
    #####

    @testset "Testing Insertion" begin
        tree = ReuseDistance.RedBlackTree{Int}()
        for i in 1:100
            push!(tree, i)
        end

        @test ReuseDistance.validate(tree)
        @test length(tree) == 100
        @test all(in(tree), 1:100)
        @test !any(in(tree), 101:200)
    end

    #####
    ##### Stress Testing
    #####

    @testset "Performing Stress Test" begin
        groundtruth = Set{Int}()
        tree = ReuseDistance.RedBlackTree{Int}()
        domain = 1:10000
        numswaps = 1000
        increment = 0.02

        swapcount = 0
        add = true
        meter = ProgressMeter.Progress(numswaps, 1)
        while swapcount <= numswaps
            swapchance = 0.0
            while rand() > swapchance
                numops = rand(1:10)
                for _ in Base.OneTo(numops)
                    if add
                        i = rand(domain)
                        while length(groundtruth) < length(domain) && in(i, groundtruth)
                            i = rand(domain)
                        end
                        push!(tree, i)
                        push!(groundtruth, i)
                    elseif !isempty(groundtruth)
                        i = rand(groundtruth)
                        delete!(tree, i)
                        delete!(groundtruth, i)
                    end
                end
                swapchance += increment
            end

            passed = true
            for i in domain
                if in(i, groundtruth)
                    passed &= in(i, tree)
                else
                    passed &= !in(i, tree)
                end
            end

            @test passed
            @test ReuseDistance.validate(tree)
            swapcount += 1
            add = !add
            ProgressMeter.next!(meter)
        end
    end
end


text:
const _pkg_root = dirname(dirname(@__FILE__))
const _pkg_deps = joinpath(_pkg_root,"deps")
const _pkg_assets = joinpath(_pkg_root,"assets")

!isdir(_pkg_assets) && mkdir(_pkg_assets)

deps = [
    "https://raw.githubusercontent.com/piever/BulmaThemes/master/public/styles/main.css",
    "https://cdnjs.cloudflare.com/ajax/libs/bulma/0.7.1/css/bulma.min.css",
    "https://cdn.jsdelivr.net/npm/bulma-extensions@1.0.14/bulma-slider/dist/bulma-slider.min.css",
    "https://cdn.jsdelivr.net/npm/bulma-extensions@1.0.14/bulma-switch/dist/bulma-switch.min.css",
    "https://cdn.jsdelivr.net/npm/bulma-extensions@1.0.14/bulma-checkradio/dist/bulma-checkradio.min.css",
    "https://cdn.jsdelivr.net/npm/bulma-extensions@1.0.14/bulma-accordion/dist/bulma-accordion.min.css",
    "https://cdn.jsdelivr.net/npm/bulma-extensions@1.0.14/bulma-tooltip/dist/bulma-tooltip.min.css",
]
for dep in deps
    download(dep, joinpath(_pkg_assets, splitdir(dep)[2]))
end


text:
"""
	$(SIGNATURES)

Histogram.
"""
function histogram_plot(dataV; fig = blank_plot(), pos = (1,1), kwargs...)
    args = merge(histogram_defaults(), kwargs);
    ax = fig[pos...] = Axis(fig; args...);
    hist!(ax, dataV; args...);
    return fig, ax
end

histogram_defaults() = Dict([
    :strokewidth => 0,
    :color => main_color()
]);

"""
	$(SIGNATURES)

Contour plot.
"""
function contour_plot(x, y, z; fig = blank_plot(), pos = (1,1), kwargs...)
    args = merge(contour_defaults(), kwargs);
    ax = make_axis(fig, pos; args...);
    contour!(ax, x, y, z; args...);
    return fig, ax
end

contour_defaults() = Dict([
    :linewidth => 3,
]);

# ----------------


text:
using Test
include("./../src/Lar.jl")

include("./arrangement/planar_arrangement.jl")
include("./model.jl")


text:
using LazyCSV
using Test

DO_BENCHMARK = false
DEBUG_ISSUE = false

# creates a IO instance using the given CSV string
function csv_io(csv::AbstractString)
    csv_io = IOBuffer()
    print(csv_io, csv)
    seekstart(csv_io)
    csv_io
end
csv_io(csv::IO) = csv

function csv_count_lines(csv; header_exists::Bool=false)
    csv_file = LazyCSV.csvread(csv; header_exists=header_exists, eager_parse_fields=false)
    count_lines(csv_file)
end

function csv_count_fields(csv; delim=',', header_exists::Bool=false, quotechar=LazyCSV.DEFAULT_QUOTE, escapechar=quotechar)
    csv_file = LazyCSV.csvread(csv; delim=delim, header_exists=header_exists,
                               eager_parse_fields=true, quotechar=quotechar, escapechar=escapechar)
    count_fields(csv_file)
end

function csv_string(csv, num_fields::Int = -1;
                    delim=',', header_exists::Bool=false, quotechar=LazyCSV.DEFAULT_QUOTE, escapechar=quotechar)
    csv_file = LazyCSV.csvread(csv; delim=delim, header_exists=header_exists,
                               eager_parse_fields=true, quotechar=quotechar, escapechar=escapechar)
    LazyCSV.csv_string(csv_file, num_fields)
end

function typed_csv_string(csv, field_types::Type{T}, fallback_type::LazyCSV.FieldType=LazyCSV.NO_TYPE,
                          other_type_handlers=Tuple{}();
                          delim=',', header_exists::Bool=false,
                          quotechar=LazyCSV.DEFAULT_QUOTE, escapechar=quotechar) where {T}
    csv_file = LazyCSV.csvread(csv; delim=delim, header_exists=header_exists,
                               eager_parse_fields=true, quotechar=quotechar, escapechar=escapechar)
    LazyCSV.typed_csv_string(csv_file, field_types, fallback_type, other_type_handlers)
end

function csv_equals(base_csv, to_csv; delim=',')
    base_csv_io = csv_io(base_csv)
    to_csv_io = csv_io(to_csv)
    for (base_line, to_line) in zip(eachline(base_csv_io), eachline(to_csv_io))
        @test csv_line_equals(base_line, to_line; delim=delim) ||
            error("$(strip_csv_line(base_line, delim)) != $(strip_csv_line(to_line, delim))")
    end
end

isspecial_in_regex(delim) = delim == '|'

escape_char(delim) = if isspecial_in_regex(delim)
        "\\$delim"
    else
        "$delim"
    end

function strip_csv_line(line, delim)
    escaped_delim = escape_char(delim)

    line = replace(line, r"^\s*\"" => "")
    line = replace(line, r"^\s*" => "")
    line = replace(line, r"\"\s*$" => "")
    line = replace(line, r"\s*$" => "")
    line = replace(line, Regex("$(escaped_delim)\\s*\\\"") => delim)
    line = replace(line, Regex("\\\"\\s*$(escaped_delim)") => delim)
    line = replace(line, Regex("\\s*$(escaped_delim)\\s*") => delim)
    line = replace(line, "\"\"" => "\"")
    line
end

function csv_line_equals(base_line, to_line; delim=',')
    if base_line != to_line
        strip_csv_line(base_line, delim) == strip_csv_line(to_line, delim)
    else
        true
    end
end

function simple_csv_test(csv_str, num_lines, num_fields; delim=',', quotechar='"',
                         escapechar=quotechar, header_exists::Bool=false)
    computed_num_lines = csv_count_lines(csv_io(csv_str))
    @test computed_num_lines == num_lines || error("$computed_num_lines != $num_lines in \n----------------\n$csv_str\n----------------")
    computed_num_fields = csv_count_fields(csv_io(csv_str); delim=delim, quotechar=quotechar, escapechar=escapechar)
    @test computed_num_fields == num_fields || error("$computed_num_fields != $num_fields in \n----------------\n$csv_str\n----------------")
    generated_csv = csv_string(csv_io(csv_str);
       delim=delim, quotechar=quotechar, escapechar=escapechar,
       header_exists=header_exists)
    csv_equals(replace(csv_str, "\r" => "\r\n"), generated_csv; delim=delim)
end

count_errors(output) = count_rec_errors(output) + count_field_errors(output)
count_rec_errors(output) = count(x -> true, eachmatch(r"ERROR", output))
count_field_errors(output) = count(x -> true, eachmatch(r"FIELD_ERR", output))

@testset "LazyCSV tests" begin
@testset "Untyped Tests" begin
    lineitem_sample = """
    1|155190|7706|1|17|21168.23|0.04|0.02|N|O|1996-03-13|1996-02-12|1996-03-22|DELIVER IN PERSON|TRUCK|egular courts above the
    1|67310|7311|2|36|45983.16|0.09|0.06|N|O|1996-04-12|1996-02-28|1996-04-20|TAKE BACK RETURN|MAIL|ly final dependencies: slyly bold
    1|63700|3701|3|8|13309.60|0.10|0.02|N|O|1996-01-29|1996-03-05|1996-01-31|TAKE BACK RETURN|REG AIR|riously. regular, express dep
    1|2132|4633|4|28|28955.64|0.09|0.06|N|O|1996-04-21|1996-03-30|1996-05-16|NONE|AIR|lites. fluffily even de
    1|24027|1534|5|24|22824.48|0.10|0.04|N|O|1996-03-30|1996-03-14|1996-04-01|NONE|FOB| pending foxes. slyly re
    1|15635|638|6|32|49620.16|0.07|0.02|N|O|1996-01-30|1996-02-07|1996-02-03|DELIVER IN PERSON|MAIL|arefully slyly ex
    2|106170|1191|1|38|44694.46|0.00|0.05|N|O|1997-01-28|1997-01-14|1997-02-02|TAKE BACK RETURN|RAIL|ven requests. deposits breach a
    3|4297|1798|1|45|54058.05|0.06|0.00|R|F|1994-02-02|1994-01-04|1994-02-23|NONE|AIR|ongside of the furiously brave acco
    3|19036|6540|2|49|46796.47|0.10|0.00|R|F|1993-11-09|1993-12-20|1993-11-24|TAKE BACK RETURN|RAIL| unusual accounts. eve
    3|128449|3474|3|27|39890.88|0.06|0.07|A|F|1994-01-16|1993-11-22|1994-01-23|DELIVER IN PERSON|SHIP|nal foxes wake.
    """

    simple_csv_test(lineitem_sample, 10, 160; delim='|')

    quoted_csv1 = """
    John,Doe,120 jefferson st.,Riverside, NJ, 08075
    Jack,McGinnis,220 hobo Av.,Phila, PA,09119
    "John ""Da Man""\",Repici,120 Jefferson St.,Riverside, NJ,08075
    Stephen,Tyler,"7452 Terrace ""At the Plaza"" road",SomeTown,SD, 91234
    ,Blankman,,SomeTown, SD, 00298
    ,Blankman,"",SomeTown, SD, 00298
    "Joan ""\""the bone"", Anne",Jet,"9th, at Terrace plc",Desert City,CO,00123
    """

    simple_csv_test(quoted_csv1, 7, 42)

    quoted_csv2 = """
    John,Doe,120 jefferson st.,Riverside, NJ, 08075
    Jack,McGinnis,220 hobo Av.,Phila, PA,09119
    "John %"Da Man%"",Repici,120 Jefferson St.,Riverside, NJ,08075
    Stephen,Tyler,"7452 Terrace %"At the Plaza%" road",SomeTown,SD, 91234
    ,Blankman,,SomeTown, SD, 00298
    ,Blankman,"",SomeTown, SD, 00298
    "Joan %"%"the bone%", Anne",Jet,"9th, at Terrace plc",Desert City,CO,00123
    """

    simple_csv_test(quoted_csv2, 7, 42; escapechar='%')

    airtravel_csv = """
    "Month", "1958", "1959", "1960"
    "JAN",  340,  360,  417
    "FEB",  318,  342,  391
    "MAR",  362,  406,  419
    "APR",  348,  396,  461
    "MAY",  363,  420,  472
    "JUN",  435,  472,  535
    "JUL",  491,  548,  622
    "AUG",  505,  559,  606
    "SEP",  404,  463,  508
    "OCT",  359,  407,  461
    "NOV",  310,  362,  390
    "DEC",  337,  405,  432
    """

    simple_csv_test(airtravel_csv, 13, 52; header_exists=true)

    biostats_csv = """
    "Name",     "Sex", "Age", "Height (in)", "Weight (lbs)"
    "Alex",       "M",   41,       74,      170
    "Bert",       "M",   42,       68,      166
    "Carl",       "M",   32,       70,      155
    "Dave",       "M",   39,       72,      167
    "Elly",       "F",   30,       66,      124
    "Fran",       "F",   33,       66,      115
    "Gwen",       "F",   26,       64,      121
    "Hank",       "M",   30,       71,      158
    "Ivan",       "M",   53,       72,      175
    "Jake",       "M",   32,       69,      143
    "Kate",       "F",   47,       69,      139
    "Luke",       "M",   34,       72,      163
    "Myra",       "F",   23,       62,       98
    "Neil",       "M",   36,       75,      160
    "Omar",       "M",   38,       70,      145
    "Page",       "F",   31,       67,      135
    "Quin",       "M",   29,       71,      176
    "Ruth",       "F",   28,       65,      131
    """

    simple_csv_test(biostats_csv, 19, 95; header_exists=true)

    quoted_csv3 = """
    John,Doe,120 jefferson st.,Riverside, NJ, 08075
    Jack,McGinnis,220 hobo Av.,Phila, PA,09119,error1
    "John ""Da Man""\",Repici,120 Jefferson St.,Riverside, NJ,08075
    Stephen,Tyler,"7452 Terrace ""At the Plaza"" road",SomeTown,SD, 91234
    ,Blankman,,SomeTown, SD, 00298
    ,Blankman,"",SomeTown, error2
    "Joan ""\""the bone"", Anne",Jet,"9th, at Terrace plc",Desert City,CO,00123
    """

    quoted_csv3_out = csv_string(csv_io(quoted_csv3), 6)
    @test count_rec_errors(quoted_csv3_out) == 2

    simple_csv_test(airtravel_csv, 13, 52; header_exists=true)

end

@testset "Typed Tests" begin
    biostats_csv1 = """
    "Name",     "Sex", "Age", "Height (in)", "Weight (lbs)"
    "Alex",       "M",   41,       74,      170
    "Bert",       "M",   42,       68,      166
    "Carl",       "M",   32,       70,      155
    """
    biostats_csv1_out = typed_csv_string(csv_io(biostats_csv1), Tuple{String,Char,Int,Int,Int}; header_exists=true)
    @test "\"Name (string)\",\"Sex (char)\",\"Age (int)\",\"Height (in) (int)\",\"Weight (lbs) (int)\"" == first(eachline(csv_io(biostats_csv1_out)))
    @test count_errors(biostats_csv1_out) == 0

    biostats_csv2 = """
    "Name",     "Sex", "Age", "Height (in)", "Weight (lbs)"
    "Alex",       "M",   41,       74,      170
    "Bert",       "M",   42.1,     68,      166
    "Carl",       "M",   32,       70,      155
    """
    biostats_csv2_out = typed_csv_string(csv_io(biostats_csv2), Tuple{String,Char,Int,Int,Int}; header_exists=true)
    @test "\"Name (string)\",\"Sex (char)\",\"Age (int)\",\"Height (in) (int)\",\"Weight (lbs) (int)\"" == first(eachline(csv_io(biostats_csv2_out)))
    @test count_rec_errors(biostats_csv2_out) == 0
    @test count_field_errors(biostats_csv2_out) == 1

    biostats_csv3_out = typed_csv_string(csv_io(biostats_csv2), Tuple{String,Char,Union{Int,Float64},Int,Int}; header_exists=true)
    @test "\"Name (string)\",\"Sex (char)\",\"Age (int,float)\",\"Height (in) (int)\",\"Weight (lbs) (int)\"" == first(eachline(csv_io(biostats_csv3_out)))
    @test count_errors(biostats_csv3_out) == 0

    biostats_csv4 = """
    "Name",     "Sex", "Age", "Height (in)", "Weight (lbs)"
    "Alex",       "M",   41,       74,      170
    "Bert",          ,   42.1,     68,      166
    "Carl",       "M",   32,       70,      155
    """
    biostats_csv4_out = typed_csv_string(csv_io(biostats_csv4), Tuple{String,Char,Union{Int,Float64},Int,Int}; header_exists=true)
    @test "\"Name (string)\",\"Sex (char)\",\"Age (int,float)\",\"Height (in) (int)\",\"Weight (lbs) (int)\"" == first(eachline(csv_io(biostats_csv4_out)))
    @test count_rec_errors(biostats_csv4_out) == 0
    @test count_field_errors(biostats_csv4_out) == 1

    biostats_csv5_out = typed_csv_string(csv_io(biostats_csv4), Tuple{String,Union{Char,Missing},Union{Int,Float64},Int,Int}; header_exists=true)
    @test "\"Name (string)\",\"Sex (char,missing)\",\"Age (int,float)\",\"Height (in) (int)\",\"Weight (lbs) (int)\"" == first(eachline(csv_io(biostats_csv5_out)))
    @test count_errors(biostats_csv5_out) == 0
end

additional_fields = Dict("taxables.csv" => 4, "deniro.csv" => 3, "oscar_age_male.csv" => 2,
                         "oscar_age_female.csv" => 2, "freshman_lbs.csv" => 2)

files_without_header = Dict{String, Int}()

@testset "Untyped File Tests" begin
    csv_dir = abspath(joinpath(dirname(@__FILE__), "sample_csv"))
    for csv_file in readdir(csv_dir)
        endswith(csv_file, ".csv") || continue
        DEBUG_ISSUE && println("$csv_file:")
        csv_file_path = joinpath(csv_dir, csv_file)
        csv_file_io = open(csv_file_path, "r")
        csv_file_str = read(csv_file_io, String)
        line_count = count(x -> x == '\n' || x == '\r', csv_file_str)
        delim_count = count(x -> x == ',', csv_file_str)
        field_count = delim_count + line_count
        simple_csv_test(csv_file_str, line_count, field_count-get(additional_fields,
                        csv_file, 0); header_exists=!haskey(files_without_header, csv_file))
    end
end

end


text:
__precompile__()

module OrdinalMultinomialModels

using LinearAlgebra
using Distributions, Reexport, StatsModels
@reexport using StatsBase
@reexport using GLM
@reexport using Ipopt
@reexport using NLopt
using MathOptInterface
import StatsBase: coef, coeftable, deviance, dof, fit, modelmatrix, nobs,
response, score, stderror, weights
import StatsModels: drop_intercept
import LinearAlgebra: BlasReal

const MOI = MathOptInterface

export
    # types
    AbstractOrdinalMultinomialModel,
    OrdinalMultinomialModel,
    OrdinalMultinomialLrtTest,
    OrdinalMultinomialScoreTest,
    # functions
    coef,
    coeftable,
    confint,
    cor,
    deviance,
    dof,
    dof_residual,
    drop_intercept,
    loglikelihood,
    modelmatrix,
    nobs,
    polr,
    polrtest,
    loglikelihood,
    loglikelihood!,
    fit,
    fitted,
    predict,
    response,
    rpolr,
    score,
    stderror,
    vcov,
    weights,
    set_optimizer_attributes

abstract type AbstractOrdinalMultinomialModel <: RegressionModel end
drop_intercept(::Type{AbstractOrdinalMultinomialModel}) = true

"""
    OrdinalMultinomialModel

Data, parameters, and various derived variables for the proportional odds
logistic regression model.
"""
struct OrdinalMultinomialModel{TY<:Integer, T<:BlasReal, TL<:GLM.Link} <: MOI.AbstractNLPEvaluator
    # dimensions
    "`n`: number of observations"
    n::Int
    "`p`: number of covariates, excluding intercept"
    p::Int
    "`J`: number of categories"
    J::Int
    "`npar`: number of parameters"
    npar::Int
    # data
    "`Y`: response vector"
    Y::Vector{TY}
    "`X`: n-by-p covariate matrix"
    X::Matrix{T}
    "`wts`: prior observation weights, can be empty"
    wts::Vector{T}
    # parameters
    "`θ`: intecept parameters, satisfying `θ[1]≤...≤θ[J-1]`"
    θ::Vector{T}
    "`α`: unconstrained parameterization of `θ`"
    α::Vector{T}
    "`β`: regression coefficients"
    β::Vector{T}
    "`link`: link function"
    link::TL
    # working parameters
    "`η`: linear systematic component `Xβ`"
    η::Vector{T}
    "`dθdα`: Jacobian dθdα[i, j] = dθj/dαi"
    dθdα::Matrix{T}
    "`∇`: gradient wrt (θ, β)"
    ∇::Vector{T}
    "`FIM`: Fisher information matrix at (θ, β)"
    FIM::Matrix{T}
    "`vcov`:  `inv(FIM)`"
    vcov::Matrix{T}
    "`wtwk`: working weights"
    wtwk::Vector{T}
    "`wt∂β``: weight vector for computing derivative, `n x 1`"
    wt∂β::Vector{T}
    "`wt∂θ∂β`: weight matrix for computing Hessian, n x (J - 1)"
    wt∂θ∂β::Matrix{T}
    "`wt∂β∂β`: weight vector for computing Hessian, n x 1"
    wt∂β∂β::Vector{T}
    "`scratchm1`: scratch matrix of same size as `X`"
    scratchm1::Matrix{T}
end

# Constructor
function OrdinalMultinomialModel(
    X::Matrix{T},
    y::Vector{TY},
    wts::Vector{T} = similar(X, 0),
    link::GLM.Link = LogitLink()) where {TY <: Integer, T <: BlasReal}
    # check y has observations in each category
    yct = counts(y)
    J   = length(yct)
    J < 2 && throw(ArgumentError("Response must have 3 or more levels"))
    for j in 1:J
        yct[j] < 1 && throw(ArgumentError("No observations in category $j"))
    end
    if !isempty(wts)
        lw, ly = length(wts), length(y)
        lw ≠ ly && throw(ArgumentError("wts has length $lw, should be 0 or $ly"))
        minw = minimum(wts)
        minw < 0 && throw(ArgumentError("wts should be nonnegative, found entry $minw"))
    end
    n, p   = size(X)
    npar   = J - 1 + p
    η      = zeros(T, n)
    θ      = zeros(T, J - 1)
    α      = zeros(T, J - 1)
    β      = zeros(T, p)
    dθdα   = zeros(T, J - 1, J - 1)
    ∇      = zeros(T, J - 1 + p)
    FIM    = zeros(T, J - 1 + p, J - 1 + p)
    vcov   = zeros(T, J - 1 + p, J - 1 + p)
    wtwk   = zeros(T, n)
    wt∂β   = zeros(T, n)
    wt∂θ∂β = zeros(T, n, J - 1)
    wt∂β∂β = zeros(T, n)
    scratchm1 = zero(X)
    OrdinalMultinomialModel{eltype(y), eltype(X), typeof(link)}(n, p, J, npar, y, X, wts,
        θ, α, β, link, η, dθdα, ∇, FIM, vcov, wtwk, wt∂β, wt∂θ∂β, wt∂β∂β, scratchm1)
end
OrdinalMultinomialModel(X, y, link) = OrdinalMultinomialModel(X, y, similar(X, 0), link)

coef(m::OrdinalMultinomialModel) = [m.θ; m.β]
deviance(m::OrdinalMultinomialModel) = -2loglikelihood!(m, false, false)
dof(m::OrdinalMultinomialModel) = m.npar
dof_residual(m::OrdinalMultinomialModel) = m.n - m.npar
fitted(m::OrdinalMultinomialModel) = nothing # TODO
loglikelihood(m::OrdinalMultinomialModel) = loglikelihood!(m, false, false)
modelmatrix(m::OrdinalMultinomialModel) = m.X
nobs(m::OrdinalMultinomialModel) = m.n
predict(m::OrdinalMultinomialModel) = nothing # TODO
response(m::OrdinalMultinomialModel) = m.Y
score(m::OrdinalMultinomialModel) = m.∇
stderror(m::OrdinalMultinomialModel) = sqrt.(diag(m.vcov))
vcov(m::OrdinalMultinomialModel) = m.vcov
weights(m::OrdinalMultinomialModel) = m.wts

function coeftable(m::OrdinalMultinomialModel)
    cc = coef(m)
    se = stderror(m)
    tt = cc ./ se
    CoefTable(hcat(cc, se, tt, ccdf.(FDist(1, dof_residual(m)), abs2.(tt))),
              ["Estimate", "Std.Error", "t value", "Pr(>|t|)"],
              [["θ$i" for i = 1:m.J-1]; ["β$i" for i = 1:m.p]], 4)
end

confint(m::OrdinalMultinomialModel, level::Real) = hcat(coef(m), coef(m)) +
    stderror(m) * quantile(Normal(), (1. - level) / 2.) * [1. -1.]
confint(m::OrdinalMultinomialModel) = confint(m, 0.95)

function cor(m::OrdinalMultinomialModel)
    Σ = copy(vcov(m))
    invstd = similar(Σ, size(Σ, 1))
    for i in eachindex(invstd)
        invstd[i] = 1 / sqrt(Σ[i, i])
    end
    lmul!(Diagonal(invstd), rmul!(Σ, Diagonal(invstd)))
end

function set_optimizer_attributes(m::MathOptInterface.AbstractOptimizer, pairs::Pair...)
    for (name, value) in pairs
        MathOptInterface.set(
            m,
            MathOptInterface.RawOptimizerAttribute(name),
            value
        )
    end
end

include("ordmnrand.jl")
include("ordmnfit.jl")
include("ordmntest.jl")

end # module


text:
# --- User Parameters

# Number of samples to use to estimate probability distributions
sample_size = 10000;

# --- Imports

import Distributions
using Distributions: MvNormal, ContinuousUnivariateDistribution, cdf
using HypothesisTests
using LinearAlgebra
using Plots
using Random
using StatsBase: Histogram, fit

# --- Internal Parameters

# ------ Uniformity test parameters

# p-value cutoff
p_value_cutoff = 0.1;

# maximum number of attempts to generate a sample of vectors
max_attempts = 10;

# --- Probability distribution for θ (angle between vector and x_1-axis)

struct ThetaDistribution <: ContinuousUnivariateDistribution
    n::Integer
end

function Distributions.cdf(dist::ThetaDistribution, x::Real)
    # Preparations
    k = dist.n-2

    # Compute pdf
    value = 0
    if k % 2 == 0
        # Precompute re-used coefficients
        sgn = (-1)^(k÷2)
        k_choose_k_over_2 = binomial(big(k), k÷2)

        for j = 0:(k÷2 - 1)
            value += (-1)^j * binomial(big(k), j) * sin((k - 2*j) * x) / (k - 2*j)
        end
        value *= 2
        value += sgn * k_choose_k_over_2 * x  # contribution from i = k÷2 term
        value *= sgn / π / k_choose_k_over_2

    else
        norm_ = 0
        for j = 0:(k-1)÷2
            coef = (-1)^j * binomial(big(k), j) / (k - 2*j)
            value += coef * (1 - cos((k - 2*j) * x))
            norm_ += coef
        end
        value /= 2 * norm_

    end

    return value
end

# --- Tests

# 2D: cdf(x) = x / π
theta_dist_2d = ThetaDistribution(2)
for x = range(0, π; length=10)
    computed_cdf = cdf(theta_dist_2d, x)
    expected_cdf = x / π
    @assert (computed_cdf ≈ expected_cdf) "Expected: $(expected_cdf). Got: $(computed_cdf)"
end

# 3D: cdf(x) = 0.5 * (1 - cos(x))
theta_dist_3d = ThetaDistribution(3)
for x = range(0, π; length=10)
    computed_cdf = cdf(theta_dist_3d, x)
    expected_cdf = 0.5 * (1 - cos(x))
    @assert (computed_cdf ≈ expected_cdf) "Expected: $(expected_cdf). Got: $(computed_cdf)"
end

function generate_blades(n::Integer, s::Integer, sample_size::Integer)::Tuple
    """
    Generate a sample of s-dimensional blades in an n-dimensional space where each blade
    is defined by unit vectors drawn from a uniform distribution over the surface of an
    unit n-sphere.

    Parameters
    ----------
    n: dimension of space

    s: grade of blades

    sample_size: number of blades to include in sample

    Return values
    -------------
    vectors: sample of vectors

    uniformity_stats: results of tests for uniformity of sample on unit n-sphere
    """
    # --- Set up uniformity tests

    theta_dist = ThetaDistribution(n)
    AD_p_value = 0
    KS_p_value = 0

    # --- Generate sample of vectors

    num_attempts = 0
    vectors = undef
    while ( ((AD_p_value < p_value_cutoff) || (KS_p_value < p_value_cutoff)) &&
            (num_attempts < max_attempts) )

        # Generate new sample
        dist = MvNormal(zeros(n), ones(n))
        num_vectors = sample_size * s
        vectors = rand(dist, num_vectors)
        for i = 1:num_vectors
            vectors[:, i] /= norm(vectors[:, i])
        end

        # Perform test for uniformity on unit n-sphere
        thetas = map(i -> atan(norm(vectors[:, i][2:end]), vectors[:, i][1]), 1:num_vectors)

        # Anderson-Darling Test
        AD_p_value = pvalue(OneSampleADTest(thetas, theta_dist))

        # Kolmogorov-Smirnov Test
        KS_p_value = pvalue(ExactOneSampleKSTest(thetas, theta_dist))

        # Update attempt count
        num_attempts += 1
    end

    # --- Group vectors into blades

    blades = []
    for i = 1:sample_size
        push!(blades, vectors[:, (i-1)*s+1:i*s])
    end

    # --- Save uniformity test results

    uniformity_stats = (
        success=(num_attempts < max_attempts),
        num_attempts=num_attempts,
        AD_p_value=AD_p_value,
        KS_p_value=KS_p_value,
    )

    return blades, uniformity_stats
end

# --- Generate blades

# Parameters
n_2d = 2;  # dimension of space
s = 2;  # grade of blades

# Generate sample of blades
blades, uniformity_stats = generate_blades(n_2d, s, sample_size)
if !uniformity_stats[:success]
    println("FAILED to generated sample with sufficient uniformity on unit circle.")
    display(uniformity_stats)
end

# --- Compute blade volumes

signed_volumes = Vector{Float64}()
for i = 1:sample_size
    F = qr(blades[i])
    signed_volume = det(F.R)
    push!(signed_volumes, signed_volume)
end

volumes = abs.(signed_volumes);

# --- Compute distributions

# Histogram parameters
num_hist_bins = 50

# Plot histograms
hist_bins_signed = range(-1, 1; length=num_hist_bins)
plt_signed = plot(histogram(signed_volumes; bins=hist_bins_signed, normalize=true))

hist_bins = range(0, 1; length=num_hist_bins)
plt = plot(histogram(volumes; bins=hist_bins, normalize=true))

display(plt_signed)
display(plt)

# Compute empirical volume distributions
hist_signed = fit(Histogram, signed_volumes; nbins=num_hist_bins)
edf_signed = normalize(hist_signed; mode=:pdf);

hist = fit(Histogram, volumes; nbins=num_hist_bins)
edf = normalize(hist; mode=:pdf);

# --- Generate blades

# Parameters
n_3d = 3;  # dimension of space
s = 2;  # grade of blades

# Generate sample of blades
blades, uniformity_stats = generate_blades(n_3d, s, sample_size)
if !uniformity_stats[:success]
    println("FAILED to generated sample with sufficient uniformity on unit sphere.")
    display(uniformity_stats)
end

# --- Compute blade volumes

signed_volumes = Vector{Float64}()
for i = 1:sample_size
    F = qr(blades[i])
    signed_volume = det(F.R)
    push!(signed_volumes, signed_volume)
end

volumes = abs.(signed_volumes);

# --- Compute distributions

# Histogram parameters
num_hist_bins = 50

# Plot histograms
hist_bins_signed = range(-1, 1; length=num_hist_bins)
plt_signed = plot(histogram(signed_volumes; bins=hist_bins_signed, normalize=true))

hist_bins = range(0, 1; length=num_hist_bins)
plt = plot(histogram(volumes; bins=hist_bins, normalize=true))

display(plt_signed)
display(plt)

# Compute empirical volume distributions
hist_signed = fit(Histogram, signed_volumes; nbins=num_hist_bins)
edf_signed = normalize(hist_signed; mode=:pdf);

hist = fit(Histogram, volumes; nbins=num_hist_bins)
edf = normalize(hist; mode=:pdf);

# --- Generate blades

# Parameters
n = 100;  # dimension of space
s = 10;  # grade of blades

# Generate sample of vectors
blades, uniformity_stats = generate_blades(n, s, sample_size)
if !uniformity_stats[:success]
    println("FAILED to generated sample with sufficient uniformity on unit $(n)-sphere.")
    display(uniformity_stats)
end

# --- Compute blade volumes

signed_volumes = Vector{Float64}()
for i = 1:sample_size
    F = qr(blades[i])
    signed_volume = det(F.R)
    push!(signed_volumes, signed_volume)
end

volumes = abs.(signed_volumes);

# --- Compute distributions

# Histogram parameters
num_hist_bins = 50

# Plot histograms
hist_bins_signed = range(-1, 1; length=num_hist_bins)
plt_signed = plot(histogram(signed_volumes; bins=hist_bins_signed, normalize=true))

hist_bins = range(0, 1; length=num_hist_bins)
plt = plot(histogram(volumes; bins=hist_bins, normalize=true))
display(plt_signed)
display(plt)

# Compute empirical volume distributions
hist_signed = fit(Histogram, signed_volumes; nbins=num_hist_bins)
edf_signed = normalize(hist_signed; mode=:pdf);

hist = fit(Histogram, volumes; nbins=num_hist_bins)
edf = normalize(hist; mode=:pdf);

# --- Generate blades

# Parameters
n = 500;  # dimension of space
s = 10;  # grade of blades

# Generate sample of vectors
blades, uniformity_stats = generate_blades(n, s, sample_size)
if !uniformity_stats[:success]
    println("FAILED to generated sample with sufficient uniformity on unit $(n)-sphere.")
    display(uniformity_stats)
end

# --- Compute blade volumes

signed_volumes = Vector{Float64}()
for i = 1:sample_size
    F = qr(blades[i])
    signed_volume = det(F.R)
    push!(signed_volumes, signed_volume)
end

volumes = abs.(signed_volumes);

# --- Compute distributions

# Histogram parameters
num_hist_bins = 50

# Plot histograms
hist_bins_signed = range(-1, 1; length=num_hist_bins)
plt_signed = plot(histogram(signed_volumes; bins=hist_bins_signed, normalize=true))

hist_bins = range(0, 1; length=num_hist_bins)
plt = plot(histogram(volumes; bins=hist_bins, normalize=true))
display(plt_signed)
display(plt)

# Compute empirical volume distributions
hist_signed = fit(Histogram, signed_volumes; nbins=num_hist_bins)
edf_signed = normalize(hist_signed; mode=:pdf);

hist = fit(Histogram, volumes; nbins=num_hist_bins)
edf = normalize(hist; mode=:pdf);

max_stars_repo_name: ktchu/RESEARCH-math-exploration

max_issues_repo_name: ktchu/RESEARCH-math-exploration

max_forks_repo_name: ktchu/RESEARCH-math-exploration


text:
# Load modules
using Plots, LaTeXStrings
pyplot()
using ThinFilmsTools
using Optim

##
function glass_transmittance(beam, incident, emergent)
    # Create a measured spectrum of reflection
    layers = [
        LayerTMMO(incident),
        LayerTMMO(RI.looyenga([0.5 0.5],[RIdb.air(beam.λ) RIdb.silicon(beam.λ)]); d=100.),
        LayerTMMO(RIdb.glass(beam.λ./1e3); d=300.),
        LayerTMMO(RIdb.glass(beam.λ./1e3); d=100.),
        LayerTMMO(emergent),
    ]
    sol = tmm_optics(beam, layers)[1]
    return vec(beam.p.*sol.Tp .+ (1.0 - beam.p).*sol.Ts)
end
##

# Wavelength range [nm]
λ = 400:1000
# Angle of incidence [degrees]
θ = [5.]
# Polarisation (1.0 = p, 0.0 = s, between 0.0 and 1.0 = average)
pol = 0.5
beam = PlaneWave(λ, θ; p=pol)

# Refractive indices of incident (0) and substrate (2)
incident = RIdb.air(beam.λ)
emergent = RIdb.silicon(beam.λ)

# Define the RI model to use
layers = [
    LayerTMMO(incident), # 1
    ModelFit(:looyenga; N=(ninc=incident, nhost=emergent)), # 2
    LayerTMMO(RIdb.glass(beam.λ./1e3); d=300.), # 3
    ModelFit(:cauchyurbach), # 4
    LayerTMMO(emergent), # 5
]

# Create transmittance spectrum to fit
Texp = glass_transmittance(beam, incident, emergent)
# plot(Spectrum1D(), beam.λ, Texp)
# gui()

options = Optim.Options(
    g_abstol=1e-8, g_reltol=1e-8, iterations=10^5, show_trace=true, store_trace=true,
);

# Seeds of the ModelFit that appears in layers. The order here matters and
# the elements in each seed must match the number of elements accepted in the
# ModelFit layer used
seed = [
    [100.0, 0.5], # intial guess for :looyenga model, layers[2]
    [100.0, 1.5, 0.12, 0.5, 0.0, 0.6, 1.2], # intial guess for :cauchyurbach model, layers[4]
]

# solOptim = FitTMMOptics(
#     Transmittance(Texp), seed, beam, layers;
#     options=options, alg=SAMIN(), lb=0.5.*seed, ub=1.5.*seed,
# )
solOptim = fit_tmm_optics(
    Transmittance(Texp), seed, beam, layers;
    options=options, alg=LBFGS(),
)

plot(FitSpectrum(),
    solOptim.beam.λ, solOptim.spectrumExp, solOptim.spectrumFit,
    xaxis=("Wavelength [nm]"),
    yaxis=("Transmittance"),
)
gui()

# julia> solOptim.objfunMin
# 2.0852427777269855e-7
#
# julia> solOptim.optParams
# 2-element Array{Array{Float64,1},1}:
#  [100.20261188423837, 0.5009380743144798]
#  [100.01313444247826, 1.5059167041688846, 0.1150000000000001, 0.5000000000000016, 0.0, 0.6175516492383015, 1.23581135410979]


text:
begin
pushfirst!(LOAD_PATH, "@stdlib")
import Pkg
popfirst!(LOAD_PATH)

# We need to Pkg.instantiate the package environment that this notebook worker process will launch in
local my_dir = @__DIR__
local pluto_dir = joinpath(my_dir, "..", "..")

local runner_env_dir = mktempdir()
local new_ptoml_path = joinpath(runner_env_dir, "Project.toml")

local ptoml_contents = read(joinpath(pluto_dir, "Project.toml"), String)
write(new_ptoml_path, ptoml_contents[findfirst("[deps]", ptoml_contents)[1]:end])

local pkg_ctx = Pkg.Types.Context(env=Pkg.Types.EnvCache(new_ptoml_path))

try
    Pkg.instantiate(pkg_ctx; io=devnull)
catch
    # if it failed, do it again without suppressing io
    Pkg.instantiate(pkg_ctx)
end

pushfirst!(LOAD_PATH, runner_env_dir)

#
include(joinpath(my_dir, "PlutoRunner.jl"))
#

popfirst!(LOAD_PATH)
end


text:
using MLJExampleInterface # substitute for correct interface pkg name
using Test
using MLJBase
import Distributions
using StableRNGs # for RNGs stable across all julia versions
rng = StableRNGs.StableRNG(123)

@testset "cool classifier" begin
    n = 100
    p = 3
    nclasses = 5
    X, y = make_blobs(n, p, centers=nclasses, rng=rng);
    model = MLJExampleInterface.CoolProbabilisticClassifier()
    mach = machine(model, X, y)
    fit!(mach, rows=1:2, verbosity=0)
    yhat = predict(mach, rows=3)
    @test size(Distributions.pdf(yhat, levels(y)) ) == (1, nclasses)

    e = evaluate!(mach, measure=BrierLoss(), verbosity=0)
    @test e.measurement[1] < 1.0

    @test fitted_params(mach).classes_seen_in_training == levels(y)
    @test report(mach).n_classes_seen == nclasses
end


text:
# protect minimum precision
minprecerror() = throw(DomainError("minimum precision is 24 bits"))

for T in (:Float32, :Float64, :Int16, :Int32, :Int64, :Int128, :BigFloat)
  for N in collect(0:23)
    @eval ArbFloat{$N}(x::$T) = minprecerror()
    @eval ArbReal{$N}(x::$T) = minprecerror()
    @eval ArbComplex{$N}(x::$T) = minprecerror()
    @eval ArbComplex{$N}(x::$T, y::$T) = minprecerror()
  end
end


text:
"Forms groups of components."

@inline indices(g::AbstractGroup) = g.indices

@inline function Base.iterate(g::AbstractGroup, state=1)
    state > length(g) && return nothing
    return g.indices.packed[state], state+1
end

@inline indices_iterator(g::AbstractGroup) = g

"Groups components, creating a vector of the shared entity ids."
struct BasicGroup <: AbstractGroup
    component_ids::NTuple
    indices::Indices
end

BasicGroup(cs) = BasicGroup(map(x -> component_id(eltype(x)), cs), Indices(shared_entity_ids(cs)))

@inline Base.in(c::Int, g::BasicGroup) = c ∈ g.indices
Base.length(bg::BasicGroup) = length(bg.indices)

"Groups components, and makes sure that the order of entities is the same in each of the components."
mutable struct OrderedGroup{CT,N} <: AbstractGroup
    components::CT
    component_ids::NTuple{N, Int}
    indices::Indices
    len::Int
    parent::Union{Nothing, OrderedGroup}
    child::Union{Nothing, OrderedGroup}
end

@inline Base.in(c::Entity, g::OrderedGroup) = c.id ∈ g.indices && g.indices[c.id] <= length(g)
@inline Base.in(::Type{T}, g::OrderedGroup) where {T} = component_id(T) ∈ g.component_ids

function OrderedGroup(cs, parent, child)
    valid_entities = shared_entity_ids(cs)
    for (datid, e) in enumerate(valid_entities)
        for c in cs
            ensure_entity_id!(c, e, datid)
        end
    end
    return OrderedGroup(cs, map(x -> component_id(eltype(x)), cs), cs[1].indices, length(valid_entities), parent, child)
end

Base.length(g::OrderedGroup) = g.len

Base.in(c::AbstractComponent, g::OrderedGroup) = c in g.components

group(m::AbstractLedger,  cs::Type{<:ComponentData}...; ordered=true) = group(m, map(x -> m[x], cs))

function group(m::AbstractLedger, comps)
    for g in groups(m)
        if !(g isa OrderedGroup)
            continue
        end
        t_g = group(g, comps)
        if t_g !== nothing
            return t_g
        end
    end

    parent_group = findfirst(g -> all(in(g), comps), groups(m))

    if parent_group !== nothing
        parent = find_lowest_parent(groups(m)[parent_group], comps)
        if length(parent.components) == length(comps) || (parent.child !== nothing && length(parent.child.components) == length(comps))
            throw(ArgumentError("An ordered group with at least one but not all of the components $(eltype.(comps)) already exists.\nAn ordered component group can only be a subgroup of a previously ordered group."))
        end

        child  = deepcopy(parent.child)
        g = OrderedGroup(comps, parent, child)
        parent.child = g
        return g
    else
        child_group = findfirst(g -> all(in(comps), g.components) && g isa OrderedGroup, groups(m))
        if child_group === nothing
            any_group = findfirst(g -> any(in(g.components), comps), groups(m))
            if any_group !== nothing

                throw(ArgumentError("An ordered group with at least one but not all of the components $(eltype.(comps)) already exists.\nAn ordered component group can only be a subgroup of a previously ordered group."))
            end
            g = OrderedGroup(comps, nothing, nothing)
            push!(groups(m), g)
        else
            g = OrderedGroup(comps, nothing, groups(m)[child_group])
            g.child.parent = g
            groups(m)[child_group] = g
        end
        return g
    end
end

function find_lowest_parent(g::OrderedGroup, cs)
    if length(g.components) < length(cs)
        return g.parent
    elseif g.child === nothing
        return g
    else
        return find_lowest_parent(g.child, cs)
    end
end

function find_lowest_child(g::OrderedGroup, ::Type{T}) where {T}
    if g.child !== nothing && T in g.child
        return find_lowest_child(g.child, T)
    else
        return g
    end
end

function group(g::OrderedGroup, cs)
    if all(in(g), cs) && length(g.components) == length(cs)
        return g
    else
        return group(g.child, cs)
    end
end

group(::Nothing, cs) = nothing

@generated function Base.getindex(g::OrderedGroup{CT}, ::Type{T}) where {CT,T<:ComponentData}
    tid = findfirst(x -> eltype(x) == T, CT.parameters)
    return :(g.components[$tid])
end

Base.@propagate_inbounds @inline function Base.getindex(g::OrderedGroup, e::Entity)
    @boundscheck if !(e in g)
        throw(BoundsError(g, e))
    end
    id = g.indices[e.id]
    return map(x -> x[id], g.components)
end

function register_new!(g::OrderedGroup, e::Entity)
    if all(map(x -> e in x && ensure_entity_id!(x, e.id, g.len + 1), g.components))
        g.len += 1
    end
    if g.parent !== nothing
        register_new!(g.parent, e)
    end
end

@inline function Base.setindex!(g::OrderedGroup, v::T, e::Entity) where {T<:ComponentData}
    comp = g[T]
    @boundscheck if !(e in comp)
        comp[e] = v
        register_new!(find_lowest_child(g, T), e)
        return v
    end
    return @inbounds comp[e] = v
end


text:
## Lotka-Volterra

function lotka(du,u,p,t)
  x = u[1]
  y = u[2]
  du[1] = p[1]*x - p[2]*x*y
  du[2] = -p[3]*y + p[4]*x*y
end

"""
Lotka-Voltera Equations (Non-stiff)

```math
\\frac{dx}{dt} = ax - bxy
```
```math
\\frac{dy}{dt} = -cy + dxy
```

with initial condition ``x=y=1``
"""
prob_ode_lotkavoltera = ODEProblem(lotka,[1.0,1.0],(0.0,1.0),[1.5,1.0,3.0,1.0])

## Fitzhugh-Nagumo

function fitz(du,u,p,t)
  v = u[1]
  w = u[2]
  a = p[1]
  b = p[2]
  τinv = p[3]
  l = p[4]
  du[1] = v - v^3/3 -w + l
  du[2] = τinv*(v +  a - b*w)
end
"""
Fitzhugh-Nagumo (Non-stiff)

```math
\\frac{dv}{dt} = v - \\frac{v^3}{3} - w + I_{est}
```
```math
τ \\frac{dw}{dt} = v + a -bw
```

with initial condition ``v=w=1``
"""
prob_ode_fitzhughnagumo = ODEProblem(fitz,[1.0;1.0],(0.0,1.0),(0.7,0.8,1/12.5,0.5))

#Van der Pol Equations
@parameters t μ
@variables x(t) y(t)
@derivatives D'~t
eqs = [D(x) ~ μ*((1-x^2)*y - x),
       D(y) ~ y]
de = ODESystem(eqs)
van = ODEFunction(de, [x,y], [μ], jac=true, Wfact=true)

"""
Van der Pol Equations

```math
\\frac{dx}{dt} = y
```
```math
\\frac{dy}{dt} = μ((1-x^2)y -x)
```

with ``μ=1.0`` and ``u_0=[0,\\sqrt{3}]``

Non-stiff parameters.
"""
prob_ode_vanderpol = ODEProblem(van,[0;sqrt(3)],(0.0,1.0),1.0)

"""
Van der Pol Equations

```math
\\frac{dx}{dt} = y
```
```math
\\frac{dy}{dt} = μ(1-x^2)y -x
```

with ``μ=10^6`` and ``u_0=[0,\\sqrt{3}]``

Stiff parameters.
"""
prob_ode_vanstiff = ODEProblem(van,[0;sqrt(3)],(0.0,1.0),1e6)

# ROBER
@parameters t k₁ k₂ k₃
@variables y₁(t) y₂(t) y₃(t)
@derivatives D'~t
eqs = [D(y₁) ~ -k₁*y₁+k₃*y₂*y₃,
       D(y₂) ~ k₁*y₁-k₂*y₂^2-k₃*y₂*y₃,
       D(y₃) ~ k₂*y₂^2]
de = ODESystem(eqs)
rober = ODEFunction(de, [y₁,y₂,y₃], [k₁,k₂,k₃], jac=true, Wfact=true)

"""
The Robertson biochemical reactions: (Stiff)

```math
\\frac{dy₁}{dt} = -k₁y₁+k₃y₂y₃
```
```math
\\frac{dy₂}{dt} =  k₁y₁-k₂y₂^2-k₃y₂y₃
```
```math
\\frac{dy₃}{dt} =  k₂y₂^2
```

where ``k₁=0.04``, ``k₂=3\\times10^7``, ``k₃=10^4``. For details, see:

Hairer Norsett Wanner Solving Ordinary Differential Equations I - Nonstiff Problems Page 129

Usually solved on ``[0,1e11]``
"""
prob_ode_rober = ODEProblem(rober,[1.0;0.0;0.0],(0.0,1e11),(0.04,3e7,1e4))

# Three Body
const threebody_μ = big(0.012277471); const threebody_μ′ = 1 - threebody_μ

threebody = (du,u,p,t) -> begin
  # 1 = y₁
  # 2 = y₂
  # 3 = y₁'
  # 4 = y₂'
  D₁ = ((u[1]+threebody_μ)^2 + u[2]^2)^(3/2)
  D₂ = ((u[1]-threebody_μ′)^2 + u[2]^2)^(3/2)
  du[1] = u[3]
  du[2] = u[4]
  du[3] = u[1] + 2u[4] - threebody_μ′*(u[1]+threebody_μ)/D₁ - threebody_μ*(u[1]-threebody_μ′)/D₂
  du[4] = u[2] - 2u[3] - threebody_μ′*u[2]/D₁ - threebody_μ*u[2]/D₂
end

@doc doc"""
The ThreeBody problem as written by Hairer: (Non-stiff)

```math
\frac{dy₁}{dt} = y₁ + 2\frac{dy₂}{dt} - \bar{μ}\frac{y₁+μ}{D₁} - μ\frac{y₁-\bar{μ}}{D₂}
```
```math
\frac{dy₂}{dt} = y₂ - 2\frac{dy₁}{dt} - \bar{μ}\frac{y₂}{D₁} - μ\frac{y₂}{D₂}
```
```math
D₁ = ((y₁+μ)^2 + y₂^2)^{3/2}
```
```math
D₂ = ((y₁-\bar{μ})^2+y₂^2)^{3/2}
```
```math
μ = 0.012277471
```
```math
\bar{μ} =1-μ
```

From Hairer Norsett Wanner Solving Ordinary Differential Equations I - Nonstiff Problems Page 129

Usually solved on ``t₀ = 0.0`` and ``T = 17.0652165601579625588917206249``
Periodic with that setup.
"""
prob_ode_threebody = ODEProblem(threebody,[0.994, 0.0, 0.0, big(-2.00158510637908252240537862224)],(big(0.0),big(17.0652165601579625588917206249)))

# Rigid Body Equations

@parameters t I₁ I₂ I₃
@variables y₁(t) y₂(t) y₃(t)
@derivatives D'~t
eqs = [D(y₁) ~ I₁*y₂*y₃,
       D(y₂) ~ I₂*y₁*y₃,
       D(y₃) ~ I₃*y₁*y₂]
de = ODESystem(eqs)
rigid = ODEFunction(de, [y₁,y₂,y₃], [I₁,I₂,I₃], jac=true, Wfact=true)

"""
Rigid Body Equations (Non-stiff)

```math
\\frac{dy₁}{dt}  = I₁y₂y₃
```
```math
\\frac{dy₂}{dt}  = I₂y₁y₃
```
```math
\\frac{dy₃}{dt}  = I₃y₁y₂
```

with ``I₁=-2``, ``I₂=1.25``, and ``I₃=-1/2``.

The initial condition is ``y=[1.0;0.0;0.9]``.

From Solving Differential Equations in R by Karline Soetaert

or Hairer Norsett Wanner Solving Ordinary Differential Equations I - Nonstiff Problems Page 244

Usually solved from 0 to 20.
"""
prob_ode_rigidbody = ODEProblem(rigid,[1.0,0.0,0.9],(0.0,20.0),(-2.0,1.25,-0.5))

# Pleiades Problem

pleiades = (du,u,p,t) -> begin
  x = view(u,1:7)   # x
  y = view(u,8:14)  # y
  v = view(u,15:21) # x′
  w = view(u,22:28) # y′
  du[1:7] .= v
  du[8:14].= w
  for i in 14:28
    du[i] = zero(eltype(u))
  end
  for i=1:7,j=1:7
    if i != j
      r = ((x[i]-x[j])^2 + (y[i] - y[j])^2)^(3/2)
      du[14+i] += j*(x[j] - x[i])/r
      du[21+i] += j*(y[j] - y[i])/r
    end
  end
end

@doc doc"""
Pleiades Problem (Non-stiff)

```math
\frac{d^2xᵢ}{dt^2} = \sum_{j≠i} mⱼ(xⱼ-xᵢ)/rᵢⱼ
```
```math
\frac{d^2yᵢ}{dt^2} = \sum_{j≠i} mⱼ(yⱼ-yᵢ)/rᵢⱼ
```

where

```math
rᵢⱼ = ((xᵢ-xⱼ)^2 + (yᵢ-yⱼ)^2)^{3/2}
```

and initial conditions are

```math
x₁(0) = 3
```
```math
x₂(0) = 3
```
```math
x₃(0) = -1
```
```math
x₄(0) = -3
```
```math
x₅(0) = 2
```
```math
x₆(0) = -2
```
```math
x₇(0) = 2
```
```math
y₁(0) = 3
```
```math
y₂(0) = -3
```
```math
y₃(0) = 2
```
```math
y₄(0) = 0
```
```math
y₅(0) = 0
```
```math
y₆(0) = -4
```
```math
y₇(0) = 4
```

and with ``\frac{dxᵢ(0)}{dt}=\frac{dyᵢ(0)}{dt}=0`` except for

```math
\frac{dx₆(0)}{dt} = 1.75
```
```math
\frac{dx₇(0)}{dt} = -1.5
```
```math
\frac{dy₄(0)}{dt} = -1.25
```
```math
\frac{dy₅(0)}{dt} = 1
```

From Hairer Norsett Wanner Solving Ordinary Differential Equations I - Nonstiff Problems Page 244

Usually solved from 0 to 3.
"""
prob_ode_pleiades = ODEProblem(pleiades,[3.0,3.0,-1.0,-3.0,2.0,-2.0,2.0,3.0,-3.0,2.0,0,0,-4.0,4.0,0,0,0,0,0,1.75,-1.5,0,0,0,-1.25,1,0,0],(0.0,3.0))

Random.seed!(100)
const mm_A = rand(4,4)
mm_linear = function (du,u,p,t)
  mul!(du,mm_A,u)
end
const MM_linear =Matrix(Diagonal(0.5ones(4)))
mm_f = ODEFunction(mm_linear;analytic = (u0,p,t) -> exp(inv(MM_linear)*mm_A*t)*u0, mass_matrix=MM_linear)
prob_ode_mm_linear = ODEProblem(mm_f,rand(4),(0.0,1.0))

@parameters t p1 p2 p3 p4 p5 p6 p7 p8 p9 p10 p11 p12
@variables y1(t) y2(t) y3(t) y4(t) y5(t) y6(t) y7(t) y8(t)
@derivatives D'~t
eqs = [D(y1) ~ -p1*y1 + p2*y2 + p3*y3 + p4,
       D(y2) ~ p1*y1 - p5*y2,
       D(y3) ~ -p6*y3 + p2*y4 + p7*y5,
       D(y4) ~ p3*y2 + p1*y3 - p8*y4,
       D(y5) ~ -p9*y5 + p2*y6 + p2*y7,
       D(y6) ~ -p10*y6*y8 + p11*y4 + p1*y5 -
                p2*y6 + p11*y7,
       D(y7) ~  p10*y6*y8 - p12*y7,
       D(y8) ~ -p10*y6*y8 + p12*y7]
de = ODESystem(eqs)
hires = ODEFunction(de, [y1,y2,y3,y4,y5,y6,y7,y8],
                        [p1,p2,p3,p4,p5,p6,p7,p8,p9,p10,p11,p12],
                        jac=true, Wfact=true)

u0 = zeros(8)
u0[1] = 1
u0[8] = 0.0057

"""
Hires Problem (Stiff)

It is in the form of

```math
\\frac{dy}{dt} = f(y)
```

 with

```math
 y(0)=y_0, \\quad y \\in ℝ^8, \\quad 0 ≤ t ≤ 321.8122
```

where ``f`` is defined by

``f(y) = \\begin{pmatrix} −1.71y_1 & +0.43y_2 & +8.32y_3 & +0.0007y_4 & \\\\ 1.71y_1 & −8.75y_2 & & & \\\\ −10.03y_3 & +0.43y_4 & +0.035y_5 & & \\\\ 8.32y_2 & +1.71y_3 & −1.12y_4 & & \\\\ −1.745y_5 & +0.43y_6 & +0.43y_7 & & \\\\ −280y_6y_8 & +0.69y_4 & +1.71y_5 & −0.43y_6 & +0.69y_7 \\\\ 280y_6y_8 & −1.81y_7 & & & \\\\ −280y_6y_8 & +1.81y_7 & & & \\end{pmatrix}``

Reference: [demohires.pdf](http://www.radford.edu/~thompson/vodef90web/problems/demosnodislin/Demos_Pitagora/DemoHires/demohires.pdf)
Notebook: [Hires.ipynb](http://nbviewer.jupyter.org/github/JuliaDiffEq/DiffEqBenchmarks.jl/blob/master/StiffODE/Hires.ipynb)
"""
prob_ode_hires = ODEProblem(hires,u0,(0.0,321.8122), (1.71, 0.43, 8.32, 0.0007, 8.75,
                                                      10.03, 0.035, 1.12, 1.745, 280.0,
                                                      0.69, 1.81))

@parameters t p1 p2 p3
@variables y1(t) y2(t) y3(t)
@derivatives D'~t
eqs = [D(y1) ~ p1*(y2+y1*(1-p2*y1-y2)),
       D(y2) ~ (y3-(1+y1)*y2)/p1,
       D(y3) ~ p3*(y1-y3)]
de = ODESystem(eqs)
jac = calculate_jacobian(de)
ModelingToolkit.calculate_factorized_W(de)
orego = ODEFunction(de, [y1,y2,y3], [p1,p2,p3], jac=true, Wfact=true)

"""
Orego Problem (Stiff)

It is in the form of ``\\frac{dy}{dt}=f(y), \\quad y(0)=y0,`` with

```math
y \\in ℝ^3, \\quad 0 ≤ t ≤ 360
```

where ``f`` is defined by

``f(y) = \\begin{pmatrix} s(y_2 - y_1(1-qy_1-y_2)) \\\\ (y_3 - y_2(1+y_1))/s \\\\ w(y_1-y_3) \\end{pmatrix}``

where ``s=77.27``, ``w=0.161`` and ``q=8.375⋅10^{-6}``.

Reference: [demoorego.pdf](http://www.radford.edu/~thompson/vodef90web/problems/demosnodislin/Demos_Pitagora/DemoOrego/demoorego.pdf)
Notebook: [Orego.ipynb](http://nbviewer.jupyter.org/github/JuliaDiffEq/DiffEqBenchmarks.jl/blob/master/StiffODE/Orego.ipynb)
"""
prob_ode_orego = ODEProblem(orego,[1.0,2.0,3.0],(0.0,30.0),[77.27,8.375e-6,0.161])


text:
# Traits
"Hard inference procedure that expects Boolean predicate"
struct IsHard{ALG} end

"Soft Inference procedure"
struct IsSoft{ALG} end


text:
module ReachabilityModels

using Reexport
@reexport using ReachabilityAnalysis

include("utils.jl")

# name alias
fetch_model = load_model
fetch_meta = load_meta

export fetch_model, fetch_meta,
       load_model, load_meta,
       list_models,
       @relpath

end # module

max_stars_repo_name: mforets/ReachabilityLibrary

max_issues_repo_name: mforets/ReachabilityLibrary

max_forks_repo_name: mforets/ReachabilityLibrary


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

using KernelMethods
using SimilaritySearch: l2_distance
# import KernelMethods.Nets: KlusterClassifier
using Test

include("loaddata.jl")

@testset "KNN Regression" begin
    X, y = loadlinearreg()
    nnr = NearNeighborRegression(l2_distance, X, y)
    @test optimize!(nnr, pearson, folds=3)[1][1] > 0.95
    @test optimize!(nnr, spearman, folds=3)[1][1] > 0.95
    @show optimize!(nnr, isqerror, folds=5)
    #@test optimize!(nnr, accuracy, runs=5, trainratio=0.3, testratio=0.3)[1][1] > 0.9
    #@test sum([maximum(x) for x in predict_proba(nnr, X, smoothing=0)])/ length(X) > 0.9 ## close to have all ones, just in case
    #@test sum([maximum(x) for x in predict_proba(nnr, X, smoothing=0.01)])/ length(X) > 0.9 ## close to have all ones, just in case
end


text:
#region AbstractModel Deconstruction & Reconstruction
#########################################################
######### MODEL DECONSTRUCTION & RECONSTRUCTION #########
# Deconstruct a model (which is naturally hierarchical)
# into an array representation.
# Alter array representation.
# Reconstruct hierarchical model.
#########################################################

abstract type AbstractDeconstruction{OBJ} end

struct Deconstruction{OBJ} <: AbstractDeconstruction{OBJ}
    source::OBJ
    deconstruction::Tuple{Type,Array}
end

"""
    set_deep_dx!(model_deconstruction, dxs, val)

Within an object deconstructed into an array, set one element to `val`,
where `dxs` indicate that element by successive indexing into nested arrays.
"""
function set_deep_dx!(model_deconstruction::Tuple{Type,Array}, dxs, val)
    tmp = model_deconstruction
    for dx in dxs[1:end-1]
        tmp = tmp[2][dx]
    end
    tmp[2][dxs[end]] = val
end

abstract type AbstractDeconstructor end
struct Deconstructor <: AbstractDeconstructor end

function deconstruct(val::Union{AbstractString,Number}, ad::AbstractDeconstructor=Deconstructor())
    return (typeof(val), val)
end

function deconstruct(arr::AA, ad::AbstractDeconstructor=Deconstructor()) where {AA<:AbstractArray}
    deconstruction = map(ad, arr)
    return (base_type(typeof(arr)), [deconstruction...]) # remade incase static
end

function deconstruct(tup::Tuple, ad::AbstractDeconstructor=Deconstructor())
    deconstruction = map(ad, tup)
    return (base_type(typeof(tup)), [deconstruction...])
end

function base_type(T::Int)
    T
end

function base_type(::Type{T}) where T <: Real
    return T
end

function base_type(::Type{Array{T,N}}) where {T, N}
    BT = base_type(T)
    return Array{BT,N}
end

function base_type(::Type{Tuple{T,S}}) where {T,S}
    BT = base_type(T)
    BS = base_type(S)
    return Tuple{BT,BS}
end

# function base_type(::Type{SA}) where {N,M,T,TUP,SA<:Union{SArray{TUP,T,N,M},MArray{TUP,T,N,M},SizedArray{TUP,T,N,M}}}
#     BT = base_type(T)
#     return SArray{TUP,BT,N,M}
# end

function reconstruct(tup::Tuple{Type,<:Union{Number,AbstractString}})
    return tup[2]
end

function reconstruct(tup::Tuple{Type,Array})
    typ, arr = tup
    base_typ = base_type(typ)
    if typ <: Union{AbstractArray, Tuple}
        return base_typ(reconstruct.(arr))
    else
        return base_typ(reconstruct.(arr)...)
    end
end

############################################
##endregion


text:
# run simple tests on the filter and updates

using SubmodularMaximization
using LinearAlgebra
using HistogramFilters

grid_size = 2
grid = Grid(grid_size, grid_size)

robot_state = div.((grid_size, grid_size), 2)
sensor = RangingSensor()

grid_states = get_states(grid)
histogram_filter = Filter(grid)

@show prior = HistogramFilters.get_data(histogram_filter)

if !all(prior .== 1/grid_size^2)
  println("Problem with uniform prior")
end

# Test neighbors
@show ns = neighbors(grid, (1,1))
real_ns = Set([(1,1), (1,2), (2,1)])
if real_ns != Set(ns)
  println("Neighbors do not match")
end

@show transition = Matrix(transition_matrix(grid))

# true 4x4 transition matrix for a 2x2 grid noting that the indices are
# [1, 3]
# [2, 4]
real_transition = vcat([1/3  1/3  1/3  0  ],
                       [1/3  1/3    0  1/3],
                       [1/3    0  1/3  1/3],
                       [  0  1/3  1/3  1/3])

if transition != real_transition
  println("Problem with transition matrix")
end

max_stars_repo_name: mcorah/MultiAgentSensing

max_issues_repo_name: mcorah/MultiAgentSensing

max_forks_repo_name: mcorah/MultiAgentSensing