text:
Eye4f0() = Mat4f0(Matrix(1.0f0I, 4, 4))
#Came from GLAbstraction/GLMatrixmath.jl
function scalemat(s::Vec{3,T}) where {T}
    T0, T1 = zero(T), one(T)
    return Mat{4}(s[1], T0, T0, T0, T0, s[2], T0, T0, T0, T0, s[3], T0, T0, T0, T0, T1)
end

function scalemat(s::T) where {T}
    T0, T1 = zero(T), one(T)
    return Mat{4}(s, T0, T0, T0, T0, s, T0, T0, T0, T0, s, T0, T0, T0, T0, T1)
end

scalemat(s::Point{3,T}) where {T} = scalemat(convert(Vec{3,T}, s))

translmat_x(x::T) where {T} = translmat(Vec{3,T}(x, 0, 0))
translmat_y(y::T) where {T} = translmat(Vec{3,T}(0, y, 0))
translmat_z(z::T) where {T} = translmat(Vec{3,T}(0, 0, z))

function translmat(t::Vec{3,T}) where {T}
    T0, T1 = zero(T), one(T)
    return Mat{4}(T1, T0, T0, T0, T0, T1, T0, T0, T0, T0, T1, T0, t[1], t[2], t[3], T1)
end
translmat(t::Point{3,T}) where {T} = translmat(convert(Vec{3,T}, t))

function rotate(a::Number, axis::Vec{3,T}) where {T}
    axis = normalize(axis)
    u, v, w = axis
    return Mat4{T}(u^2 + (1 - u^2) * cos(a), u * v * (1 - cos(a)) + w * sin(a),
                   u * w * (1 - cos(a)) - v * sin(a), 0.0,
                   u * v * (1 - cos(a)) - w * sin(a), v^2 + (1 - v^2) * cos(a),
                   v * w * (1 - cos(a)) + u * sin(a), 0.0,
                   u * w * (1 - cos(a)) + v * sin(a), v * w * (1 - cos(a)) - u * sin(a),
                   w^2 + (1 - w^2) * cos(a), 0.0, 0.0, 0.0, 0.0, 1.0)
end

function rotate(::Type{T}, angle::Number, axis::Vec{3}) where {T}
    return rotate(T(angle), convert(Vec{3,T}, axis))
end

function rotate(v1::T, angle::Number, axis) where {T<:StaticArray{Tuple{3}}}
    return T((rotate(angle, axis)*SVector{4}(v1..., 1.0))[1:3]...)
end
function rotate(v1::T, angle::Number, origin, axis) where {T<:StaticArray{Tuple{3}}}
    return T((translmat(origin)*rotate(angle, axis)*translmat(-origin)*SVector{4}(v1...,
                                                                                1.0))[1:3]...)
end

function rotate(v1::Vec3{T}, v2::Vec3{T}) where {T}
    vr = v2 - v1
    l = norm(vr)
    angle = acos(vr[3] / l)
    axis = normalize(cross(Vec3(0.0000000001f0, 0.0f0, 1.0f0), vr))
    if length(axis) > 0.0001
        return rotate(angle, axis)
    end
    return Mat4(Eye4f0(T, 4))
end

function rotmat_x(angle::T) where {T}
    T0, T1 = zero(T), one(T)
    return Mat{4}(T1, T0, T0, T0, T0, cos(angle), sin(angle), T0, T0, -sin(angle),
                  cos(angle), T0, T0, T0, T0, T1)
end
function rotmat_y(angle::T) where {T}
    T0, T1 = zero(T), one(T)
    return Mat{4}(cos(angle), T0, -sin(angle), T0, T0, T1, T0, T0, sin(angle), T0,
                  cos(angle), T0, T0, T0, T0, T1)
end
function rotmat_z(angle::T) where {T}
    T0, T1 = zero(T), one(T)
    return Mat{4}(cos(angle), sin(angle), T0, T0, -sin(angle), cos(angle), T0, T0, T0, T0,
                  T1, T0, T0, T0, T0, T1)
end
"""
    Create view frustum

    Parameters
    ----------
        left : float
         Left coordinate of the field of view.
        right : float
         Left coordinate of the field of view.
        bottom : float
         Bottom coordinate of the field of view.
        top : float
         Top coordinate of the field of view.
        znear : float
         Near coordinate of the field of view.
        zfar : float
         Far coordinate of the field of view.

    Returns
    -------
        M : array
         View frustum matrix (4x4).
"""
function frustum(left::T, right::T, bottom::T, top::T, znear::T, zfar::T) where {T}
    (right == left || bottom == top || znear == zfar) && return Mat4{T}(I)
    T0, T1, T2 = zero(T), one(T), T(2)
    return Mat4{T}(T2 * znear / (right - left), T0, T0, T0, T0, T2 * znear / (top - bottom),
                   T0, T0, (right + left) / (right - left), (top + bottom) / (top - bottom),
                   -(zfar + znear) / (zfar - znear), -T1, T0, T0,
                   -(T2 * znear * zfar) / (zfar - znear), T0)
end

"""
`proj = projmatpersp([T], fovy, aspect, znear, zfar)` defines
a projection matrix with a given angular field-of-view `fovy` along
the y-axis (measured in degrees), the specified `aspect` ratio, and
near and far clipping planes `znear`, `zfar`. Optionally specify the
element type `T` of the matrix.
"""
function projmatpersp(fovy::T, aspect::T, znear::T, zfar::T) where {T}
    (znear == zfar) && error("znear ($znear) must be different from tfar ($zfar)")
    h = T(tan(fovy / 180.0 * pi) * znear)
    w = T(h * aspect)
    return frustum(-w, w, -h, h, znear, zfar)
end

function projmatpersp(w::Integer, h::Integer, fovy::T, znear::T, zfar::T) where {T}
    return projmatpersp(fovy, T(w / h), znear, zfar)
end

function projmatpersp(::Type{T}, fovy::Number, aspect::Number, znear::Number,
                      zfar::Number) where {T}
    return projmatpersp(T(fovy), T(aspect), T(znear), T(zfar))
end
"""
`proj = projmatpersp([T], rect, fov, near, far)` defines the
projection ratio in terms of the rectangular view size `rect` rather
than the aspect ratio.
"""
function projmatpersp(wh::Area, fov::T, near::T, far::T) where {T}
    return projmatpersp(fov, T(wh.w / wh.h), near, far)
end
function projmatpersp(::Type{T}, wh::Area, fov::Number, near::Number, far::Number) where {T}
    return projmatpersp(T(fov), T(wh.w / wh.h), T(near), T(far))
end

"""
`view = lookatmat(eyeposition, lookatmat, up)` creates a view matrix with
the eye located at `eyeposition` and looking at position `lookatmat`,
with the top of the window corresponding to the direction `up`. Only
the component of `up` that is perpendicular to the vector pointing
from `eyeposition` to `lookatmat` will be used.  All inputs must be
supplied as 3-vectors.
"""
function lookatmat(eyePos::Union{Vec{3,T},Point{3,T}}, lookAt::Union{Vec{3,T},Point{3,T}},
                   up::Union{Vec{3,T},Point{3,T}}) where {T}
    zaxis  = normalize(eyePos - lookAt)
    xaxis  = normalize(cross(up, zaxis))
    yaxis  = normalize(cross(zaxis, xaxis))
    T0, T1 = zero(T), one(T)
    return Mat{4}(xaxis[1], yaxis[1], zaxis[1], T0, xaxis[2], yaxis[2], zaxis[2], T0,
                  xaxis[3], yaxis[3], zaxis[3], T0, T0, T0, T0, T1) * translmat(-eyePos)
end
function lookatmat(eyePos, lookAt::Union{Vec{3,T},Point{3,T}},
                   up::Union{Vec{3,T},Point{3,T}}) where {T}
    return lookatmat(Vec{3,T}(eyePos), lookAt, up)
end
function lookatmat(::Type{T}, eyePos, lookAt::Vec{3}, up::Vec{3}) where {T}
    return lookatmat(Vec{3,T}(eyePos), Vec{3,T}(lookAt), Vec{3,T}(up))
end

"""
    projmatortho(left::T, right::T, bottom::T, top::T) where T

2D orthographic projection.
"""
function projmatortho(left::T, right::T, bottom::T, top::T) where {T}
    out = zeros(T, 4, 4)
    out[1, 1] = 2 / (right - left)
    out[2, 2] = 2 / (top - bottom)
    out[3, 3] = -1
    out[4, 1] = -(right + left) / (right - left)
    out[4, 2] = -(top + bottom) / (top - bottom)
    out[4, 4] = 1
    return Mat4{T}(out')
end

function projmatortho(::Type{T}, left::Number, right::Number, bottom::Number,
                      top::Number) where {T}
    return projmatortho(T(left), T(right), T(bottom), T(top))
end

function projmatortho(wh::Area, near::T, far::T) where {T}
    return projmatortho(zero(T), T(wh.w), zero(T), T(wh.h), near, far)
end
function projmatortho(::Type{T}, wh::Area, near::Number, far::Number) where {T}
    return projmatortho(wh, T(near), T(far))
end
function projmatortho(w::Integer, h::Integer, near::T, far::T) where {T}
    return projmatortho(zero(T), T(w), zero(T), T(h), near, far)
end

"""
    projmatortho(left::T, right::T, bottom::T, top::T, znear::T, zfar::T) where T

3D orthographic projection.
"""
function projmatortho(left::T, right::T, bottom::T, top::T, znear::T, zfar::T) where {T}
    if right == left || bottom == top || znear == zfar
        return Mat4{T}(I)
    else
        out = zeros(T, 4, 4)
        out[1, 1] = 2 / (right - left)
        out[2, 2] = 2 / (top - bottom)
        out[3, 3] = -2 / (zfar - znear)
        out[4, 1] = -(right + left) / (right - left)
        out[4, 2] = -(top + bottom) / (top - bottom)
        out[4, 3] = -(zfar + znear) / (zfar - znear)
        out[4, 4] = 1
        return Mat4{T}(out')
    end
end

function projmatortho(::Type{T}, left::Number, right::Number, bottom::Number, top::Number,
                      znear::Number, zfar::Number) where {T}
    return projmatortho(T(left), T(right), T(bottom), T(top), T(znear), T(zfar))
end

function Base.:(*)(q::Quaternions.Quaternion{T}, v::Vec{3,T}) where {T}
    t = T(2) * cross(Vec(q.v1, q.v2, q.v3), v)
    return v + q.s * t + cross(Vec(q.v1, q.v2, q.v3), t)
end
function Quaternions.qrotation(axis::Vec{3,T}, theta) where {T<:Real}
    u = normalize(axis)
    s = sin(theta / 2)
    return Quaternions.Quaternion{T}(cos(theta / 2), s * u[1], s * u[2], s * u[3], true)
end

mutable struct Pivot{T}
    origin      :: Vec{3,T}
    xaxis       :: Vec{3,T}
    yaxis       :: Vec{3,T}
    zaxis       :: Vec{3,T}
    rotation    :: Quaternions.Quaternion
    translation :: Vec{3,T}
    scale       :: Vec{3,T}
end

GeometryBasics.origin(p::Pivot) = p.origin

rotmat4(q::Quaternions.Quaternion{T}) where {T} = Mat4{T}(q)

function (::Type{M})(q::Quaternions.Quaternion) where {M<:Mat4}
    T = eltype(M)
    sx, sy, sz = 2q.s * q.v1, 2q.s * q.v2, 2q.s * q.v3
    xx, xy, xz = 2q.v1^2, 2q.v1 * q.v2, 2q.v1 * q.v3
    yy, yz, zz = 2q.v2^2, 2q.v2 * q.v3, 2q.v3^2
    T0, T1 = zero(T), one(T)
    return Mat{4}(T1 - (yy + zz), xy + sz, xz - sy, T0, xy - sz, T1 - (xx + zz), yz + sx,
                  T0, xz + sy, yz - sx, T1 - (xx + yy), T0, T0, T0, T0, T1)
end

function (::Type{M})(q::Quaternions.Quaternion) where {M<:Mat3}
    T = eltype(M)
    sx, sy, sz = 2q.s * q.v1, 2q.s * q.v2, 2q.s * q.v3
    xx, xy, xz = 2q.v1^2, 2q.v1 * q.v2, 2q.v1 * q.v3
    yy, yz, zz = 2q.v2^2, 2q.v2 * q.v3, 2q.v3^2
    T0, T1 = zero(T), one(T)
    return Mat{3}(T1 - (yy + zz), xy + sz, xz - sy, xy - sz, T1 - (xx + zz), yz + sx,
                  xz + sy, yz - sx, T1 - (xx + yy))
end
function transfmat(p::Pivot)
    return translmat(p.origin) * #go to origin
           rotmat4(p.rotation) * #apply rotation
           translmat(-p.origin) * # go back to origin
           translmat(p.translation) #apply translation
end

function transfmat(translation, scale)
    T = eltype(translation)
    T0, T1 = zero(T), one(T)
    return Mat{4}(scale[1], T0, T0, T0, T0, scale[2], T0, T0, T0, T0, scale[3], T0,
                  translation[1], translation[2], translation[3], T1)
end

function transfmat(translation, scale, rotation::Quaternions.Quaternion)
    T = eltype(translation)
    trans_scale = transfmat(translation, scale)
    rotation = Mat4f0(rotation)
    return trans_scale * rotation
end
function transfmat(translation, scale, rotation::Vec{3,T}, up = Vec{3,T}(0, 0, 1)) where {T}
    q = rotation(rotation, up)
    return transfmat(translation, scale, q)
end

#Calculate rotation between two vectors
function rotation(u::Vec{3,T}, v::Vec{3,T}) where {T}
    # It is important that the inputs are of equal length when
    # calculating the half-way vector.
    u, v = normalize(u), normalize(v)
    # Unfortunately, we have to check for when u == -v, as u + v
    # in this case will be (0, 0, 0), which cannot be normalized.
    if (u == -v)
        # 180 degree rotation around any orthogonal vector
        other = (abs(dot(u, Vec{3,T}(1, 0, 0))) < 1.0) ? Vec{3,T}(1, 0, 0) :
                Vec{3,T}(0, 1, 0)
        return Quaternions.qrotation(normalize(cross(u, other)), T(180))
    end
    half = normalize(u + v)
    return Quaternions.Quaternion(dot(u, half), cross(u, half)...)
end


text:
@testset "set difference" begin
    comvals = Set([1])
    addvals = Set([2, 3])
    remvals = Set([4, 5])
    a = SetDifference(comvals, addvals, remvals)
    b = SetDifference(comvals, addvals, remvals)
    c = SetDifference(comvals, addvals, remvals)

    @testset "constructors" begin
        @test isa(a, SetDifference)
        @test a.comvals == comvals && a.addvals == addvals && a.remvals == remvals
    end

    @testset "equality operator" begin
        @test a == a
        @test a == b && b == a
        @test a == b && b == c && a == c
    end

    @testset "hash" begin
        @test hash(a) == hash(a)
        @test a == b && hash(a) == hash(b)
    end

    @testset "accessors" begin
        @test common(a) == comvals
        @test added(a) == addvals
        @test removed(a) == remvals
    end

    @testset "difference" begin
        a = Set([1, 2, 3, 3])
        b = Set([4, 2, 1])
        ei = Int[]
        @test diff(a, a) == SetDifference(Set([1, 2, 3]), Set(ei), Set(ei))
        @test diff(a, b) == SetDifference(Set([1, 2]), Set([4]), Set([3]))
        @test diff(b, a) == SetDifference(Set([2, 1]), Set([3]), Set([4]))
    end

    @testset "math" begin
        a = Set([1, 5, 3])
        b = Set([1, 4, 2, 3])
        d = SetDifference(Set([1, 3]), Set([4, 2]), Set([5]))
        @test a + d == b
        @test d + a == b
    end
end

ext: jl


text:
"""
    LogicalCombo <: Any

A LogicalCombo stores the variable names, domains, as well as an binary
representation of feasible values given constraints.

**Fields**

* keys : Stores the names of variables.
* domain : Stores the range of possible matches variables can have.
* logical : A binary vector marking feasibility of length equal to every
    possible combination of value which the matrix could take on.
* commandlist : An array collection of the strings input to generate the current
    state of the the object.

**Indexing**

* [x,y] Where x is numeric or : and y is numeric, colon, or symbol will
    return matrix point values regardless of feasibility.
* [:,:] Will collect the entire possible domain and is the same as `collect()`.
* [x,0] Will return the logical vector values
* [x,:,true] Will return the xth feasible value of the LogicalCombo set.
* [x,:,false] Will return the xth infeasible value of the LogicalCombo set.
"""
struct LogicalCombo
  keys::Array{Symbol}
  domain::AbstractArray
  logical::Array{Bool}
  lookup::Function
  type::String
  commands::Array{String}
end

LogicalCombo(keys, domain, logical) =
    LogicalCombo(keys, domain, logical, permutationlookup, "Permutation", String[])

LogicalCombo(keys, domain, logical, permutationlookup, type) =
    LogicalCombo(keys, domain, logical, permutationlookup, type, String[])

LogicalCombo() = LogicalCombo(Symbol[],[], Bool[])

function LogicalCombo(; kwargs...)
    if isempty(kwargs)
        return LogicalCombo(Symbol[],[], Bool[])
    else
        keys = []; domain = []
        for (kw, val) in kwargs;
            push!(keys, kw)
            push!(domain, val)
        end
    end
    LogicalCombo(keys, domain, fill(true,*(length.(domain)...)))
end

function LogicalCombo(x::Union{AbstractArray{Pair{Symbol, Any}}, Array{Pair{Symbol,UnitRange{Int64}},1}})
    # if length(x) == 0
    #     return LogicalCombo(Symbol[],[], Bool[0])
    # else
        mykeys = []; mydomain = []
        for (kw, val) in x;
            push!(mykeys, kw)
            push!(mydomain, val)
        end
    # end
    LogicalCombo(mykeys, mydomain, fill(true,*(length.(mydomain)...)))
end

Base.keys(x::LogicalCombo)     = x.keys
domain(x::LogicalCombo)        = x.domain
Base.size(x::LogicalCombo)     = [length(x.logical), length(x.keys)]
Base.size(x::LogicalCombo, y::Integer) = size(x)[y]

function permutationlookup(x::LogicalCombo, row::Integer, col::Integer)
    # Divisor is calculating based on how many values remains how many times to repeat the current value
    divisor = (col+1 > size(x)[2] ? 1 : x.domain[(col+1):size(x)[2]] .|> length |> prod)
    # indexvalue finds the index to select from the column of interest
    indexvalue = (Integer(ceil(row / divisor)) - 1) % length(domain(x)[col]) + 1
    domain(x)[col][indexvalue]
end

function permutationuniquelookup(x::LogicalCombo, row::Integer, col::Integer)
    n = size(x,2)
    fset = 1:n
    y = row
    z = 0
    setout = Integer[]
    for j in 1:col
        K = factorial(n-j)
        z = Integer(ceil(y/ K))
        y = y - (z-1) * K
        push!(setout, fset[z])
        fset = fset[[(1:(z-1))..., ((z+1):length(fset))...]]
    end
    domain(x)[setout[col]]
end

function Base.getindex(x::LogicalCombo, row::Integer, col::Integer)
    (col==0) && (row==0) && return :keys
    (row==0) && return keys(x)[col]
    (col==0) && return x.logical[row]
    x.lookup(x, row, col)
end

Base.collect(x::LogicalCombo; bool::Bool=true, varnames::Bool=true) =
  [x[i,j] for i in !varnames:size(x)[1], j in !bool:size(x)[2]]

function Base.getindex(x::LogicalCombo, row::Integer, col::Symbol)
  keymatch = findall(y -> y == col , x.keys)
  length(keymatch)==0 && throw("Symbol :$col not found")
  x[row, keymatch...]
end

Base.getindex(x::LogicalCombo, ::Colon) =  x.logical
Base.getindex(x::LogicalCombo, y::Union{Int64,UnitRange}) =  x.logical[y]
# Base.getindex(x::LogicalCombo, y::Union{Array}) =
#   [x[i,j] for i in (1:size(x)[1])[x[:] .& y], j in 1:size(x)[2]]

###############################################################################
## Logical Combo [x,y] two index

Base.getindex(x::LogicalCombo, row::Int, ::Colon; bool::Bool=false) =
  [ x[row,i] for i = !bool:size(x)[2] ]

Base.getindex(x::LogicalCombo, y::Union{Integer,Array}, col::String) = x[y, Symbol(col)]

Base.getindex(x::LogicalCombo, ::Colon, col::String) = x[:, Symbol(col)]

Base.getindex(x::LogicalCombo, y::Union{Integer,Array{Int64}}, ::Colon) =
  hcat([ x[i,:] for i in collect(y) ]...)

Base.getindex(x::LogicalCombo, row::UnitRange, ::Colon) = x[collect(row),:]

Base.getindex(x::LogicalCombo, ::Colon, col::Union{Int64,Symbol}; bool::Bool=false) =
  [ x[i,col] for i = 1:size(x)[1]]

# Base.getindex(x::LogicalCombo, ::Colon, col::Union{Int64,Symbol}) =
#   [ x[i,col] for i = 1:size(x)[1] ]

###############################################################################
## Logical Combo [x,y,z] three index
Base.getindex(x::LogicalCombo, ::Colon, ::Colon, ::Colon)   =
  [x[i,j] for i in (1:size(x)[1])[x[:]], j in 1:size(x)[2]]

Base.getindex(x::LogicalCombo, ::Colon, ::Colon, z::Bool)   =
  [x[i,j] for i in (1:size(x)[1])[x[:] .== z], j in 1:size(x)[2]]

Base.getindex(x::LogicalCombo, ::Colon, y::Union{Int64,Symbol,String}, ::Colon) =
  [x[i,y] for i in (1:size(x)[1])[x[:]]]

Base.getindex(x::LogicalCombo, ::Colon, y::Union{Int64,Symbol,String}, z::Bool) =
  [x[i,y] for i in (1:size(x)[1])[x[:] .== z]]

Base.getindex(x::LogicalCombo, y::Int64 , ::Colon, ::Colon) =
  x[(1:size(x)[1])[x[:]][y],:]

Base.getindex(x::LogicalCombo, y::Union{Int64, Array}, ::Colon, z::Bool) =
  x[(1:size(x)[1])[x[:] .== z][y],:]

Base.getindex(x::LogicalCombo, j::Int64, y::Union{Int64,Symbol,String}, z::Bool) =
  x[(1:size(x,1))[x[:] .== z][j],y]

###############################################################################
## Logical Combo [x,y;] two index with named arguments

Base.getindex(x::LogicalCombo, ::Colon, ::Colon; bool=false, varnames=false) =
  collect(x,bool=bool,varnames=varnames)

###############################################################################
# Set index!
# Base.setindex!(x::LogicalCombo, y::Union{Int64,UnitRange}) =  x.logical[y]

Base.setindex!(x::LogicalCombo, y::Bool, z::Integer)   = x.logical[z] = y
Base.setindex!(x::LogicalCombo, y::Bool, z::Union{UnitRange, AbstractArray}) =
  x.logical[z] .= y
Base.setindex!(x::LogicalCombo, y::Union{Array{Bool},Array{Bool,1},BitArray{1}}, ::Colon) = x.logical[:] .= y

Base.setindex!(x::LogicalCombo, y::Union{Array{Bool},Array{Bool,1}}, z::Union{UnitRange, AbstractArray}) =
  x.logical[z] = y

###############################################################################
# Extra Functions

Base.fill(v; each::Integer) = collect(Iterators.flatten([fill(x, each) for x in v]))

function Base.range(x::LogicalCombo; feasible=true)
  p = Dict()
  for i in 1:size(x,2); p[x.keys[i]] = sort(unique(x[:,i,feasible])); end
  p
end

Base.range(x::LogicalCombo, y::String; feasible=true) =
  Dict(Symbol(y) => x[:, Symbol(y), feasible] |> unique |> sort)

using StatsBase

function StatsBase.sample(x::LogicalCombo; n=1, feasible=true)
  if nfeasible(x) == 0
    return []
  elseif nfeasible(x) == 1
    return x[:,:,:]
  else
    y = x[sample(1:nfeasible(x, feasible), n, replace=false),:,feasible]
    (n == 1) && return y
    (n > 1) && return sort(y, dims=1)
 end
end

StatsBase.sample(x::LogicalCombo, n; feasible=true) = sample(x, n=n, feasible=feasible)
StatsBase.sample(x::LogicalCombo, n, feasible) = sample(x, n=n, feasible=feasible)


text:
function bisection(fun, a, b; threshold = 1e-10, maxiter = 100)
    root = (a+b)/2
    iter = 0
    while (abs(fun(root)) > threshold) && (iter < maxiter)
        if fun(a)*fun(root) < 0
            b = root
        else
            a = root
        end
        root = (a+b)/2
        iter += 1
    end
    if iter == maxiter
        println("Maximal number of iterations reached in bisection method.")
    end
    root
end

####################
# Newton iterations
#####################

"The abstract NonlinearSystem type represents a system of non-linear equations in n variables."
abstract type NonlinearSystem
end

# Default numeric type is Float64
# eltype(::Type{S}) where {S <: NonlinearSystem} = Float64

function residual(sys::NonlinearSystem, x)
    result = Array{eltype(sys)}(length(sys))
    residual!(result, sys, x)
end

function jacobian(sys::NonlinearSystem, x)
    J = Array{eltype(sys)}(size(sys))
    jacobian!(J, sys, x)
end

"""
Apply Newton iterations with the given starting values to find the roots of the given
system of equations. A convergence threshold and maximum number of iterations are optional.

The maximal step size in each iteration can be limited using `maxstep`.

Damped Newton iterations can be achieved by supplying damping parameter 0 < alpha < 1.
"""
function newton(sys::NonlinearSystem, x0; threshold = 1e-6, alpha = 1, maxiter = 100, maxstep = 1e5, verbose = false)
    n = length(x0)
    S = Array{eltype(sys)}(undef, n)
    J = Array{eltype(sys)}(undef, n, n)

    residual!(S, sys, x0)
    jacobian!(J, sys, x0)
    c = cond(J)
    threshold = max(threshold, eps(eltype(x0))*c)
    u = J \ S
    if norm(alpha*u) > maxstep
        u = u * maxstep/norm(u)
    end
    x1 = x0 - alpha*u

    iter = 1
    while norm(x1-x0) > threshold && iter < maxiter
        iter += 1
        x0 = x1
        residual!(S, sys, x0)
        jacobian!(J, sys, x0)
        u = J \ S
        if norm(alpha*u) > maxstep
            u = u * maxstep/norm(u)
        end
        x1 = x0 - alpha*u
    end
    if verbose
        if iter == maxiter
            println("Newton: maximum number of iterations, no convergence")
        else
            println("Newton: Converged to solution in ", iter, " iterations.")
        end
    end
    x1,iter,norm(x1-x0)
end

"Invoke Newton but restart with smaller maximal step size if an error is thrown."
function newton_with_restart(args...; maxstep = 1e5, maxiter = 100, verbose = false, options...)
    local nx,iter,normx
    nb_decreased = 0
    try
    nx,iter,normx = newton(args...; maxstep = maxstep, options...)
    catch e
        print("Error thrown by newton: ")
        println(e)
        iter = maxiter
    end
    while iter == maxiter
        verbose && println("Decreasing step size.")
        maxstep /= 2
        try
            nx,iter,normx = newton(args...; maxstep = maxstep, verbose=verbose, options...)
        catch e
            print("Error thrown by newton: ")
            println(e)
            iter = maxiter
        end
        nb_decreased += 1
        if nb_decreased > 20
            error("Too many step decreases in newton_with_restart. Aborting.")
        end
    end
    nx,iter,normx
end


text:
# test suite that can be used for all packages inheriting from GPUArrays
#
# use this by either calling `test`, or inspecting the dictionary `tests`

module TestSuite

export supported_eltypes

using GPUArrays

using LinearAlgebra
using Random
using Test

using FFTW
using Adapt
using FillArrays

struct ArrayAdaptor{AT} end
Adapt.adapt_storage(::ArrayAdaptor{AT}, xs::AbstractArray) where {AT} = AT(xs)

function compare(f, AT::Type{<:AbstractGPUArray}, xs...; kwargs...)
    # copy on the CPU, adapt on the GPU, but keep Ref's
    cpu_in = map(x -> isa(x, Base.RefValue) ? x[] : deepcopy(x), xs)
    gpu_in = map(x -> isa(x, Base.RefValue) ? x[] : adapt(ArrayAdaptor{AT}(), x), xs)

    cpu_out = f(cpu_in...; kwargs...)
    gpu_out = f(gpu_in...; kwargs...)

    collect(cpu_out) ≈ collect(gpu_out)
end

function supported_eltypes()
    (Float32, Float64, Int32, Int64, ComplexF32, ComplexF64)
end

# list of tests
const tests = Dict()
macro testsuite(name, ex)
    safe_name = lowercase(replace(name, " "=>"_"))
    fn = Symbol("test_$(safe_name)")
    quote
        $fn(AT) = $(esc(ex))(AT)

        @assert !haskey(tests, $name)
        tests[$name] = $fn
    end
end

include("testsuite/construction.jl")
include("testsuite/gpuinterface.jl")
include("testsuite/indexing.jl")
include("testsuite/io.jl")
include("testsuite/base.jl")
#include("testsuite/vector.jl")
include("testsuite/mapreduce.jl")
include("testsuite/broadcasting.jl")
include("testsuite/linalg.jl")
include("testsuite/math.jl")
include("testsuite/fft.jl")
include("testsuite/random.jl")
include("testsuite/uniformscaling.jl")

"""
Runs the entire GPUArrays test suite on array type `AT`
"""
function test(AT::Type{<:AbstractGPUArray})
    for (name, fun) in tests
        code = quote
            $fun($AT)
        end
        @eval @testset $name $code
    end
end

end


text:
m = [1+1im 2+2im 3+3im 4+4im;
     5+5im 6+6im 7+7im 8+8im;
     9+9im 10+10im 11+11im 12+12im;
     13+13im 14+14im 15+15im 16+16im]

qm = QuArray(m, LabelBasis(4,1), LabelBasis(2,2))
qm_qm = tensor(qm,qm)

@assert qm[(2,0),(1,0)] == qm[3,2]
@assert qm'[(1,0),(2,0)] == qm[3,2]'
@assert QuBase.bases(qm_qm) == (LabelBasis(4,1,4,1), LabelBasis(2,2,2,2))

@assert qm_qm[(2,0,1,0),(1,1,0,1)] == qm_qm[7,12]
qm_qm[(2,0,1,0),(1,1,0,1)] = 123
@assert qm_qm[(2,0,1,0),(1,1,0,1)] == 123

x2subspace = [filter((x...)->sum(x...)==2, QuBase.bases(qm_qm,2))...]
@assert x2subspace == [(1,1,0,0),(1,0,1,0),(0,1,1,0),(1,0,0,1),(0,1,0,1),(0,0,1,1)]

@assert qm_qm[(3,0,3,0), x2subspace] == [0+416im 0+392im 0+420im 0+420im 0+450im 0+416im]
qm_qm[(3,0,3,0), x2subspace] = [1 2 3 4 5 6]
@assert qm_qm[(3,0,3,0), x2subspace] == [1 2 3 4 5 6]


text:
"""
    dmrg_step!(state::SweepState, set::SweepSettings)

Update `state` by doing a two-site DMRG step with the parameters in `set`.
"""
function dmrg_step!(state::SweepState{L,T}, set::SweepSettings) where {L,T}
    psi, H, site = state.psi, state.H, state.site

    left = state.H_cntrctns[site-1].tnsr
    right = state.H_cntrctns[site+2].tnsr

    # The previous two-site wavefunction will be used as the starting vector:
    #
    #        m     n
    #        |     |
    #     a -#- b -#- c
    @tensor prev[a, m, n, c] :=
        psi.tnsrs[site][a, m, b] * psi.tnsrs[site+1][b, n, c]

    # Application of the two-site-projected Hamiltonian onto a vector:
    #
    #     #- c           i -#
    #     #     n     p     #
    #     #     |     |     #
    #     #- b -#- d -#- h -#
    #     #     |     |     #
    #     #     m     o     #
    #     #     |     |     #
    #     #- a -#######- g -#
    function H_f(dest::AbstractVector, src::AbstractVector)
        dest_ = reshape(dest, size(prev)...)
        src_ = reshape(src, size(prev)...)
        @tensor dest_[c, n, p, i] =
            left[a, b, c] * src_[a, m, o, g] *
            H.tnsrs[site][b, m, n, d] * H.tnsrs[site+1][d, o, p, h] *
            right[g, h, i]
        dest
    end
    H_map = LinearMap{T}(H_f, prod(size(prev)); issymmetric=true)

    # Improve our guess for the ground state of the two-site-projected
    # Hamiltonian by using a sparse diagonalizer. Since there's no way to
    # request a specific number of iterations, we just say we're converged
    # after a single iteration and repeat several times.
    F = eigs(H_map; nev=1, which=:SR, tol=Inf, maxiter=1, v0=vec(prev))
    wf = vec(F[2])
    for _ in 1:(set.num_iters-1)
        F = eigs(H_map; nev=1, which=:SR, tol=Inf, maxiter=1, v0=wf)
        wf = vec(F[2])
    end

    wf_mat = reshape(wf, size(prev, 1) * size(prev, 2),
                         size(prev, 3) * size(prev, 4))
    U, S, V = svd(wf_mat)
    trunc_len = length(S)
    trunc_cutoff = 0.0
    while trunc_len > set.bond_max ||
          (set.bond_min < trunc_len &&
           trunc_cutoff + S[trunc_len]^2 <= set.cutoff_max)

        trunc_cutoff += S[trunc_len]^2
        trunc_len -= 1
    end
    # Because the diagonalization results in a normalized eigenvector, we are
    # guaranteed that sum(S.^2) == 1. However, when we throw away some singular
    # values, we should renormalize.
    S ./= sqrt(sum(S[1:trunc_len].^2))
    # Reduced density matrix eigenvalues.
    eigvals = S[1:trunc_len].^2

    # Update the MPS.
    if state.dir == Right
        A = U[:, 1:trunc_len]
        B = diagm(S[1:trunc_len]) * V[:, 1:trunc_len]'
    elseif state.dir == Left
        A = U[:, 1:trunc_len] * diagm(S[1:trunc_len])
        B = V[:, 1:trunc_len]'
    end
    psi.tnsrs[site] = reshape(A, size(prev, 1), size(prev, 2), trunc_len)
    psi.tnsrs[site+1] = reshape(B, trunc_len, size(prev, 3), size(prev, 4))

    # Update the contractions.
    if state.dir == Right
        state.H_cntrctns[site] = state.H_cntrctns[site-1] * (psi, H)
    elseif state.dir == Left
        state.H_cntrctns[site+1] = state.H_cntrctns[site+2] * (psi, H)
    end
    state.H2_cntrctn = state.H2_cntrctn * (psi, H, H)

    # Update the maximum bond dimension.
    state.mps_max_bond_dim = maximum(size(tnsr, 3) for tnsr in psi.tnsrs)

    eigvals
end

"""
    dmrg!(psi::MPS{L,T}, H::MPO{L,T}, sch::SweepSchedule;
          outputs::Vector{<:SweepOutput}=SweepOutput[])

Perform two-site DMRG on the state `psi` using the Hamiltonian `H` with
parameters specified in `sch`.
"""
function dmrg!(psi::MPS{L,T}, H::MPO{L,T}, sch::SweepSchedule;
               outputs::Vector{<:SweepOutput}=SweepOutput[]) where {L,T}
    L >= 2 || throw(DomainError(L, "At least 2 sites."))

    state = SweepState(psi, H)
    sweep_details = SweepDetails[]

    init.(outputs)

    converged = false
    for n in 1:sch.max_sweeps
        sweep_start_time = time()

        set = sch[n]

        state.sweep_num = n
        if L == 2
            state.dir = Right
        else
            state.dir = flip(state.dir)
        end
        state.H2_cntrctn = cap_contraction(T, state.dir, L, 2)
        if L == 2
            # 1 2
            # ^ ^
            range = 1:1
        else
            if state.dir == Right
                # 1 2 3 4 5
                # ^ ^
                #   ^ ^
                #     ^ ^
                range = 1:(L-2)
            elseif state.dir == Left
                # 1 2 3 4 5
                #       ^ ^
                #     ^ ^
                #   ^ ^
                range = (L-1):-1:2
            end
        end

        # Sweep!
        for i in range
            state.site = i
            eigvals = dmrg_step!(state, set)
            i == div(L, 2) && (state.middle_eigvals = eigvals)

            step.(outputs, state)
        end

        if L == 2
            energy = state.H_cntrctns[1] * (psi, H) * state.H_cntrctns[L+1]
            H2 = state.H2_cntrctn * (psi, H, H) *
                 cap_contraction(T, Left, L, 2)
        else
            if state.dir == Right
                energy = state.H_cntrctns[L-2] * (psi, H) * state.H_cntrctns[L]
                H2 = state.H2_cntrctn * (psi, H, H) * (psi, H, H) *
                     cap_contraction(T, Left, L, 2)
            elseif state.dir == Left
                energy = state.H_cntrctns[1] * (psi, H) * state.H_cntrctns[3]
                H2 = cap_contraction(T, Right, L, 2) *
                     (psi, H, H) * (psi, H, H) * state.H2_cntrctn
            end
        end

        state.energy = realize(energy)
        state.dH2 = realize(H2)/energy^2 - 1.0
        state.duration = time() - sweep_start_time

        push!(sweep_details, SweepDetails(set, state))

        sweep.(outputs, state)

        if abs(state.dH2) <= sch.tolerance
            converged = true
            break
        end
    end

    done.(outputs)

    SweepHistory(sweep_details, converged)
end


text:
cfl = 0.25
dx = 0.05

T_per_Gi = 8 #how many folding times
L_per_l = 2

NxMesh = 20

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

u2gam(u) = sqrt(1+u^2)
u2v(u) = u/u2gam(u)

Gi(gam) = 1/sqrt(8*gam^2)

#--------------------------------------------------
function setup_twostream(u)

    gam = u2gam(u)
    v = u2v(u)

    println("u:   $u")
    println("v:   $v")

    println("gam: $gam")

    G = Gi(gam)
    println("G:   $G")

    dt = cfl*dx
    Nt = (T_per_Gi / G)/dt
    println("Nt:  $Nt")

    km = sqrt( (3/8) / (gam^3 * v^2) )
    lm = 2pi/km
    println("km:  $km")
    println("lm:  $lm")

    L = lm * L_per_l
    Nx = L/dx

    println("L:   $L")
    println("Nxf: $Nx")
    println("Nx:  $(Nx/NxMesh)")

    Nxsug = signif(Nx/NxMesh,1)
    println("suggested size ")
    println("Nx:$Nxsug NxMesh:$NxMesh dx:$(L/(NxMesh*Nxsug))")

    optimal_size(L, Nxsug, NxMesh)

    return
end

# Brute force search for best grid size
function optimal_size(L, Nx, NxMesh)

    err_min = Inf
    dxn_min = 0
    Nx2  = 0
    NxM2 = 0

    for NxM1 in range(Int(NxMesh) - 50, Int(NxMesh) + 50)
        (NxM1 < 5) && continue
        for Nx1 in range(Int(Nx) - 50, Int(Nx) + 50)
            (Nx1 < 5) && continue

            dxn = L/(NxM1*Nx1)
            err = abs(dxn - dx)/dx

            #println("NxM1: $NxM1 Nx1:$Nx1 dxn:$dxn err:$err")

            if err < err_min
                Nx2     = Nx1
                NxM2    = NxM1
                err_min = err
                dxn_min = dxn
                #println("better!")
            end
        end
    end

    println("optimal")
    println("Nx2:$Nx2 NxM2: $NxM2 dxn:$dxn_min err:$err_min")

end

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

println("--------------------------------------------------")
setup_twostream(1.0)

println("--------------------------------------------------")
setup_twostream(2.0)

println("--------------------------------------------------")
setup_twostream(3.0)

println("--------------------------------------------------")
setup_twostream(4.0)

max_stars_repo_name: Krissmedt/imprunko

max_issues_repo_name: Krissmedt/imprunko

max_forks_repo_name: Krissmedt/imprunko


text:
using DimensionalData

@inline gaussian(t, t₀, w) = exp(-((t - t₀) / w)^2)

const KDM = 4.148808e3 # MHz^2 pc^-1 cm^3 s

"""
Δt(DM, ν₁, ν₂)
Calculates the time delay corresponding to a dispersed pulse at dispersion measure `DM` in pc/cc between
frequencies `ν₁` and `ν₂`
"""
Δt(DM, ν₁, ν₂) = KDM * DM * (ν₁^-2 - ν₂^-2)

"""
fake_pulse(DM,f_low,f_high)

Generates a `Filterbank` file corresponding to a fake pulse of dispersion measure `DM`
from frequencies `f_low` to `f_high` (MHz).

# Optional Arguments
- `channels`: Number of frequency channels
- `samples`: Number of time samples
- 'δt`: The time step (s) represented by one sample
- `w`: Width of the pulse in time samples
- `A`: Amplitude of pulse in SNR
- `α`: Spectral index, i.e. attenuation of pulse over frequency
- `start`: Index of start of pulse
- `noise_floor`: Value of the top of the noise floor. This is 1 for `dtype` of Floats and one tenth `typemax` for integers by default.
- `dtype`: Data type of data
"""
function fake_pulse(DM, f_low, f_high;
                    channels=1024,
                    samples=1024,
                    δt=1e-3,
                    w=8,
                    A=2,
                    α=4,
                    start=1,
                    noise_floor::Union{T,Nothing}=nothing,
                    dtype::Type{T}=Float32) where {T}
    @assert start < samples "Starting sample must be less than the number of samples"

    if isnothing(noise_floor)
        if dtype <: Integer
            # Default to one tenth the range of the integer type
            noise_floor = typemax(T) ÷ 0xA
        else
            noise_floor = one(T)
        end
    end

    freqs = range(; start=f_high, stop=f_low, length=channels)
    time = range(; start=0, step=δt, length=samples)

    t_start = start * δt

    # Caclulate shifts
    shifts = @. Δt(DM, freqs', f_high)

    # Generate the raw pulse, these are floats at this point
    raw_pulse = @. gaussian(time, t_start + shifts, w * δt) * A * noise_floor * (freqs' / f_high)^α

    # Prep the dynamic spectrum with the background noise
    dyn_spec = rand(zero(T):noise_floor, samples, channels)

    # Add in the pulse
    if dtype <: Integer
        dyn_spec += round.(dtype,raw_pulse)
    else
        dyn_spec += raw_pulse
    end

    # Build Filterbank
    return Filterbank(DimArray(dyn_spec,
                               (Ti(time),
                                Freq(freqs))),
                      Dict("tsamp" => δt,
                           "nbits" => sizeof(dtype) * 8,
                           "nsamples" => samples,
                           "nchans" => channels,
                           "fch1" => f_high,
                           "foff" => step(freqs)))
end

export fake_pulse


text:
using Documenter, FlatBuffers

makedocs(
    modules = [FlatBuffers],
)

deploydocs(
    deps = Deps.pip("mkdocs", "mkdocs-material", "python-markdown-math"),
    repo = "github.com/dmbates/FlatBuffers.jl.git"
)


text:
# tensor product filtering for dimension D discrete signal.
function ranges2collection( x_ranges::Vector{LinRange{T,L}},
                            ::Val{D}) where {T,D,L}
    # set up.
    @assert !isempty(x_ranges)
    @assert length(x_ranges) == D
    N_array = collect( length(x_ranges[d]) for d = 1:D )
    N = prod(N_array)
    sz_N = tuple(N_array...)

    # Position.
    X_nD = Array{Vector{T},D}(undef,sz_N)
    for 𝑖 in CartesianIndices(sz_N)
        X_nD[𝑖] = Vector{T}(undef,D)

        for d = 1:D
            X_nD[𝑖][d] = x_ranges[d][𝑖[d]]
        end
    end

    return X_nD
end

"""
    convertcompactdomain(x::T, a::T, b::T, c::T, d::T)::T

converts compact domain x ∈ [a,b] to compact domain out ∈ [c,d].
"""
function convertcompactdomain(x::T, a::T, b::T, c::T, d::T)::T where T <: Real

    return (x-a)*(d-c)/(b-a)+c
end

function convertcompactdomain(x::Vector{T}, a::Vector{T}, b::Vector{T}, c::Vector{T}, d::Vector{T})::Vector{T} where T <: Real

    return collect( convertcompactdomain(x[i], a[i], b[i], c[i], d[i]) for i = 1:length(x) )
end

max_stars_repo_name: RoyCCWang/PatchMixtureKriging

max_issues_repo_name: RoyCCWang/PatchMixtureKriging


text:
# implementation/indices.jl
#
# Implements the index calculations, i.e. converting the tensor labels into
# indices specifying the operations

# Extract index information
#---------------------------
function add_indices(IA, IC)
    indCinA::Vector{Int}=indexin(collect(IC), collect(IA))
    isperm(indCinA) || throw(IndexError("invalid index specification: $IA to $IC"))
    return indCinA
end

function trace_indices(IA, IC)
    indCinA::Vector{Int}=indexin(collect(IC), collect(IA))

    cI=unique(setdiff(IA, IC))
    cindA1=Array(Int, length(cI))
    cindA2=Array(Int, length(cI))
    for i=1:length(cI)
        cindA1[i] = findfirst(IA, cI[i])
        cindA2[i] = findnext(IA, cI[i], cindA1[i]+1)
        findnext(IA, cI[i], cindA2[i]+1)==0 || throw(IndexError("invalid trace specification: $IA to $IC"))
    end
    pA = vcat(indCinA, cindA1, cindA2)
    isperm(pA) || throw(IndexError("invalid trace specification: $IA to $IC"))
    return indCinA, cindA1, cindA2
end

function contract_indices(IA, IB, IC)
    # Compute contraction indices and check for valid permutation
    NA = length(IA)
    NB = length(IB)
    NC = length(IC)

    NA == length(unique(IA)) || throw(IndexError("handle partial trace first: $IA"))
    NB == length(unique(IB)) || throw(IndexError("handle partial trace first: $IB"))
    NC == length(unique(IC)) || throw(IndexError("handle partial trace first: $IC"))

    cI = intersect(IA, IB)
    cN = length(cI)
    oIA = intersect(IC, IA)
    oNA = length(oIA)
    oIB = intersect(IC, IB)
    oNB = length(oIB)

    if cN+oNA != NA || cN+oNB != NB || oNA+oNB != NC
        throw(IndexError("invalid contraction pattern: $IA * $IB to $IC"))
    end

    cindA::Vector{Int} = indexin(cI, collect(IA))
    oindA::Vector{Int} = indexin(oIA, collect(IA))
    cindB::Vector{Int} = indexin(cI, collect(IB))
    oindB::Vector{Int} = indexin(oIB, collect(IB))
    indCinoAB::Vector{Int} = indexin(collect(IC), vcat(oIA, oIB))

    if !isperm(vcat(oindA, cindA)) || !isperm(vcat(oindB, cindB)) || !isperm(indCinoAB)
        throw(IndexError("invalid contraction pattern: $IA and $IB to $IC"))
    end

    return oindA, cindA, oindB, cindB, indCinoAB
end


text:
module Vikor

import GreyDecision.GreyNumbers.GreyNumber
import GreyDecision.Utility: unitize, apply_columns, weightize, reverseminmax

function vikor(decisionMat::Array{GreyNumber{T}, 2}, weights::Array{Float64,1}, fns::Array{Function,1};  v::Float64=0.5) where {T <: Real}

    w = unitize(weights)

    nalternatives, ncriteria = size(decisionMat)

    col_max = apply_columns(fns, decisionMat)
    col_min = apply_columns(reverseminmax(fns), decisionMat)

    A = similar(decisionMat)

    for i in 1:nalternatives
        for j in 1:ncriteria
            if fns[j] == maximum
                 A[i, j] = abs((col_max[j] - decisionMat[i, j]) / (col_max[j] - col_min[j]))
            elseif fns[j] == minimum
                 A[i, j] = abs((decisionMat[i, j] - col_min[j]) / (col_max[j] - col_min[j]))
            else
                @warn fns[j]
                error("Function must be either maximum or minimum.")
            end
        end
    end

    weightedA = weightize(A, w)

    s = Array{GreyNumber{Float64},1}(undef, nalternatives)
    r = similar(s)
    q = similar(s)

    for i in 1:nalternatives
        s[i] = sum(weightedA[i,:])
        r[i] = maximum(weightedA[i,:])
    end

    smin = minimum(s)
    smax = maximum(s)
    rmin = minimum(r)
    rmax = maximum(r)
    ranges = (smax - smin)
    ranger = (rmax - rmin)
    for i in 1:nalternatives
        #q = v .* ((s .- smin ./ (smax - smin))) .+ (1 - v) .* ((r .- rmin ./ (rmax - rmin)))
        q[i] = v * (s[i] - smin / ranges) + (1.0-v) * (r[i] - rmin / ranger)
    end

    best_index = sortperm(q) |> first

    return Dict(
        :scores => q,
        :bestindex => best_index
    )
end # end of Vikor function

end # End of module Vikor


text:
@testset "Environment" begin

    @info "Set, get and delete ENV keys"
    ENV["CONFIGS_FOO"] = "test"
    @test Configs.parseenvkey("CONFIGS_FOO", false) === "test"
    @test ENV["CONFIGS_FOO"] === "test"
    @test_nowarn Configs.deleteenvkey!("CONFIGS_FOO")
    @test_nowarn Configs.deleteenvkey!("CONFIGS_NOTTHERE")
    @test !haskey(ENV, "CONFIGS_FOO")
    @test Configs.parseenvkey("CONFIGS_FOO", false) === false

    @info "Initialises, sets and gets"

    @test_nowarn initconfig(; configs_directory = normpath(joinpath(pwd(), "configs")))
    @test_throws Configs.Configserror initconfig()
    @test_nowarn setconfig!("a.test.int", 100)
    @test_nowarn setconfig!("a.test.float", 100.00)
    @test_nowarn setconfig!("a.test.true", true)
    @test_nowarn setconfig!("a.test.false", false)
    @test_nowarn setconfig!("a.test.val", "test")
    @test_nowarn setconfig!("otherstuff.defaultmessage", "test")
    @test_nowarn setconfig!("otherstuff.empty", """{
        "test": {
            "newpath": "test",
            "newarray": ["one", 1]
        }
    }""")
    @test_nowarn setconfig!("otherstuff.tuples", (; foo = (1,2,3)))
    @test_nowarn setconfig!("otherstuff.tuples.bar", [1,2,3])
    @test_nowarn setconfig!("otherstuff.tuples.foobar", (1,2,3))
    @test_nowarn setconfig!("otherstuff.dict", Dict(:foo => "bar"))
    @test getconfig("a.test.val") === "test"
    @test getconfig("database.credentials.password") === "supersecret"
    @test getconfig("otherstuff.defaultmessage") === "test"
    @test getconfig("otherstuff.empty.test.newpath") === "test"
    @test getconfig("otherstuff.dict.foo") === "bar"
    @test isequal(getconfig("otherstuff.tuples.bar"), (1,2,3))
    @test isequal(getconfig("otherstuff.tuples.foo"), (1,2,3))
    @test isequal(getconfig("otherstuff.tuples.foobar"), (1,2,3))

    @info "Errors if no key"
    @test_throws Configs.Configserror getconfig("nogo.error")

    @info "Checks for hasconfig"
    @test !hasconfig("nogo.error")
    @test hasconfig("a.test.val")
end

@testset "imutability" begin
    @info "Cannot set after get"
    @test_throws Configs.Configserror setconfig!("nogo", "error")
    @info "The config is immutable"
    parselib(value) = (value isa Number || value isa String)
    function parselib(branch::NamedTuple)
        for value in branch
            parselib(value)
        end
    end
    function parselib(branch::Tuple)
        for value in branch
            parselib(value)
        end
    end
    parselib(Configs.state.immutable)
end

max_stars_repo_name: citkane/Config

max_issues_repo_name: citkane/Config

max_forks_repo_name: citkane/Config


text:
#####
#inteval censored
#####
function make_ic(rng::AbstractRNG, dt::UnivariateDistribution, N::Int)
    ER = zeros(N)
    EL = zeros(N)
    S = zeros(N)
    at = 0.0
    for i in 1:N
        ue = rand(rng)
        at -= log(rand(rng))
        tau = -log(rand(rng))
        y = rand(rng, dt)
        S[i] = at + y
        ud = rand(rng)
        EL[i] = max(0.0, at - tau * (1.0-ue))
        ER[i] = min(S[i], at + tau * ue)
    end
    return EL, ER, S
end

#inteval censored (infinite)
function make_ic(rng::AbstractRNG, dt::UnivariateDistribution, N::Int, p)
    ER = zeros(N)
    EL = zeros(N)
    S = zeros(N)
    at = 0.0
    for i in 1:N
        ue = rand(rng)
        at -= log(rand(rng))
        tau = -log(rand(rng))
        y = rand(rng,dt)
        S[i] = at + y
        ud = rand(rng)
        if ud <= p
            EL[i] = max(0.0, at - tau * (1.0-ue))
        else
            EL[i] = -Inf
        end
        ER[i] = min(S[i], at + tau * ue)
    end
    return EL, ER, S
end

#inteval censored with right truncated
function make_icrt(rng::AbstractRNG, dt::ContinuousUnivariateDistribution, N::Int, Tmax)
    L = zeros(N)
    R = zeros(N)
    S = zeros(N)
    d = trues(N)
    for i in 1:N
        ue = rand(rng)
        y = rand(rng,dt)
        at = rand(rng) * Tmax
        tau = -log(rand(rng))
        y = rand(rng,dt)
        S[i] = at + y
        L[i] = max(0.0, at - tau * (1.0-ue))
        R[i] = min(S[i], at + tau * ue)
        d[i] = (S[i] <= Tmax)
        d[i] = (S[i] <= Tmax)
    end
    return L[d], R[d], S[d]
end

#inteval censored (infinite) with right truncated
function make_icrt(rng::AbstractRNG, dt::ContinuousUnivariateDistribution, N::Int, Tmax, p)
    L = zeros(N)
    R = zeros(N)
    S = zeros(N)
    d = trues(N)
    for i in 1:N
        ue = rand(rng)
        y = rand(rng,dt)
        at = rand(rng) * Tmax
        tau = -log(rand(rng))
        ud = rand(rng)
        if ud <= p
            L[i] = max(0.0, at - tau * (1.0-ue))
        else
            L[i] = -Inf
        end
        R[i] = min(S[i], at + tau * ue)
        d[i] = (S[i] <= Tmax)
    end
    return L[d], R[d], S[d]
end

#####
#doubly inteval censored
#####
function make_dic(rng::AbstractRNG, dt::ContinuousUnivariateDistribution, N::Int)
    EL = zeros(N)
    ER = zeros(N)
    SL = zeros(N)
    SR = zeros(N)
    at = 0.0
    for i in 1:N
        ue = rand(rng)
        y = rand(rng, dt)
        at -= log(rand(rng))
        S = y + at
        ue = rand(rng)
        us = rand(rng)
        tau_e = -log(rand(rng))
        tau_s = -log(rand(rng))
        EL[i] = max(0.0, at - tau_e * (1.0 - ue))
        SR[i] = S + tau_s * us
        SL[i] = max(EL[i], S - tau_s * (1.0 - us))
        ER[i] = min(SR[i], at + tau_e * ue)
    end
    return EL, ER, SL, SR
end

#infinite
function make_dic(rng::AbstractRNG, dt::ContinuousUnivariateDistribution, N::Int, p)
    EL = zeros(N)
    ER = zeros(N)
    SL = zeros(N)
    SR = zeros(N)
    at = 0.0
    for i in 1:N
        at -= log(rand(rng))
        y = rand(rng,dt)
        S = y + at
        tau_e = -log(rand(rng))
        tau_s = -log(rand(rng))
        ue = rand(rng)
        ud = rand(rng)
        if ud <= p
            EL[i] = max(0.0, at - tau_e * (1.0 - ue))
        else
            EL[i] = -Inf
        end
        us = rand(rng)
        SR[i] = S + tau_s * us
        SL[i] = max(EL[i], S - tau_s * (1.0 - us))
        ER[i] = min(SR[i], at + tau_e * ue)
    end
    return EL, ER, SL, SR
end

#right truncated
function make_dicrt(rng::AbstractRNG, dt::ContinuousUnivariateDistribution, N::Int, Tmax)
    EL = zeros(N)
    ER = zeros(N)
    SL = zeros(N)
    SR = zeros(N)
    d = trues(N)
    for i in 1:N
        ue = rand(rng)
        y = rand(rng,dt)
        at = rand(rng) * Tmax
        S = y + at
        tau_e = -log(rand(rng))
        tau_s = -log(rand(rng))
        ue = rand(rng)
        ud = rand(rng)
        if ud <= p
            EL[i] = max(0.0, at - tau_e * (1d - ue))
        else
            EL[i] = -Inf
        end
        us = rand(rng)
        SR[i] = S + tau_s * us
        SL[i] = max(EL[i], S - tau_s * (1 - us))
        ER[i] = min(SR[i], at + tau_e * ue)
        d[i] = SR <= Tmax
    end
    return EL[d], ER[d], SL[d], SR[d]
end

#doubly inteval censored (infinite) with right truncated
function make_dicrt(rng::AbstractRNG, dt::ContinuousUnivariateDistribution, N::Int, Tmax, p)
    EL = zeros(N)
    ER = zeros(N)
    SL = zeros(N)
    SR = zeros(N)
    d = trues(N)
    for i in 1:N
        ue = rand(rng)
        y = rand(rng,dt)
        at = rand(rng) * Tmax
        S = y + at
        tau_e = -log(rand(rng))
        tau_s = -log(rand(rng))
        ue = rand(rng)
        ud = rand(rng)
        if ud <= p
            EL[i] = max(0.0, at - tau_e * (1.0 - ue))
        else
            EL[i] = -Inf
        end
        us = rand(rng)
        SR[i] = S + tau_s * us
        SL[i] = max(EL[i], S - tau_s * (1.0 - us))
        ER[i] = min(SR[i], at + tau_e * ue)
        d[i] = SR <= Tmax
    end
    return EL[d], ER[d], SL[d], SR[d]
end

max_stars_repo_name: abikoushi/IntervalCensored

max_issues_repo_name: abikoushi/IntervalCensored

max_forks_repo_name: abikoushi/IntervalCensored


text:
using Oceananigans.Operators: Δzᵃᵃᶜ

"""
    HorizontallyCurvilinearAnisotropicDiffusivity{N, K}

Holds viscosity and diffusivities for models with prescribed isotropic diffusivities.
"""
struct HorizontallyCurvilinearAnisotropicDiffusivity{NH, NZ, KH, KZ} <: AbstractTurbulenceClosure
    νh :: NH
    νz :: NZ
    κh :: KH
    κz :: KZ
end

const HCAD = HorizontallyCurvilinearAnisotropicDiffusivity

"""
    HorizontallyCurvilinearAnisotropicDiffusivity(; ν=0, κ=0)

Returns parameters for an horizontal diffusivity model on
curvilinear grids with viscosity `ν` and diffusivities `κ` for each tracer
field in `tracers`. `ν` and the fields of `κ` may be constants or functions
of `(x, y, z, t)`, and may represent molecular diffusivities in cases that all flow
features are explicitly resovled, or turbulent eddy diffusivities that model the effect of
unresolved, subgrid-scale turbulence.
`κ` may be a `NamedTuple` with fields corresponding
to each tracer, or a single number to be a applied to all tracers.
"""
function HorizontallyCurvilinearAnisotropicDiffusivity(FT=Float64; νh=0, κh=0, νz=0, κz=0)
    νh = convert_diffusivity(FT, νh)
    νz = convert_diffusivity(FT, νz)
    κh = convert_diffusivity(FT, κh)
    κz = convert_diffusivity(FT, κz)
    return HorizontallyCurvilinearAnisotropicDiffusivity(νh, νz, κh, κz)
end

function with_tracers(tracers, closure::HorizontallyCurvilinearAnisotropicDiffusivity)
    κh = tracer_diffusivities(tracers, closure.κh)
    κz = tracer_diffusivities(tracers, closure.κz)
    return HorizontallyCurvilinearAnisotropicDiffusivity(closure.νh, closure.νz, κh, κz)
end

calculate_diffusivities!(K, arch, grid, closure::HorizontallyCurvilinearAnisotropicDiffusivity, args...) = nothing

viscous_flux_ux(i, j, k, grid, clock, closure::HCAD, U, args...) = + ν_δᶜᶜᶜ(i, j, k, grid, clock, closure.νh, U.u, U.v)
viscous_flux_uy(i, j, k, grid, clock, closure::HCAD, U, args...) = - ν_ζᶠᶠᶜ(i, j, k, grid, clock, closure.νh, U.u, U.v)
viscous_flux_uz(i, j, k, grid, clock, closure::HCAD, U, args...) = + ν_uzᶠᶜᶠ(i, j, k, grid, clock, closure.νz, U.u)

viscous_flux_vx(i, j, k, grid, clock, closure::HCAD, U, args...) = ν_ζᶠᶠᶜ(i, j, k, grid, clock, closure.νh, U.u, U.v)
viscous_flux_vy(i, j, k, grid, clock, closure::HCAD, U, args...) = ν_δᶜᶜᶜ(i, j, k, grid, clock, closure.νh, U.u, U.v)
viscous_flux_vz(i, j, k, grid, clock, closure::HCAD, U, args...) = ν_vzᶜᶠᶠ(i, j, k, grid, clock, closure.νz, U.v)

@inline function diffusive_flux_x(i, j, k, grid, clock, closure::HCAD, c, ::Val{tracer_index}, args...) where tracer_index
    @inbounds κh = closure.κh[tracer_index]
    return diffusive_flux_x(i, j, k, grid, clock, κh, c)
end

@inline function diffusive_flux_y(i, j, k, grid, clock, closure::HCAD, c, ::Val{tracer_index}, args...) where tracer_index
    @inbounds κh = closure.κh[tracer_index]
    return diffusive_flux_y(i, j, k, grid, clock, κh, c)
end

@inline function diffusive_flux_z(i, j, k, grid, clock, closure::HCAD, c, ::Val{tracer_index}, args...) where tracer_index
    @inbounds κz = closure.κz[tracer_index]
    return diffusive_flux_z(i, j, k, grid, clock, κz, c)
end

Base.show(io::IO, closure::HorizontallyCurvilinearAnisotropicDiffusivity) =
    print(io, "AnisotropicDiffusivity: " *
              "(νh=$(closure.νh), νz=$(closure.νz)), " *
              "(κh=$(closure.κh), κz=$(closure.κz))")


text:
# Basic utilities for meshing parametric surfaces. Produces a `GenericMesh`

"""
    meshgen(Ω::Domain,num_elements)

Generate a `GenericMesh` for the domain `Ω` with `num_elements` per entity. To
specify a different number of elements per entity, `num_elements` should be a
vector with as many elements as there are entities in `Ω`.

Requires the entities forming `Ω` to `ParametricEntity`.
"""
function meshgen(Ω::Domain,num_elements::Vector)
    # extract the ambient dimension for these entities (i.e. are we in 2d or
    # 3d). Only makes sense if all entities have the same ambient dimension.
    N  = ambient_dimension(first(Ω))
    @assert all(p->ambient_dimension(p)==N,entities(Ω))
    mesh = GenericMesh{N,Float64}()
    meshgen!(mesh,Ω,num_elements) # fill in
    return mesh
end
function meshgen(Ω::Domain,num_elements::Union{Int,Tuple{Int},Tuple{Int,Int}})
    n = length(entities(Ω))
    meshgen(Ω,[num_elements for _ in 1:n])
end

"""
    meshgen!(mesh,Ω,sz)

Similar to [`meshgen`](@ref), but append entries to `mesh`.
"""
function meshgen!(mesh::GenericMesh,Ω::Domain,num_elements)
    @assert length(entities(Ω)) == length(num_elements)
    for (ent,sz) in zip(Ω,num_elements)
        ent isa ParametricEntity || error("meshgen! only works on parametric entites")
        _meshgen!(mesh,ent,sz)
    end
    return mesh
end

function _meshgen!(mesh::GenericMesh,ent::ParametricEntity,sz)
    # extract relevant fields and mesh the entity
    f = parametrization(ent)
    d = domain(ent)
    els    = _meshgen(f,d,sz)
    # push related information to mesh
    E      = eltype(els)
    vals   = get!(mesh.elements,E,Vector{E}())
    istart = length(vals) + 1
    append!(vals,els)
    iend   = length(vals)
    haskey(mesh.ent2tags,ent) && @debug "skipping entity $(key(ent)): already present in mesh"
    mesh.ent2tags[ent] = Dict(E=>collect(istart:iend)) # add key
    return mesh
end

"""
    _meshgen(f,d,sz)

Create a UniformCartesianMesh` of `d` push-forward map. The cartesian mesh has size `sz`, and is uniform in parameter
space.

!!! warning
    For parametrizations with large `jacobian`, the mesh size in physical space can
    be far from uniform.
"""
function _meshgen(f,d,sz)
    grid = UniformCartesianMesh(d,sz)
    iter = ElementIterator(grid)
    els = [ParametricElement(f,d) for d in iter]
    return els
end

# Element iterator interface
function Base.size(iter::ElementIterator{<:ParametricElement,<:GenericMesh})
    E = eltype(iter)
    M = mesh(iter)
    els::Vector{E} = M.elements[E]
    return (length(els),)
end

function Base.getindex(iter::ElementIterator{<:ParametricElement,<:GenericMesh},i::Int)
    E = eltype(iter)
    M = mesh(iter)
    els::Vector{E} = M.elements[E]
    return els[i]
end

function Base.iterate(iter::ElementIterator{<:ParametricElement,<:GenericMesh}, state=1)
    state > length(iter) && (return nothing)
    iter[state], state + 1
end

max_stars_repo_name: WaveProp/ParametricSurfaces

max_issues_repo_name: WaveProp/ParametricSurfaces

max_forks_repo_name: WaveProp/ParametricSurfaces


text:
############################################################################
### bootstrap based standard errors
############################################################################

function se_bootstrap(semfit::SemFit; n_boot = 3000, data = nothing, kwargs...)

    if isnothing(data)
        data = get_data(observed(model(semfit)))
    end

    data = prepare_data_bootstrap(data)

    start = solution(semfit)

    new_solution = zero(start)
    sum = zero(start)
    squared_sum = zero(start)

    n_failed = 0.0

    converged = true

    for _ = 1:n_boot

        sample_data = bootstrap_sample(data)
        new_model = swap_observed(model(semfit); data = sample_data, kwargs...)

        new_solution .= 0.0

        try
            new_solution = get_solution(sem_fit(new_model; start_val = start))
        catch
            n_failed += 1
        end

        @. sum += new_solution
        @. squared_sum += new_solution^2

        converged = true
    end

    n_conv = n_boot - n_failed
    sd = sqrt.(squared_sum/n_conv - (sum/n_conv).^2)
    print(n_failed, "\n")
    return sd

end

function prepare_data_bootstrap(data)
    return Matrix(data)
end

function bootstrap_sample(data)
    nobs = size(data, 1)
    index_new = rand(1:nobs, nobs)
    data_new = data[index_new, :]
    return data_new
end


text:
"""
     IndexedMarkovChain{T}

Representation of a Markov Chain type with "arbitrary" support (`T`). Assumes
that **rows** add up to one.

"""
abstract type IndexedMarkovChain{T} <: AbstractArray{T, 1} end
Base.getindex(mc::IndexedMarkovChain, i)=supp(mc)[i]
Base.length(mc::IndexedMarkovChain)=length(supp(mc))
Base.size(mc::IndexedMarkovChain)=(length(mc), )
Base.IndexStyle(::Type{<:IndexedMarkovChain}) = IndexLinear()
supp(mc::IndexedMarkovChain, i)=error("Method `supp` must be implemented for $(typeof(MC)) ")
transition(mc::IndexedMarkovChain, i, j)=transition(mc)[i,j]

" Unconditional mean of a MarkovChain "
Statistics.mean(mc::IndexedMarkovChain)=dot(markov_invariant(mc), supp(mc))
Statistics.mean(mc::IndexedMarkovChain, v::AbstractVector)=dot(markov_invariant(mc), v)

"Markov Chain with pre-stored cumulative transition matrix"
struct MarkovChain{T, S} <: IndexedMarkovChain{T}
    support::Vector{T}
    transition::Matrix{S}
    transition_cumu::Matrix{S}
end
supp(mc::MarkovChain)=mc.support
transition(MC::MarkovChain)=MC.transition # Make rows add up to one
transition_cumu(MC::MarkovChain)=MC.transition_cumu
Base.eltype(mc::MarkovChain{T, S}) where {T, S}=T

" Simples constructor for any iterable `v`. Note: no sanity check for dimensions."
MarkovChain(v, M::Matrix)=MarkovChain([v[i] for i in 1:length(v)], M,
                                      cumsum(M, dims=2))

" Same as `MarkovChain`, but pre-store invariant distribution "
struct MarkovChainExtended{T, S} <: IndexedMarkovChain{T}
    mc::MarkovChain{T, S}
    invariant::Array{Float64, 1}
end

MarkovChainExtended(mc::MarkovChain{T, S}) where {T,S}=
    MarkovChainExtended{T, S}(mc, markov_invariant(mc))
MarkovChainExtended(v, M::Matrix, tol=1e-12)=
    MarkovChainExtended(MarkovChain(v, M))

supp(mc::MarkovChainExtended)=supp(mc.mc)
transition(mc::MarkovChainExtended)=transition(mc.mc)
transition_cumu(mc::MarkovChainExtended)=transition_cumu(mc.mc)

# Note: will also work with `Aggregate` (& anything for which `transition` is implemented)
function expect_markov(v::AbstractArray, current_state, MC)
    e = 0.

    for next_state in 1:length(MC)
        e += v[next_state] * transition(MC, current_state, next_state)
    end

    e
end

" TODO: unrolled version for arbitrary tuple of arrays "
function expect_markov(v1::AbstractArray, v2::AbstractArray, current_state, MC)
    e1 = e2 = 0.

    for next_state in 1:length(MC)
        e1 += v1[next_state] * transition(MC, current_state, next_state)
        e2 += v2[next_state] * transition(MC, current_state, next_state)
    end

    e1, e2
end

markov_invariant(mc::MarkovChainExtended)=mc.invariant
markov_invariant(mc::MarkovChain, tol=1e-12)=markov_invariant(mc.transition,tol)

" Returns eigenvector associated with A's last eigenvalue "
markov_invariant(A, tol)=begin
    ed = eigen(transpose(A))

    !within_tolerance(ed.values[end], 1, tol) &&
        error("No invariant distribution found")

    evec = Float64.(ed.vectors[:, end])

    makeprob!(evec)

    return evec
end

makeprob!(u)=(u ./= sum(u))

"""

Simulate a Markov Chain  for `len` periods, starting with index `i0`

The `rand_tmp` vector should contain `len` draws from a uniform distribution.

Overwrites vector `path_prealloc`.

"""
function markov_genpath!(path_prealloc, mc, i0, len, rand_tmp)
    for t in 1:len
        i0 = draw_next(mc, i0, rand_tmp[t])
        path_prealloc[t] = i0
    end
end

function markov_genpath_seeded(mc, i0, len, seed=1234)
    rand_tmp=rand(len)
    path_prealloc = fill(0, len)
    markov_genpath!(path_prealloc, mc, i0, len, rand_tmp)
    return path_prealloc
end

function markov_genpath(mc, i0, len)
    rand_tmp=rand(len)
    path_prealloc = fill(0, len)
    markov_genpath!(path_prealloc, mc, i0, len, rand_tmp)
    return path_prealloc
end

" Guaranteed to not be `nothing`; use with EXTREME caution "
function findfirstg(bvec)
    idx = 1

    for i in 1:length(bvec)
        if bvec[i]
            idx = i
            break
        end
    end

    return idx
end

" Guaranteed to not be `nothing`; use with EXTREME caution "
function findfirstpositive(nvec)
    j = 0
    while (nvec[j+1] < 0) && (j+1 < length(nvec))
        j += 1
    end
    j+1
end

draw_next(mc::IndexedMarkovChain, i0, r)=begin
    @views findfirstpositive(transition_cumu(mc)[i0, :] .- r)
end

draw_next(mc::IndexedMarkovChain, i0)=draw_next(mc, i0, rand())

indeptensor(mc0::IndexedMarkovChain, mc1::IndexedMarkovChain)=MarkovChainExtended(vec(collect(Base.product(mc0, mc1))),
                                                                                Base.kron(transition(mc1), transition(mc0)))

function makemean1!(mc::IndexedMarkovChain{T}) where {T <: Number}
    mm = Statistics.mean(mc)
    foreach(1:length(mc)) do i
        supp(mc)[i] /= mm
    end
end


text:
module MetadataSyntax

"No meta-syntax."
no_meta_syntax = ()

"\!!var(one_backslash_escape:One backslash is skipped and treated like a normal meta syntax.)"
one_backslash_escape = ()

"!! var(space_between_backslash_metaname:Space between double ! and metaname is not a metasyntax.)"
space_between_backslash_metaname = ()

"!!var (space_between_metaname_bracket:Space between metaname and bracket is not a metasyntax.)"
space_between_metaname_bracket = ()

"!!var(license:[MIT](https://github.com/MichaelHatherly/Lexicon.jl/blob/master/LICENSE.md))"
brackets_within_meta = ()

"!!var(笔者:所以不多说了)"
chinese_unicode = ()

"\\!!var(russian:бежал мета)"
backslash_escaped_meta = ()

"\\!!var(unicode_meta_in_meta:所以不多说了 !!var(笔者_inner:бежал мета))"
backslash_escaped_nested_meta = ()

end


text:
## temporary place for all the stuff that should go to Yota.jl

function Base.reshape(dy::ZeroTangent, ::Tuple)
    return dy
end

# import ChainRulesCore

# function ChainRulesCore.rrule(::Colon, a::Int, b::Int)
#     y = a:b
#     function colon_pullback(dy)
#         return NoTangent(), NoTangent(), NoTangent()
#     end
#     return y, colon_pullback
# end

# test_rrule(Colon(), 1, 2; ouput_tangent=o)


text:
using StatsPlots
using CSV

jt_s = CSV.read(joinpath(@__DIR__, "julia_circle_apertures_nt=1.csv"))
jt_4 = CSV.read(joinpath(@__DIR__, "julia_circle_apertures_nt=4.csv"))
jt = CSV.read(joinpath(@__DIR__, "julia_circle_apertures_nt=8.csv"))
pt = CSV.read(joinpath(@__DIR__, "python_circle_apertures.csv"))

plot(markerstrokealpha=0, xscale=:log10, legend=:topleft)

@df pt plot!(:N, :time, c=2, label="photutils")
@df jt_s plot!(:N, :time, markerstrokealpha=0, c=1, label="Photometry.jl - 1 thread", yscale=:log10)
@df jt_4 plot!(:N, :time, c=3,  label="Photometry.jl - 4 threads")
@df jt plot!(:N, :time, c=3, ls=:dash, label="Photometry.jl - 8 threads")
title!("size(512, 512) - CircularAperture - exact")
ylabel!("time [s]")
xlabel!("number of apertures")

savefig(joinpath(@__DIR__, "..", "..", "docs", "src", "assets", "circle_apertures_benchmark.png"))


text:
using Flux
using Base: @get!
using MacroTools: @forward

const ϵ = 1e-8

# TODO: should use weak refs

"""
    Descent(η)

Classic gradient descent optimiser with learning rate `η`.
For each parameter `p` and its gradient `δp`, this runs `p -= η*δp`.
"""
mutable struct Descent
  eta::Float64
end

Descent() = Descent(0.1)

function apply!(o::Descent, x, Δ)
  Δ .*= o.eta
end

"""
    Momentum(params, η = 0.01; ρ = 0.9)

Gradient descent with learning rate `η` and momentum `ρ`.
"""
mutable struct Momentum
  eta::Float64
  rho::Float64
  velocity::IdDict
end

Momentum(η = 0.01, ρ = 0.9) = Momentum(η, ρ, IdDict())

function apply!(o::Momentum, x, Δ)
  η, ρ = o.eta, o.rho
  v = get!(o.velocity, x, zero(x))::typeof(data(x))
  @. v = ρ * v - η * Δ
  @. Δ = -v
end

"""
    Nesterov(eta, ρ = 0.9)

Gradient descent with learning rate  `η` and Nesterov momentum `ρ`.
"""
mutable struct Nesterov
  eta::Float64
  rho::Float64
  velocity::IdDict
end

Nesterov(η = 0.001, ρ = 0.9) = Nesterov(η, ρ, IdDict())

function apply!(o::Nesterov, x, Δ)
  η, ρ = o.eta, o.rho
  v = get!(o.velocity, x, zero(x))::typeof(data(x))
  d = @. ρ^2 * v - (1+ρ) * η * Δ
  @. v = ρ*v - η*Δ
  @. Δ = -d
end

"""
    RMSProp(η = 0.001, ρ = 0.9)

[RMSProp](https://www.cs.toronto.edu/~tijmen/csc321/slides/lecture_slides_lec6.pdf)
optimiser. Parameters other than learning rate don't need tuning. Often a good
choice for recurrent networks.
"""
mutable struct RMSProp
  eta::Float64
  rho::Float64
  acc::IdDict
end

RMSProp(η = 0.001, ρ = 0.9) = RMSProp(η, ρ, IdDict())

function apply!(o::RMSProp, x, Δ)
  η, ρ = o.eta, o.rho
  acc = get!(o.acc, x, zero(x))::typeof(data(x))
  @. acc = ρ * acc + (1 - ρ) * Δ^2
  @. Δ *= η / (√acc + ϵ)
end

"""
    ADAM(η = 0.001, β = (0.9, 0.999))

[ADAM](https://arxiv.org/abs/1412.6980v8) optimiser.
"""
mutable struct ADAM
  eta::Float64
  beta::Tuple{Float64,Float64}
  state::IdDict
end

ADAM(η = 0.001, β = (0.9, 0.999)) = ADAM(η, β, IdDict())

function apply!(o::ADAM, x, Δ)
  η, β = o.eta, o.beta
  mt, vt, βp = get!(o.state, x, (zero(x), zero(x), β))
  @. mt = β[1] * mt + (1 - β[1]) * Δ
  @. vt = β[2] * vt + (1 - β[2]) * Δ^2
  @. Δ =  mt / (1 - βp[1]) / (√(vt / (1 - βp[2])) + ϵ) * η
  o.state[x] = (mt, vt, βp .* β)
  return Δ
end

"""
    AdaMax(params, η = 0.001; β1 = 0.9, β2 = 0.999, ϵ = 1e-08)

[AdaMax](https://arxiv.org/abs/1412.6980v9) optimiser. Variant of ADAM based on
the ∞-norm.
"""
mutable struct AdaMax
  eta::Float64
  beta::Tuple{Float64,Float64}
  state::IdDict
end

AdaMax(η = 0.001, β = (0.9, 0.999)) = AdaMax(η, β, IdDict())

function apply!(o::AdaMax, x, Δ)
  η, β = o.eta, o.beta
  mt, ut, βp = get!(o.state, x, (zero(x), zero(x), β))
  @. mt = β[1] * mt + (1 - β[1]) * Δ
  @. ut = max(β[2] * ut, abs(Δ))
  @. Δ = (η/(1 - βp[1])) * mt/(ut + ϵ)
  o.state[x] = (mt, ut, βp .* β)
  return Δ
end

"""
    ADAGrad(η = 0.1; ϵ = 1e-8)

[ADAGrad](http://www.jmlr.org/papers/volume12/duchi11a/duchi11a.pdf) optimiser.
Parameters don't need tuning.
"""
mutable struct ADAGrad
  eta::Float64
  acc::IdDict
end

ADAGrad(η = 0.1) = ADAGrad(η, IdDict())

function apply!(o::ADAGrad, x, Δ)
  η = o.eta
  acc = get!(o.acc, x, fill(ϵ, size(x)))::typeof(data(x))
  @. acc += Δ^2
  @. Δ *= η / (√acc + ϵ)
end

"""
    ADADelta(ρ = 0.9, ϵ = 1e-8)

[ADADelta](https://arxiv.org/abs/1212.5701) optimiser. Parameters don't need
tuning.
"""
mutable struct ADADelta
  rho::Float64
  state::IdDict
end

ADADelta(ρ = 0.9) = ADADelta(ρ, IdDict())

function apply!(o::ADADelta, x, Δ)
  ρ = o.rho
  acc, Δacc = get!(o.state, x, (zero(x), zero(x)))
  @. acc = ρ * acc + (1 - ρ) * Δ^2
  @. Δ *= √Δacc/ (√acc + ϵ)
  @. Δacc = ρ * Δacc + (1 - ρ) * Δ^2
  return Δ
end

"""
    AMSGrad(η = 0.001, β = (0.9, 0.999))

[AMSGrad](https://openreview.net/forum?id=ryQu7f-RZ) optimiser. Parameters don't need
tuning.
"""
mutable struct AMSGrad
  eta::Float64
  beta::Tuple{Float64, Float64}
  state::IdDict
end

AMSGrad(η = 0.001, β = (0.9, 0.999)) = AMSGrad(η, β, IdDict())

function apply!(o::AMSGrad, x, Δ)
  η, β = o.eta, o.beta
  mt, vt, v̂t = get!(o.state, x, (fill(ϵ, size(x)), fill(ϵ, size(x)), fill(ϵ, size(x))))
  @. mt = β[1] * mt + (1 - β[1]) * Δ
  @. vt = β[2] * vt + (1 - β[2]) * Δ ^ 2
  @. v̂t = max.(v̂t, vt)
  @. Δ = η * mt / (√v̂t + ϵ)
end

"""
    NADAM(η = 0.001, β = (0.9, 0.999))

[NADAM](http://cs229.stanford.edu/proj2015/054_report.pdf) optimiser. Parameters don't need
tuning.
"""
mutable struct NADAM
  eta::Float64
  beta::Tuple{Float64, Float64}
  state::IdDict
end

NADAM(η = 0.001, β = (0.9, 0.999)) = NADAM(η, β, IdDict())

function apply!(o::NADAM, x, Δ)
  η, β = o.eta, o.beta
  β1p, β2p = o.beta
  mt, vt = get!(o.state, x, (zero(x), zero(x)))
  @. mt = β[1] * mt + (1 - β[1]) * Δ
  @. vt = β[2] * vt + (1 - β[2]) * Δ^2
  @. Δ = (β[1] * mt / (1 - β[1] * β1p) + (1 - β[1]) * Δ / (1 - β1p)) / (√(vt * β[2] / (1 - β2p)) + ϵ) * η
  o.state[x] = (mt, vt, (β1p * β[1], β2p * β[2]))
  return Δ
end

"""
    ADAMW((η = 0.001, β = (0.9, 0.999), decay = 0)

[ADAMW](https://arxiv.org/abs/1711.05101) fixing weight decay regularization in Adam.
"""
ADAMW(η = 0.001, β = (0.9, 0.999), decay = 0) =
  Optimiser(ADAM(η, β), WeightDecay(decay))

# Compose optimizers

"""
    Optimiser(a, b, c...)

Combine several optimisers into one; each optimiser produces a modified gradient
that will be fed into the next, and this is finally applied to the parameter as
usual.
"""
mutable struct Optimiser
  os::Vector{Any}
end

Optimiser(o...) = Optimiser(Any[o...])

@forward Optimiser.os Base.getindex, Base.first, Base.last, Base.lastindex, Base.push!, Base.setindex!
@forward Optimiser.os Base.iterate

Base.getindex(c::Optimiser, i::AbstractArray) = Optimiser(c.os[i]...)

function apply!(o::Optimiser, x, Δ)
  for opt in o.os
    Δ = apply!(opt, x, Δ)
  end
  return Δ
end

"""
`InvDecay(γ)`

Apply inverse time decay to an optimiser
```julia
  Optimiser(InvDecay(..), Opt(..))
```
"""
mutable struct InvDecay
  gamma::Float64
  state::IdDict
end

InvDecay(γ = 0.001) = InvDecay(γ, IdDict())

function apply!(o::InvDecay, x, Δ)
  γ = o.gamma
  n = get!(o.state, x, 1)
  Δ .*= 1 / (1 + γ * n)
  o.state[x] = n + 1
  return Δ
end

"""
`ExpDecay(eta, decay, decay_step, clip)`

Schedule the learning rate `eta` by `decay` every `decay_step` till a minimum of `clip`.

To apply exponential decay to an optimiser:
```julia
  Optimiser(ExpDecay(..), Opt(..))
```
"""
mutable struct ExpDecay
  eta::Float64
  decay::Float64
  step::Int64
  clip::Float64
  current::IdDict
end

ExpDecay(opt = 0.001, decay = 0.1, decay_step = 1000, clip = 1e-4) = ExpDecay(opt, decay, decay_step, clip, IdDict())

function apply!(o::ExpDecay, x, Δ)
  η, s, decay = o.eta, o.step, o.decay
  n = o.current[x] = get(o.current, x, 0) + 1
  if o.current[x]%s == 0 && count(x -> x%s == 0, values(o.current)) == 1
    η = max(η * decay^(s / n), o.clip)
    o.eta = η
  end
  @. Δ *= η
end

"""
`WeightDecay(wd)`

Decay the weight parameter by `wd`
"""
mutable struct WeightDecay
  wd::Real
end

WeightDecay() = WeightDecay(0)

function apply!(o::WeightDecay, x, Δ)
  wd = o.wd
  @. Δ += wd * data(x)
end


text:
const PyJuliaArrayValue_Type__ref = Ref(PyNULL)
PyJuliaArrayValue_Type() = begin
    ptr = PyJuliaArrayValue_Type__ref[]
    if isnull(ptr)
        c = []
        base = PyJuliaAnyValue_Type()
        isnull(base) && return PyNULL
        t = fill(
            PyType_Create(
                c,
                name = "juliaaa.ArrayValue",
                base = base,
                as_mapping = (
                    subscript = pyjlarray_getitem,
                    ass_subscript = pyjlarray_setitem,
                ),
                as_buffer = (
                    get = pyjlarray_get_buffer,
                    release = pyjlarray_release_buffer,
                ),
                getset = [
                    (name = "ndim", get = pyjlarray_ndim),
                    (name = "shape", get = pyjlarray_shape),
                    (name = "__array_interface__", get = pyjlarray_array_interface),
                ],
                methods = [
                    (name = "copy", flags = Py_METH_NOARGS, meth = pyjlarray_copy),
                    (name = "reshape", flags = Py_METH_O, meth = pyjlarray_reshape),
                    (name = "__array__", flags = Py_METH_NOARGS, meth = pyjlarray_array),
                ],
            ),
        )
        ptr = PyPtr(pointer(t))
        err = PyType_Ready(ptr)
        ism1(err) && return PyNULL
        for abcf in (CONFIG.version < v"3.6" ? (PyContainerABC_Type, PyIterableABC_Type, PySizedABC_Type) : (PyCollectionABC_Type,))
            abc = abcf()
            isnull(abc) && return PyNULL
            ism1(PyABC_Register(ptr, abc)) && return PyNULL
        end
        PYJLGCCACHE[ptr] = push!(c, t)
        PyJuliaArrayValue_Type__ref[] = ptr
    end
    ptr
end

PyJuliaArrayValue_New(x::AbstractArray) = PyJuliaValue_New(PyJuliaArrayValue_Type(), x)
PyJuliaValue_From(x::AbstractArray) = PyJuliaArrayValue_New(x)

pyjl_getaxisindex(x::AbstractUnitRange{<:Integer}, ko::PyPtr) = begin
    if PySlice_Check(ko)
        ao, co, bo = PySimpleObject_GetValue(ko, Tuple{PyPtr,PyPtr,PyPtr})
        # start
        r = PyObject_TryConvert(ao, Union{Int,Nothing})
        r == -1 && return PYERR()
        r == 0 && (
            PyErr_SetString(PyExc_TypeError(), "slice components must be integers"); return PYERR()
        )
        a = takeresult(Union{Int,Nothing})
        # step
        r = PyObject_TryConvert(bo, Union{Int,Nothing})
        r == -1 && return PYERR()
        r == 0 && (
            PyErr_SetString(PyExc_TypeError(), "slice components must be integers"); return PYERR()
        )
        b = takeresult(Union{Int,Nothing})
        # stop
        r = PyObject_TryConvert(co, Union{Int,Nothing})
        r == -1 && return PYERR()
        r == 0 && (
            PyErr_SetString(PyExc_TypeError(), "slice components must be integers"); return PYERR()
        )
        c = takeresult(Union{Int,Nothing})
        # step defaults to 1
        b′ = b === nothing ? 1 : b
        if a === nothing && c === nothing
            # when neither is specified, start and stop default to the full range,
            # which is reversed when the step is negative
            if b′ > 0
                a′ = Int(first(x))
                c′ = Int(last(x))
            elseif b′ < 0
                a′ = Int(last(x))
                c′ = Int(first(x))
            else
                PyErr_SetString(PyExc_ValueError(), "step must be non-zero")
                return PYERR()
            end
        else
            # start defaults
            a′ = Int(a === nothing ? first(x) : a < 0 ? (last(x) + a + 1) : (first(x) + a))
            c′ = Int(
                c === nothing ? last(x) :
                c < 0 ? (last(x) + 1 + c - sign(b′)) : (first(x) + c - sign(b′)),
            )
        end
        r = a′:b′:c′
        if !checkbounds(Bool, x, r)
            PyErr_SetString(PyExc_IndexError(), "array index out of bounds")
            return PYERR()
        end
        r
    else
        r = PyObject_TryConvert(ko, Int)
        ism1(r) && return PYERR()
        r == 0 && (
            PyErr_SetString(
                PyExc_TypeError(),
                "index must be slice or integer, got $(PyType_Name(Py_Type(ko)))",
            );
            return PYERR()
        )
        k = takeresult(Int)
        k′ = k < 0 ? (last(x) + k + 1) : (first(x) + k)
        checkbounds(Bool, x, k′) || (
            PyErr_SetString(PyExc_IndexError(), "array index out of bounds"); return PYERR()
        )
        k′
    end
end

pyjl_getarrayindices(x::AbstractArray, ko::PyPtr) = begin
    kos = PyTuple_Check(ko) ? [PyTuple_GetItem(ko, i - 1) for i = 1:PyTuple_Size(ko)] : [ko]
    length(kos) == ndims(x) || (
        PyErr_SetString(
            PyExc_TypeError(),
            "expecting exactly $(ndims(x)) indices, got $(length(kos))",
        );
        return PYERR()
    )
    ks = []
    for (i, ko) in enumerate(kos)
        k = pyjl_getaxisindex(axes(x, i), ko)
        k === PYERR() && return PYERR()
        push!(ks, k)
    end
    ks
end

pyjlarray_getitem(xo::PyPtr, ko::PyPtr) = begin
    x = PyJuliaValue_GetValue(xo)::AbstractArray
    k = pyjl_getarrayindices(x, ko)
    k === PYERR() && return PyNULL
    try
        if all(x -> x isa Int, k)
            PyObject_From(x[k...])
        else
            PyObject_From(view(x, k...))
        end
    catch err
        if err isa BoundsError && err.a === x
            PyErr_SetStringFromJuliaError(PyExc_IndexError(), err)
        else
            PyErr_SetJuliaError(err)
        end
        PyNULL
    end
end

pyjlarray_setitem(xo::PyPtr, ko::PyPtr, vo::PyPtr) = begin
    x = PyJuliaValue_GetValue(xo)::AbstractArray
    k = pyjl_getarrayindices(x, ko)
    k === PYERR() && return Cint(-1)
    try
        if isnull(vo)
            deleteat!(x, k...)
            Cint(0)
        else
            ism1(PyObject_Convert(vo, eltype(x))) && return Cint(-1)
            v = takeresult(eltype(x))
            if all(x -> x isa Int, k)
                x[k...] = v
                Cint(0)
            else
                PyErr_SetString(PyExc_TypeError(), "multiple assignment not supported")
                Cint(-1)
            end
        end
    catch err
        if err isa BoundsError && err.a === x
            PyErr_SetStringFromJuliaError(PyExc_IndexError(), err)
        elseif err isa MethodError && err.f === deleteat!
            PyErr_SetStringFromJuliaError(PyExc_TypeError(), err)
        else
            PyErr_SetJuliaError(err)
        end
        Cint(-1)
    end
end

pyjlarray_ndim(xo::PyPtr, ::Ptr) = begin
    x = PyJuliaValue_GetValue(xo)::AbstractArray
    PyObject_From(ndims(x))
end

pyjlarray_shape(xo::PyPtr, ::Ptr) = begin
    x = PyJuliaValue_GetValue(xo)::AbstractArray
    PyObject_From(size(x))
end

pyjlarray_copy(xo::PyPtr, ::PyPtr) =
    try
        x = PyJuliaValue_GetValue(xo)::AbstractArray
        PyObject_From(copy(x))
    catch err
        PyErr_SetJuliaError(err)
        PyNULL
    end

pyjlarray_reshape(xo::PyPtr, arg::PyPtr) =
    try
        x = PyJuliaValue_GetValue(xo)::AbstractArray
        r = PyObject_TryConvert(arg, Union{Int,Tuple{Vararg{Int}}})
        r == -1 && return PyNULL
        r == 0 && (
            PyErr_SetString(
                PyExc_TypeError(),
                "shape must be an integer or tuple of integers",
            );
            return PyNULL
        )
        PyObject_From(reshape(x, takeresult()...))
    catch err
        PyErr_SetJuliaError(err)
        PyNULL
    end

### Buffer Protocol

isflagset(flags, mask) = (flags & mask) == mask

const PYJLBUFCACHE = Dict{Ptr{Cvoid},Any}()

pyjl_get_buffer_impl(
    o::PyPtr,
    buf::Ptr{Py_buffer},
    flags,
    ptr,
    elsz,
    len,
    ndim,
    fmt,
    sz,
    strds,
    mutable,
) = begin
    b = UnsafePtr(buf)
    c = []

    # not influenced by flags: obj, buf, len, itemsize, ndim
    b.obj[] = C_NULL
    b.buf[] = ptr
    b.itemsize[] = elsz
    b.len[] = elsz * len
    b.ndim[] = ndim

    # readonly
    if isflagset(flags, PyBUF_WRITABLE)
        if mutable
            b.readonly[] = 1
        else
            PyErr_SetString(PyExc_BufferError(), "not writable")
            return Cint(-1)
        end
    else
        b.readonly[] = mutable ? 0 : 1
    end

    # format
    if isflagset(flags, PyBUF_FORMAT)
        b.format[] = cacheptr!(c, fmt)
    else
        b.format[] = C_NULL
    end

    # shape
    if isflagset(flags, PyBUF_ND)
        b.shape[] = cacheptr!(c, Py_ssize_t[sz...])
    else
        b.shape[] = C_NULL
    end

    # strides
    if isflagset(flags, PyBUF_STRIDES)
        b.strides[] = cacheptr!(c, Py_ssize_t[strds...])
    else
        if Python.size_to_cstrides(elsz, sz...) != strds
            PyErr_SetString(PyExc_BufferError(), "not C contiguous and strides not requested")
            return Cint(-1)
        end
        b.strides[] = C_NULL
    end

    # check contiguity
    if isflagset(flags, PyBUF_C_CONTIGUOUS)
        if Python.size_to_cstrides(elsz, sz...) != strds
            PyErr_SetString(PyExc_BufferError(), "not C contiguous")
            return Cint(-1)
        end
    end
    if isflagset(flags, PyBUF_F_CONTIGUOUS)
        if Python.size_to_fstrides(elsz, sz...) != strds
            PyErr_SetString(PyExc_BufferError(), "not Fortran contiguous")
            return Cint(-1)
        end
    end
    if isflagset(flags, PyBUF_ANY_CONTIGUOUS)
        if Python.size_to_cstrides(elsz, sz...) != strds &&
           Python.size_to_fstrides(elsz, sz...) != strds
            PyErr_SetString(PyExc_BufferError(), "not contiguous")
            return Cint(-1)
        end
    end

    # suboffsets
    b.suboffsets[] = C_NULL

    # internal
    cptr = Base.pointer_from_objref(c)
    PYJLBUFCACHE[cptr] = c
    b.internal[] = cptr

    # obj
    Py_IncRef(o)
    b.obj[] = o
    Cint(0)
end

pyjlarray_isbufferabletype(::Type{T}) where {T} = T in (
    Int8,
    Int16,
    Int32,
    Int64,
    UInt8,
    UInt16,
    UInt32,
    UInt64,
    Float16,
    Float32,
    Float64,
    Complex{Float16},
    Complex{Float32},
    Complex{Float64},
    Bool,
    Ptr{Cvoid},
)
pyjlarray_isbufferabletype(::Type{T}) where {T<:Tuple} =
    isconcretetype(T) &&
    Python.allocatedinline(T) &&
    all(pyjlarray_isbufferabletype, fieldtypes(T))
pyjlarray_isbufferabletype(::Type{NamedTuple{names,T}}) where {names,T} =
    pyjlarray_isbufferabletype(T)

_pyjlarray_get_buffer(xo, buf, flags, x::AbstractArray) =
    try
        if pyjlarray_isbufferabletype(eltype(x))
            pyjl_get_buffer_impl(
                xo,
                buf,
                flags,
                Base.unsafe_convert(Ptr{eltype(x)}, x),
                sizeof(eltype(x)),
                length(x),
                ndims(x),
                Python.pybufferformat(eltype(x)),
                size(x),
                strides(x) .* Python.aligned_sizeof(eltype(x)),
                Python.ismutablearray(x),
            )
        else
            error("element type is not bufferable")
        end
    catch err
        PyErr_SetString(
            PyExc_BufferError(),
            "Buffer protocol not supported by Julia '$(typeof(x))' (details: $err)",
        )
        Cint(-1)
    end

pyjlarray_get_buffer(xo::PyPtr, buf::Ptr{Py_buffer}, flags::Cint) =
    _pyjlarray_get_buffer(xo, buf, flags, PyJuliaValue_GetValue(xo)::AbstractArray)

pyjlarray_release_buffer(xo::PyPtr, buf::Ptr{Py_buffer}) = begin
    delete!(PYJLBUFCACHE, UnsafePtr(buf).internal[!])
    nothing
end

### Array Interface

pyjlarray_isarrayabletype(::Type{T}) where {T} = T in (
    UInt8,
    Int8,
    UInt16,
    Int16,
    UInt32,
    Int32,
    UInt64,
    Int64,
    Bool,
    Float16,
    Float32,
    Float64,
    Complex{Float16},
    Complex{Float32},
    Complex{Float64},
)
pyjlarray_isarrayabletype(::Type{T}) where {T<:Tuple} =
    isconcretetype(T) &&
    Python.allocatedinline(T) &&
    all(pyjlarray_isarrayabletype, T.parameters)
pyjlarray_isarrayabletype(::Type{NamedTuple{names,types}}) where {names,types} =
    pyjlarray_isarrayabletype(types)

pyjlarray_array_interface(xo::PyPtr, ::Ptr{Cvoid}) =
    _pyjlarray_array_interface(PyJuliaValue_GetValue(xo)::AbstractArray)

PyDescrObject_From(fields) = begin
    ro = PyList_New(0)
    for (name, descr) in fields
        descro = descr isa String ? PyUnicode_From(descr) : PyDescrObject_From(descr)
        isnull(descro) && (Py_DecRef(ro); return PyNULL)
        fieldo = PyTuple_From((name, PyObjectRef(descro)))
        Py_DecRef(descro)
        isnull(fieldo) && (Py_DecRef(ro); return PyNULL)
        err = PyList_Append(ro, fieldo)
        ism1(err) && (Py_DecRef(ro); return PyNULL)
    end
    ro
end

_pyjlarray_array_interface(x::AbstractArray) =
    try
        if pyjlarray_isarrayabletype(eltype(x))
            # gather information
            shape = size(x)
            data = (UInt(Base.unsafe_convert(Ptr{eltype(x)}, x)), !Python.ismutablearray(x))
            strides = Base.strides(x) .* Python.aligned_sizeof(eltype(x))
            version = 3
            typestr, descr = Python.pytypestrdescr(eltype(x))
            isempty(typestr) && error("invalid element type")
            # make the dictionary
            d = PyDict_From(
                Dict(
                    "shape" => shape,
                    "typestr" => typestr,
                    "data" => data,
                    "strides" => strides,
                    "version" => version,
                ),
            )
            isnull(d) && return PyNULL
            if descr !== nothing
                descro = PyDescrObject_From(descr)
                err = PyDict_SetItemString(d, "descr", descro)
                Py_DecRef(descro)
                ism1(err) && (Py_DecRef(d); return PyNULL)
            end
            d
        else
            error("invalid element type")
        end
    catch err
        PyErr_SetString(PyExc_AttributeError(), "__array_interface__")
        PyNULL
    end

pyjlarray_array(xo::PyPtr, ::PyPtr) = begin
    # get an array-like version of x
    if PyObject_HasAttrString(xo, "__array_interface__") == 0
        # convert to a PyObjectArray
        x = PyJuliaValue_GetValue(xo)::AbstractArray
        y = try
            Python.PyObjectArray(x)
        catch err
            PyErr_SetJuliaError(err)
            return PyNULL
        end
        yo = PyJuliaArrayValue_New(y)
        isnull(yo) && return PyNULL
    else
        # already supports the array interface
        Py_IncRef(xo)
        yo = xo
    end
    npo = PyImport_ImportModule("numpy")
    if isnull(npo)
        if PyErr_IsSet(PyExc_ImportError())
            # numpy not installed, return an array-like object
            PyErr_Clear()
            yo
        else
            Py_DecRef(yo)
            PyNULL
        end
    else
        # numpy installed, return a proper numpy array
        aao = PyObject_GetAttrString(npo, "array")
        Py_DecRef(npo)
        isnull(aao) && (Py_DecRef(yo); return PyNULL)
        ao = PyObject_CallNice(aao, PyObjectRef(yo))
        Py_DecRef(aao)
        Py_DecRef(yo)
        ao
    end
end


text:
Base.eltype(::TurbulenceClosure{T}) where T = T

"""
    typed_keyword_constructor(T, Closure; kwargs...)

Return an object `Closure` with fields provided by `kwargs`
converted to type `T`. Mainly provided for converting between
different float types when working with constructors associated
with types defined via `Base.@kwdef`.
"""
function typed_keyword_constructor(T, Closure; kwargs...)
    closure = Closure(; kwargs...)
    names = fieldnames(Closure)
    vals = [getproperty(closure, name) for name in names]
    return Closure{T}(vals...)
end

function Base.convert(::TurbulenceClosure{T2}, closure::TurbulenceClosure{T1}) where {T1, T2}
    paramdict = Dict((p, convert(T2, getproperty(closure, p))) for p in propertynames(closure))
    Closure = typeof(closure).name.wrapper
    return Closure(T2; paramdict...)
end


text:
# This file is a part of RunStatistics.jl, licensed under the MIT License (MIT).

using RunStatistics
using Test
using Distributions

# I measured the timing for N=96 with the exact and approximate solution. See if julia is really faster. It seems some experiment with threading was done but without threads, I'd be surprised if it's faster than the C++ code built in release mode

@testset "squares_approx" begin

    T_obs = 20
    L = 800
    N = 80
    n = 10
    Ns = [N, n]
    epsp_1 = 10^(-2)
    epsp_2 = 0

    @test RunStatistics.squares_pvalue_approx(T_obs, Ns, epsp_1) ≈ 0.07082230261169509
    @test RunStatistics.squares_cdf_approx(T_obs, Ns, epsp_1) ≈ 0.9291776973883049

    @test RunStatistics.squares_cdf_approx(T_obs, L, epsp_1) ≈ 0.9291776973883049
    @test RunStatistics.squares_pvalue_approx(T_obs, L, epsp_1) ≈ 0.07082230261169509

    @test RunStatistics.squares_pvalue_approx(T_obs, Ns, epsp_2) == 0.0708223271587789
    @test RunStatistics.squares_cdf_approx(T_obs, Ns, epsp_2) == 0.9291776728412211

    @test RunStatistics.squares_cdf_approx(T_obs, L, epsp_2) == 0.9291776728412211
    @test RunStatistics.squares_pvalue_approx(T_obs, L, epsp_2) == 0.0708223271587789
end

@testset "squares_approx_split" begin

    Tobs = 15.5
    N = 12
    n = 2

    F = squares_cdf(Tobs, N)
    approx = F*F / (1.0 + RunStatistics.Delta(Tobs, N, N))

    @test abs(approx - squares_cdf_approx(Tobs, [N, n])) <= 1e-15
end

function test_on_grid(K::Integer, N::Integer, n::Real, eps::Real)

    for i in 1:K
        Tobs = 20 + 2 * i
        @test abs(squares_cdf_approx(Tobs, [N, n]) - squares_cdf(Tobs, n * N)) <= eps
    end
end

@testset "squares_approx_grid" begin

    test_on_grid(15, 40, 2, 2e-7)
    # test_on_grid(15, 40, 3, 4e-7) expensive
end

@testset "squares_approx_hH" begin

    Tobs = 15.5
    x = 3.3
    N = 12

    @test abs(RunStatistics.h(Tobs, N) - 0.000373964) <= 1e-8
    @test abs(RunStatistics.H(Tobs-x, Tobs, N) -  0.00245352) <= 1e-8
    @test abs(RunStatistics.Delta(Tobs, N, N) - 0.00175994) <= 1e-8
end

@testset "squares_approx_cdf" begin

    @test abs(cdf(Chisq(12), 15.5) - 0.784775) <= 1e-6
end

@testset "squares_approx_interpolate" begin
    Tobs = 32

    F71 = squares_cdf_approx(Tobs, [71, 5])
    F88 = squares_cdf_approx(Tobs, [88, 4])
    F89 = squares_cdf_approx(Tobs, [89, 4])
    F100 = squares_cdf_approx(Tobs, [Int64(100), 3.55])

    @test abs(F71 - F100) <= 3e-9
    @test abs(F71 - (F88 + 3*F89)/4) <= 2e-9
end

@testset "squares_approx_bound_error" begin

    Tobs = 8
    @test abs(squares_cdf(Tobs, 100) - (squares_cdf(Tobs, 50) ^ 2 - squares_cdf(Tobs, 100) * RunStatistics.Delta(Tobs, 100, 100))) <= 1e-3
end


text:
using Distributions, InvertedIndices, StatsBase, ProgressMeter, RecipesBase

export iex_mcmc_mode, iex_mcmc_gamma, aux_term_eval_mode, rand_reflect

function rand_reflect(x, ε, l, u)
    ξ = ε * (2*rand() - 1)
    y = x + ξ
    if y < l
        return 2*l - y
    elseif y > u
        return 2*u - y
    else
        return y
    end
end

# Auxiliary sampler types
function aux_term_eval_mode(
    mcmc::SpfInvolutiveMcmcEdit,
    I_curr::Path{T}, I_prop::Path{T},
    posterior::SpfPosterior{T},
    γ_curr::Float64,
    aux_model::SPF{T}
    ) where {T<:Union{Int,String}}

    @assert aux_model.mode == I_prop "Auxiliary model must have proposed value as mode."

    aux_log_lik_ratio = 0.0

    # Pointers for clear exposition
    aux_curr = mcmc.curr
    aux_prop = mcmc.prop
    ind_del = mcmc.ind_del
    ind_add = mcmc.ind_add
    vals = mcmc.vals
    ν = mcmc.ν

    # We are sampling from model with pars (I_prop, γ_curr)

    # Initialise at I_prop
    copy!(aux_curr, I_prop)
    copy!(aux_prop, aux_curr)
    # ind_del = zeros(Int, mcmc.ν)
    # ind_add = zeros(Int, mcmc.ν)
    # vals = rand(posterior.V, mcmc.ν)
    count = 0

    # burn_in
    for i in 1:mcmc.burn_in
        was_accepted = accept_reject_edit!(
            aux_curr, aux_prop,
            aux_model,
            ν,
            ind_del, ind_add, vals
            )
        count+=was_accepted
    end
    # For sample size, do lag-1 steps with no action, then accept reject and
    # increment the auxiliary term
    for i in 1:posterior.sample_size
        for j in 1:(mcmc.lag-1)
            was_accepted = accept_reject_edit!(
            aux_curr, aux_prop,
            aux_model,
            ν,
            ind_del, ind_add, vals
            )
        count+=was_accepted
        end
        # Now value TO BE USED IN AUX TERM
        was_accepted = accept_reject_edit!(
            aux_curr, aux_prop,
            aux_model,
            ν,
            ind_del, ind_add, vals
            )
        count+=was_accepted

        aux_log_lik_ratio += -γ_curr * (
            posterior.dist(aux_curr, I_curr)
            - posterior.dist(aux_curr, I_prop)
            )

    end
    return aux_log_lik_ratio

end

function iex_mcmc_mode(
    posterior::SpfPosterior{T},
    mcmc_sampler::SpfMcmcSampler,
    γ_fixed::Float64;
    I_init::Path{T}=sample_frechet_mean(posterior.data, posterior.dist)[1],
    desired_samples::Int=100,
    burn_in=0,
    lag=1,
    ν=2
    )  where {T<:Union{Int, String}}

    # Find required length of chain
    req_samples = burn_in + 1 + (desired_samples - 1) * lag

    iter = Progress(req_samples, 1, "Chain for γ = $(γ_fixed) and n = $(posterior.sample_size) (mode marginal)....")  # Loading bar. Minimum update interval: 1 second

    # Intialise
    I_sample = Vector{Path{T}}() # Storage of samples
    I_curr = copy(I_init)
    γ_curr = γ_fixed
    # aux_data = Vector{Path{T}}(undef, length(posterior.data))
    I_count = 0
    aux_data = fill(T[], posterior.sample_size)
    aux_model = SPF(I_curr, γ_curr, posterior.dist, posterior.V, posterior.K)

    # Vertex distribution for proposal
    μ = get_vertex_proposal_dist(posterior)
    if T == String
        vdist = StringCategorical(posterior.V, μ)
    elseif T == Int
        vdist = Categorical(μ)
    else
        error("Path eltype not recognised for defining vertex proposal dist.")
    end

    for i in 1:req_samples
        # Sample proposed mode
        n = length(I_curr)
        δ = rand(1:ν)
        d = rand(0:min(n-1,δ))

        I_prop, log_ratio = imcmc_prop_sample_edit_informed(I_curr, δ, d, vdist)  # Sample new mode

        # Generate auxiliary data (centered on proposal)
        aux_model = SPF(I_prop, γ_curr, posterior.dist, posterior.V, posterior.K)
        # aux_data = mcmc_sampler(aux_model, desired_samples=posterior.sample_size).sample
        sample!(mcmc_sampler, aux_model, aux_data)
        # Accept reject

        log_lik_ratio = -γ_curr * (
            sum_of_dists(posterior.data, I_prop, posterior.dist)
            - sum_of_dists(posterior.data, I_curr, posterior.dist)
        )
        aux_log_lik_ratio = -γ_curr * (
            sum_of_dists(aux_data, I_curr, posterior.dist)
            - sum_of_dists(aux_data, I_prop, posterior.dist)
        )
        # aux_log_lik_ratio = aux_term_eval_mode(
        #     mcmc_sampler,
        #     I_curr, I_prop,
        #     posterior,
        #     γ_curr,
        #     aux_model
        # )
        log_α = (
            posterior.I_prior.γ * (
                sum_of_dists(posterior.I_prior, I_curr)
                - sum_of_dists(posterior.I_prior, I_prop)
            )
            + log_lik_ratio + aux_log_lik_ratio + log_ratio
        )

        if log(rand()) < log_α
            copy!(I_curr, I_prop)
            I_count += 1
        else
            copy!(I_prop, I_curr)
        end
        push!(I_sample, copy(I_curr))
        next!(iter)

    end
    performance_measures = Dict(
        "Mode acceptance probability" => I_count/req_samples
    )

    output = SpfPosteriorModeConditionalMcmcOutput(
        γ_fixed,
        I_sample[(burn_in+1):lag:end],
        posterior.dist,
        posterior.I_prior,
        posterior.data,
        performance_measures
    )
    return output
end

function iex_mcmc_gamma(
    posterior::SpfPosterior{T},
    mcmc_sampler::SpfMcmcSampler,
    I_fixed::Path{T};
    γ_init::Float64=3.9,
    desired_samples::Int=100,
    burn_in::Int=0,
    lag::Int=1,
    ε::Float64=0.1
    )  where {T<:Union{Int, String}}

    # Find required length of chain
    req_samples = burn_in + 1 + (desired_samples - 1) * lag

    iter = Progress(req_samples, 1, "Chain for mode = $(I_fixed) and n = $(posterior.sample_size) (dispersion marginal)....")  # Loading bar. Minimum update interval: 1 second

    # Intialise
    γ_sample = Vector{Float64}(undef, req_samples) # Storage of samples
    I_curr = I_fixed
    γ_curr = γ_init
    aux_data = fill(T[], posterior.sample_size)
    count = 0

    suff_stat = sum_of_dists(posterior.data, I_curr, posterior.dist) # Sum of distance to mode, fixed throughout

    for i in 1:req_samples
        # Sample proposed γ

        γ_prop = rand_reflect(γ_curr, ε, 0.0, Inf)

        # Generate auxiliary data (centered on proposal)
        aux_model = SPF(I_curr, γ_prop, posterior.dist, posterior.V, posterior.K)
        # aux_data = mcmc_sampler(a_model, desired_samples=posterior.sample_size).sample
        sample!(mcmc_sampler, aux_model, aux_data)
        # Accept reject

        log_lik_ratio = (γ_curr - γ_prop) * suff_stat
        aux_log_lik_ratio = (γ_prop - γ_curr) * sum_of_dists(aux_data, I_curr, posterior.dist)

        log_α = (
            logpdf(posterior.γ_prior, γ_prop)
            - logpdf(posterior.γ_prior, γ_curr)
            + log_lik_ratio + aux_log_lik_ratio
        )

        if log(rand()) < log_α
            γ_sample[i] = γ_prop
            γ_curr = γ_prop
            count += 1
        else
            γ_sample[i] = γ_curr
        end
        next!(iter)

    end
    performance_measures = Dict(
        "Mode acceptance probability" => count/req_samples
    )

    output = SpfPosteriorDispersionConditionalMcmcOutput(
        I_fixed,
        γ_sample[(burn_in+1):lag:end],
        posterior.γ_prior,
        posterior.data,
        performance_measures
    )
    return output
end