text:
# Variables in Xpress model

# variable kinds

function fixInf(val::Float64)
    if val == Inf
        return XPRS_PLUSINFINITY
    elseif val == -Inf
        return XPRS_MINUSINFINITY
    else
        return val
    end
end
function fixInfs!(vals::Vector{Float64})
    for i in eachindex(vals)
        vals[i] = fixInf(vals[i])
    end
end

# add_var! single var
function add_var!(model::Model, newnz::Int, mrwind::Vector{Int},dmatval::Vector{Float64}, objx::Float64, lb::Float64, ub::Float64)
    lb = fixInf(lb)
    ub = fixInf(ub)

    ret = @xprs_ccall(addcols, Cint, (
        Ptr{Void},    # model
        Cint,         # number of new cols
        Cint,         # number of nonzeros in added cols
        Ptr{Float64}, # coefs in objective
        Ptr{Cint},     # mstart ????
        Ptr{Cint},     # mrwind: row indexes of each col
        Ptr{Float64},      # elemnts values
        Ptr{Float64},      # lb
        Ptr{Float64}    # ub
        ),
        model.ptr_model, 1, newnz, fvec(objx), ivec(0), ivec(mrwind)-1, fvec(dmatval), fvec(lb), fvec(ub))

    if ret != 0
        throw(XpressError(model))
    end
    nothing
end

function chgcoltype!(model::Model,colnum::Int,vtype)
    ret = @xprs_ccall(chgcoltype,Cint,(Ptr{Void},Cint,Ptr{Cint},Ptr{Cchar}),
        model.ptr_model,1,Cint[colnum-1],Cchar[XPRS_INTEGER])
    if  0 != ret
        throw(XpressError(model))
    end
end
function chgcoltypes!(model::Model,colnums::Vector{Int},vtypes::Vector)
    n=length(colnums)
    ret = @xprs_ccall(chgcoltype,Cint,(
        Ptr{Void},
        Cint,
        Ptr{Cint},
        Ptr{Cchar}),
        model.ptr_model,n,ivec(colnums-1),cvec(vtypes) )
    if 0 != ret
        throw(XpressError(model))
    end
end
#only add in obj
function add_var!(model::Model, vtype::Cchar, c::Float64, lb::Float64, ub::Float64)
    lb = fixInf(lb)
    ub = fixInf(ub)

    ret = @xprs_ccall(addcols, Cint, (
        Ptr{Void},    # model
        Cint,         # number of new cols
        Cint,         # number of nonzeros in added cols
        Ptr{Float64}, # coefs in objective
        Ptr{Cint},     # mstart ????
        Ptr{Cint},     # mrwind: row indexes of each col
        Ptr{Float64},      # elemnts values
        Ptr{Float64},      # lb
        Ptr{Float64}    # ub
        ),
        model.ptr_model, 1, 0, fvec(c), ivec(0), C_NULL, C_NULL, fvec(lb), fvec(ub))
    if ret != 0
        throw(XpressError(model))
    end

    if vtype == XPRS_INTEGER
        chgcoltype!(model,get_intattr(model,XPRS_COLS),XPRS_INTEGER)
    elseif vtype ==XPRS_BINARY
        chgcoltype!(model,get_intattr(model,XPRS_COLS),XPRS_BINARY)
    end

    nothing
end

function add_var!(model::Model, vtype::GChars, c::Real, lb::Real, ub::Real)
    add_var!(model, cchar(vtype), Float64(c), Float64(lb), Float64(ub))
end
add_var!(model::Model, vtype::GChars, c::Real) = add_var!(model, vtype, c, XPRS_MINUSINFINITY, XPRS_PLUSINFINITY)

add_cvar!(model::Model, c::Real, lb::Real, ub::Real) = add_var!(model, XPRS_CONTINUOUS, c, lb, ub)
add_cvar!(model::Model, c::Real) = add_cvar!(model, c, XPRS_MINUSINFINITY, XPRS_PLUSINFINITY)

add_bvar!(model::Model, c::Real) = add_var!(model, XPRS_BINARY, c, 0., 1.)

add_ivar!(model::Model, c::Real, lb::Real, ub::Real) = add_var!(model, XPRS_INTEGER, c, lb, ub)
add_ivar!(model::Model, c::Real) = add_ivar!(model, c, XPRS_MINUSINFINITY, XPRS_PLUSINFINITY)

# add_vars!
#multiple variables at once
function add_vars!(model::Model, vtypes::CVec, c::FVec, lb::FVec, ub::FVec)

    # check dimensions
    n = length(vtypes)
    _chklen(c, n)
    _chklen(lb, n)
    _chklen(ub, n)

    fixInfs!(lb)
    fixInfs!(ub)

    prev_cols = get_intattr(model,XPRS_COLS)
    # main call
    ret = @xprs_ccall(addcols, Cint, (
        Ptr{Void},    # model
        Cint,         # number of new cols
        Cint,         # number of nonzeros in added cols
        Ptr{Float64}, # coefs in objective
        Ptr{Cint},     # mstart ????
        Ptr{Cint},     # mrwind: row indexes of each col
        Ptr{Float64},      # elemnts values
        Ptr{Float64},      # lb
        Ptr{Float64}    # ub
        ),
        model.ptr_model, n, 0, c, ivec( collect(0:(n-1)) ), C_NULL, C_NULL, lb, ub)
                                #check
    if ret != 0
        throw(XpressError(model))
    end

    #vtypes = cvecx(vtypes, n)
    for i in 1:n
        if vtypes[i] == XPRS_INTEGER
            chgcoltype!(model,prev_cols+i,XPRS_INTEGER)
        elseif vtypes[i] ==XPRS_BINARY
            chgcoltype!(model,prev_cols+i,XPRS_BINARY)
        end
    end
    nothing
end

function add_vars!(model::Model, vtypes::GCharOrVec, c::Vector, lb::Bounds, ub::Bounds)
    n = length(c)
    add_vars!(model, cvecx(vtypes, n), fvec(c), fvecx(lb, n), fvecx(ub, n))
end

add_cvars!(model::Model, c::Vector, lb::Bounds, ub::Bounds) = add_vars!(model, XPRS_CONTINUOUS, c, lb, ub)
add_cvars!(model::Model, c::Vector) = add_cvars!(model, c, XPRS_MINUSINFINITY, XPRS_PLUSINFINITY)

add_bvars!(model::Model, c::Vector) = add_vars!(model, XPRS_BINARY, c, 0, 1)

add_ivars!(model::Model, c::Vector, lb::Bounds, ub::Bounds) = add_vars!(model, XPRS_INTEGER, c, lb, ub)
add_ivars!(model::Model, c::Vector) = add_ivars!(model, XPRS_INTEGER, c, XPRS_MINUSINFINITY, XPRS_PLUSINFINITY) #


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

@author: Michiel Stock
michielfmstock@gmail.com

Vec trick: multiplying vectors with Kronecker systems.
=#

import Base.ReshapedArray

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

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

Calculates the vector-matrix multiplication `K * v` and stores the result in
`x`, overwriting its existing value.
"""
function mul_vec_trick!(x::AbstractVector, A::AbstractKroneckerProduct, v::AbstractVector)
    M, N = getmatrices(A)
    a, b = size(M)
    c, d = size(N)
    e = length(v)
    f = length(x)

    V = vectrick_reshape(v, d, b)
    X = vectrick_reshape(x, c, a)
    if b * c * (a + d) < a * d * (b + c)
        mul!(X, N, V * transpose(M))
    else
        mul!(X, N * V, transpose(M))
    end
    return x
end

if VERSION < v"1.3.0-alpha.115"
    function mul_vec_trick!(x::AbstractVector, K::AbstractKroneckerSum, v::AbstractVector)
        A, B = getmatrices(K)
        a, b = size(A)
        c, d = size(B)
        e = length(v)
        f = length(x)

        V = reshape(v, d, b)
        X = reshape(x, c, a)
        mul!(X, V, transpose(A))
        X .+= B * V
        return x
    end
else # 5-arg mul! is available
    function mul_vec_trick!(x::AbstractVector, K::AbstractKroneckerSum, v::AbstractVector)
        A, B = getmatrices(K)
        a, b = size(A)
        c, d = size(B)
        e = length(v)
        f = length(x)

        V = reshape(v, d, b)
        X = reshape(x, c, a)
        mul!(X, V, transpose(A))
        mul!(X, B, V, true, true)
        return x
    end
end # VERSION

function mul_vec_trick!(X::AbstractMatrix, A::GeneralizedKroneckerProduct, V::AbstractMatrix)
    @inbounds for i in eachindex(axes(X, 2), axes(V, 2))
        @views mul_vec_trick!(X[:, i], A, V[:, i])
    end
    return X
end

# SOLVING
function ldiv_vec_trick!(x::AbstractVector, A::AbstractKroneckerProduct, v::AbstractVector)
    M, N = getmatrices(A)
    b, a = size(M)
    d, c = size(N)
    e = length(v)
    f = length(x)

    V = vectrick_reshape(v, d, b)
    X = vectrick_reshape(x, c, a)
    if b * c * (a + d) < a * d * (b + c)
        S = N \ V
        copyto!(transpose(X), M \ transpose(S))
    else
        S = M \ transpose(V)
        copyto!(X, N \ transpose(S))
    end
    return x

    # size(K, 1) != length(v) && throw(DimensionMismatch("size(K, 1) != length(v)"))
    # C = reshape(v, size(K.B, 1), size(K.A, 1)) # matricify
    # return vec((K.B \ C) / K.A') #(A ⊗ B)vec(X) = vec(C) <=> BXA' = C => X = B^{-1} C A'^{-1}
end

function ldiv_vec_trick!(X::AbstractMatrix, A::AbstractKroneckerProduct, V::AbstractMatrix)
    @inbounds for i in eachindex(axes(X, 2), axes(V, 2))
        @views ldiv_vec_trick!(X[:, i], A, V[:, i])
    end
    return X
end

function check_compatible_sizes(C::AbstractVecOrMat, A::AbstractMatrix, B::AbstractVecOrMat, mul=true)
    # when performing a division (mul=false), A acts as a matrix with reversed dimensions
    m, n = mul ? size(A) : reverse(size(A))

    if n != size(B, 1)
        throw(DimensionMismatch(
            "A has size $(size(A)), B has size $(size(B))"
        ))
    elseif size(C) != (m, size(B)[2:end]...)
        throw(DimensionMismatch(
            "A has size $(size(A)), B has size $(size(B)), C has size $(size(C))"
        ))
    end
end

function mul!(C::AbstractMatrix, A::AbstractKroneckerProduct, D::Diagonal)
    check_compatible_sizes(C, A, D)
    @inbounds for j in axes(C, 2)
        @views C[:, j] = A[:, j] * D[j, j]
    end
    return C
end
function mul!(C::AbstractMatrix, D::Diagonal, A::AbstractKroneckerProduct)
    check_compatible_sizes(C, D, A)
    @inbounds for i in axes(C, 1)
        @views C[i, :] = D[i, i] * A[i, :]
    end
    return C
end

for TC in [:AbstractVector, :AbstractMatrix],
    TB in [:($TC), :(Transpose{T, <:$TC{T}} where T), :(Adjoint{T, <:$TC{T}} where T)]

    @eval function mul!(C::$TC, A::AbstractKroneckerProduct, B::$TB)
        check_compatible_sizes(C, A, B)

        factors = getallfactors(A)

        if length(factors) == 2
            return mul_vec_trick!(C, A, B)
        elseif all(issquare, factors)
            return _kron_mul_fast_square!(C, B, factors)
        else
            return _kron_mul_fast_rect!(C, B, factors)
        end
    end

    @eval function ldiv!(C::$TC, A::AbstractKroneckerProduct, B::$TB)
        check_compatible_sizes(C, A, B, false)

        matrices = getallfactors(A)

        # TODO: should figure out how to avoid this if possible, leave
        #       factorization up to the user/Julia
        factors = ntuple(length(matrices)) do i
            m = matrices[i]
            return (m isa Factorization) ? m : factorize(m)
        end

        if length(factors) == 2
            return ldiv_vec_trick!(C, A, B)
        elseif all(issquare, factors)
            return _kron_ldiv_fast_square!(C, B, factors)
        else
            return _kron_ldiv_fast_rect!(C, B, factors)
        end
    end

    @eval function mul!(C::$TC, A::AbstractKroneckerSum, B::$TB)
        check_compatible_sizes(C, A, B)

        summands = getallsummands(A)

        if length(summands) == 2
            return mul_vec_trick!(C, A, B)
        else
            return _kronsum_mul_fast!(C, B, summands)
        end
    end
end

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

# explicit list of types that need to be disambiguated against
const MulMatTypes = [:Diagonal, :AbstractTriangular]

for T in [MulMatTypes; :AbstractMatrix]
    @eval function Base.:*(K::GeneralizedKroneckerProduct, M::$T)
        return mul!(Matrix{promote_type(eltype(M), eltype(K))}(undef, size(K, 1), size(M, 2)), K, M)
    end
end
for T in MulMatTypes
    @eval function Base.:*(M::$T, K::GeneralizedKroneckerProduct)
        return mul!(Matrix{promote_type(eltype(K), eltype(M))}(undef, size(M, 1), size(K, 2)), M, K)
    end
end

function Base.:*(v::AbstractMatrix, K::GeneralizedKroneckerProduct)
    out = Matrix{promote_type(eltype(v), eltype(K))}(undef, last(size(K)), first(size(v)))
    return transpose(mul!(out, transpose(K), transpose(v)))
end

function Base.:*(v::Adjoint{<:Number, <:AbstractVector}, K::GeneralizedKroneckerProduct)
    out = Vector{promote_type(eltype(v), eltype(K))}(undef, last(size(K)))
    return mul!(out, K', v.parent)'
end

function Base.:*(v::Transpose{<:Number, <:AbstractVector}, K::GeneralizedKroneckerProduct)
    out = Vector{promote_type(eltype(v), eltype(K))}(undef, last(size(K)))
    return transpose(mul!(out, transpose(K), v.parent))
end

# special multiplication methods for Kronecker products of Diagonal matrices
# It's usually better to convert these to Diagonal to use optimized multiplication methods
# instead of using the vec trick
const KroneckerDiagonal = Union{KronProdDiagonal, KronPowDiagonal}
Base.:*(K::KroneckerDiagonal, v::AbstractVector) = Diagonal(K) * v
Base.:*(K1::KroneckerDiagonal, K2::KroneckerDiagonal) = kronecker(map(*, getallfactors(K1), getallfactors(K2))...)
for T in [MulMatTypes; :AbstractMatrix]
    @eval Base.:*(K::KroneckerDiagonal, D::$T) = Diagonal(K) * D
    @eval Base.:*(D::$T, K::KroneckerDiagonal) = D * Diagonal(K)
end

for T in [:Adjoint, :Transpose]
    @eval Base.:*(A::$T{<:Number, <:AbstractVector}, K::KroneckerDiagonal) = A * Diagonal(K)
end

function LinearAlgebra.:\(K::AbstractKroneckerProduct, v::AbstractVector)
    return ldiv!(Vector{promote_type(eltype(v), eltype(K))}(undef, last(size(K))), K, v)
end

function LinearAlgebra.:\(K::AbstractKroneckerProduct, v::AbstractMatrix)
    return ldiv!(Matrix{promote_type(eltype(v), eltype(K))}(undef, last(size(K)), last(size(v))), K, v)
end

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

ext: jl


text:
########### "Simple" Problem Definitons ###########

@counted function simple1(x::Vector)
    return -x[1] * x[2] + 2.0 / (3.0 * sqrt(3.0))
end

@counted function simple1_gradient(x::Vector)
    return [-x[2], -x[1]]
end

@counted function simple1_constraints(x::Vector)
    return [x[1] + x[2]^2 - 1,
            -x[1] - x[2]]
end

function simple1_init()
    return rand(2) * 2.0
end

@counted function simple2(x::Vector)
    return (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2
end

@counted function simple2_gradient(x::Vector)
    storage = zeros(2)
    storage[1] = -2.0 * (1.0 - x[1]) - 400.0 * (x[2] - x[1]^2) * x[1]
    storage[2] = 200.0 * (x[2] - x[1]^2)
    return storage
end

@counted function simple2_constraints(x::Vector)
    return [(x[1]-1)^3 - x[2] + 1,
            x[1] + x[2] - 2]
end

function simple2_init()
    return rand(2) .* 2.0 .- 1.0
end

@counted function simple3(x::Vector)
    return x[1] - 2*x[2] + x[3] + sqrt(6.0)
end

@counted function simple3_gradient(x::Vector)
    return [1, -2, 1]
end

@counted function simple3_constraints(x::Vector)
    return [x[1]^2 + x[2]^2 + x[3]^2 - 1]
end

function simple3_init()
    b = 2.0 .* [1.0, -1.0, 0.0]
    a = -2.0 .* [1.0, -1.0, 0.0]
    return rand(3) .* (b-a) + a
end


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

"""
    struct BAT.FunnelDistribution <: Distribution{Multivariate,Continuous}

*Experimental feature, not part of stable public API.*

A funnel distribution (see
[Caldwell et al.](https://arxiv.org/abs/1808.08051) for definition).

Constructors:

* ```FunnelDistribution(; a::Real = 1.0, b::Real = 0.5, n::Integer = 3)```

Fields:

$(TYPEDFIELDS)
"""
struct FunnelDistribution{T<:Real,U<:Integer} <: Distribution{Multivariate,Continuous}
    "Variance of the dominant normal distribution."
    a::T

    "Variance of the supporting normal distributions."
    b::T

    "Number of dimensions."
    n::U
end

function FunnelDistribution(; a::Real = 1.0, b::Real = 0.5, n::Integer = 3)
    a_pr, b_pr = promote(a, b)
    FunnelDistribution(a_pr, b_pr, n)
end

Base.length(d::FunnelDistribution) = d.n
Base.eltype(d::FunnelDistribution) = Base.eltype(d.a)

Statistics.mean(d::FunnelDistribution) = zeros(d.n)

function Statistics.cov(dist::BAT.FunnelDistribution)
    cov(nestedview(rand(bat_determ_rng(), bat_sampler(dist), 10^5)))
end

StatsBase.params(d::FunnelDistribution) = (d.a, d.b, d.n)

function Distributions._logpdf(d::FunnelDistribution, x::AbstractArray)
    idxs = eachindex(x)
    s = logpdf(Normal(0, d.a^2), x[idxs[1]])
    @inbounds for i in idxs[2:end]
        s += logpdf(Normal(0, exp(2 * d.b * x[1])), x[i])
    end
    s
end

function Distributions._rand!(rng::AbstractRNG, d::FunnelDistribution, x::AbstractVector)
    idxs = eachindex(x)
    x[idxs[1]] = rand(rng, Normal(0, d.a^2))
    @inbounds for i in idxs[2:end]
        x[i] = rand(rng, Normal(0, exp(2 * d.b * x[1])))
    end
    return x
end


text:
function playground_window()
    problem = get_problem_desc("ZDT1")
    problem_renderer = problem_gui(problem)

    algorithm = get_algorithm_desc("MOEA/D")
    algorithm_renderer = algorithm_gui(algorithm)

    runtime_image_id = nothing
    runtime_image = nothing
    generation = 0
    igd = 0.0
    igds::Vector{Float64} = []
    igd_image_id = nothing
    igd_image = nothing
    hv = 0.0
    hvs::Vector{Float64} = []
    hv_image_id = nothing
    hv_image = nothing

    function ()
        CImGui.Begin("Playground")

        problem_renderer()

        algorithm_renderer()

        if CImGui.Button("Evaluate")
            igds = []
            hvs = []
            constraints = BoxConstraints(problem.bounds)
            fig_size = get_figure_size()
            state_callback = function (state)
                try
                    generation = state.iteration

                    fig = Plots.plot(
                        objectives(pfront(state))[:, 1],
                        objectives(pfront(state))[:, 2],
                        seriestype = :scatter,
                        legend = false,
                    )
                    runtime_image = fit_image(plot_image(fig), fig_size)

                    igd = MOEA.igd(pfront(state), problem.pfront)
                    push!(igds, igd)
                    fig = Plots.plot(
                        1:length(igds),
                        igds,
                        legend = false,
                    )
                    igd_image = fit_image(plot_image(fig), fig_size)

                    hv = hypervolume(objectives(pfront(state)), ones(size(objectives(pfront(state)))[2]))
                    push!(hvs, hv)
                    fig = Plots.plot(
                        1:length(hvs),
                        hvs,
                        legend = false,
                    )
                    hv_image = fit_image(plot_image(fig), fig_size)
                catch e
                    @error "Error in state callback!" exception=e
                    Base.show_backtrace(stderr, catch_backtrace())
                end
            end
            # TODO: initial population
            population_size = get_algorithm_param(algorithm, "N")
            Threads.@spawn optimize(problem.fn,
                                    algorithm.method,
                                    constraints = constraints,
                                    population = [Individual(rand(30)) for i in 1:population_size],
                                    options = Options(state_callback = state_callback))
        end

        @cstatic v = Cint(0) begin
            v = Int32(generation)
            @c CImGui.InputInt("Generation", &v)
        end
        @cstatic v = Cfloat(0.0) begin
            v = Float32(igd)
            @c CImGui.InputFloat("IGD", &v)
        end
        @cstatic v = Cfloat(0.0) begin
            v = Float32(hv)
            @c CImGui.InputFloat("HV", &v)
        end

        if runtime_image !== nothing
            if runtime_image_id !== nothing
                ImGui_ImplOpenGL3_DestroyImageTexture(runtime_image_id)
            end
            img_width, img_height = size(runtime_image)
            runtime_image_id = ImGui_ImplOpenGL3_CreateImageTexture(img_width, img_height)
            ImGui_ImplOpenGL3_UpdateImageTexture(runtime_image_id, image_uint8(runtime_image), img_width, img_height)
            CImGui.Image(Ptr{Cvoid}(runtime_image_id), (img_width, img_height))
        end

        if igd_image !== nothing && length(igds) > 0
            if igd_image_id !== nothing
                ImGui_ImplOpenGL3_DestroyImageTexture(igd_image_id)
            end
            img_width, img_height = size(igd_image)
            igd_image_id = ImGui_ImplOpenGL3_CreateImageTexture(img_width, img_height)
            ImGui_ImplOpenGL3_UpdateImageTexture(igd_image_id, image_uint8(igd_image), img_width, img_height)
            CImGui.Image(Ptr{Cvoid}(igd_image_id), (img_width, img_height))
        end

        if hv_image !== nothing && length(hvs) > 0
            if hv_image_id !== nothing
                ImGui_ImplOpenGL3_DestroyImageTexture(hv_image_id)
            end
            img_width, img_height = size(hv_image)
            hv_image_id = ImGui_ImplOpenGL3_CreateImageTexture(img_width, img_height)
            ImGui_ImplOpenGL3_UpdateImageTexture(hv_image_id, image_uint8(hv_image), img_width, img_height)
            CImGui.Image(Ptr{Cvoid}(hv_image_id), (img_width, img_height))
        end

        CImGui.End()
    end
end


text:
@testset "Right Division" begin
    mat = diagm([1,2,3])
    v = [2,3,4]
    tv = v.'

    @test (tv/mat)::TransposedVector ≈ [2/1  3/2  4/3]

    @test (v.'/mat)::TransposedVector ≈ [2/1  3/2  4/3]
    @test (v.'/mat.')::TransposedVector ≈ [2/1  3/2  4/3]
    @test (tv/mat.')::TransposedVector ≈ [2/1  3/2  4/3]

    @test (v'/mat)::TransposedVector ≈ [2/1  3/2  4/3]
    @test (v'/mat')::TransposedVector ≈ [2/1  3/2  4/3]
    @test (tv/mat')::TransposedVector ≈ [2/1  3/2  4/3]
end


text:
# Author: Ritchie Lee, ritchie.lee@sv.cmu.edu
# Date: 04/14/2015
# Implements Lincoln Lab scripted dynamics model

module LLADMImpl

export
addObserver,

initialize,
update,

initializeDynamicModel,
simulateDynamicModel,

LLADM,
LLADMState,
LLADMOutputState,
LLADMCommand

using AbstractDynamicModelImpl
using AbstractDynamicModelInterfaces
using CommonInterfaces
using RLESUtils, Observers

using Base.Test

import CommonInterfaces.addObserver
import CommonInterfaces.initialize
import CommonInterfaces.update
import AbstractDynamicModelInterfaces.initializeDynamicModel
import AbstractDynamicModelInterfaces.simulateDynamicModel

# Note: Angles are all in radians in this module unless otherwise indicated
# CASSATT_TIMING is false

type LLADMState

  t::Float64       #time (s)
  v::Float64       #airspeed (ft/s)
  N::Float64       #position north (ft)
  E::Float64       #position east (ft)
  h::Float64       #position altitude (ft)
  psi::Float64     #heading angle (rad), N is zero, clockwise is positive
  theta::Float64   #pitch angle (rad)
  phi::Float64     #bank angle (rad)
  a::Float64       #airspeed acceleration (ft/s^2)
end

type LLADMOutputState #mitcas state in scripted_dynamics.h

  t::Float64
  x::Float64
  y::Float64
  h::Float64
  vx::Float64   #DX in scripted_dynamics.h
  vy::Float64   #DY in scripted_dynamics.h
  vh::Float64   #DH in scripted_dynamics.h
  psi::Float64
  theta::Float64
  phi::Float64
end

LLADMOutputState() = LLADMOutputState(0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0)

type LLADMCommand

  t::Float64        #time of command (s)
  hdot::Float64     #commanded vertical rate (ft/s)
  psidot::Float64   #commanded turn rate (rad/s)
  a::Float64        #commanded acceleration (ft/s^2)
  dh_min::Float64          #dh_min of advisory being followed (ft/s)
  dh_max::Float64          #dh_max of advisory being followed (ft/s)
  target_rate::Float64     #target_rate of advisory being followed (ft/s)
  ddh::Float64             #ddh of advisory being followed (ft/s^2)

end

immutable LLADMConsts

  _dt::Float64           # minor time step
  _K::Float64            # integrator gain
  _g::Float64            # acceleration of gravity (ft/s^2)
  _qmax::Float64         # as in CASSATT (rad/s)
  _rmax::Float64         # as in CASSATT GA_psidotMAX = 1e6
  _phi_max::Float64
  _phidotmax::Float64
  _v_high::Float64       # airspeed limits
  _v_low::Float64
  LIMIT_DDH::Bool
  _min_vertical_rate::Float64
  _max_vertical_rate::Float64
  _min_vertical_accel::Float64
  _max_vertical_accel::Float64

end

type LLADM <: AbstractDynamicModel

  state::LLADMState
  update::Union{LLADMCommand, Void}
  timestep::Float64
  params::LLADMConsts
  output_state::LLADMOutputState

  observer::Observer

  function LLADM(; timestep::Float64 = 1.0, #major time step
                 _dt::Float64 = 0.1, #minor time step
                 _K::Float64 = 1.0,
                 _g::Float64 = 32.2,
                 _qmax::Float64 = deg2rad(3.0),
                 _rmax::Float64 = 1000000.0,
                 _phi_max::Float64 = deg2rad(75.0),
                 _phidotmax::Float64 = 0.524,
                 _v_high::Float64 = 100000.0,
                 _v_low::Float64 = 1.7,
                 LIMIT_DDH::Bool = true,
                 _min_vertical_rate::Float64 = -9999.0,
                 _max_vertical_rate::Float64 = 9999.0,
                 _min_vertical_accel::Float64 = -9.66,
                 _max_vertical_accel::Float64 = 9.66,
                 )

    obj = new()

    obj.update = nothing
    obj.timestep = timestep

    obj.params = LLADMConsts(_dt,
                               _K,
                               _g,
                               _qmax,
                               _rmax,
                               _phi_max,
                               _phidotmax,
                               _v_high,
                               _v_low,
                               LIMIT_DDH,
                               _min_vertical_rate,
                               _max_vertical_rate,
                               _min_vertical_accel,
                               _max_vertical_accel
                               )

    obj.output_state = LLADMOutputState()

    obj.observer = Observer()

    return obj
  end
end

addObserver(adm::LLADM, f::Function) = add_observer(adm.observer, f)
addObserver(adm::LLADM, tag::AbstractString, f::Function) = add_observer(adm.observer, tag, f)

initialize(adm::LLADM, state) = initialize(adm,convert(LLADMState, state))

function initialize(adm::LLADM, state::LLADMState)
  adm.state = state

  # Intialize the mitcas state
  # Note: x is E and y is N in scripted_dynamics.h
  adm.output_state.t     = state.t
  adm.output_state.h     = state.h
  adm.output_state.vh    = sin(state.theta) * state.v
  #adm.output_state.x     = state.E
  #adm.output_state.vx    = sin(state.psi) * cos(state.theta) * state.v
  #adm.output_state.y     = state.N
  #adm.output_state.vy    = cos(state.psi) * cos(state.theta) * state.v
  adm.output_state.psi   = state.psi
  adm.output_state.theta = state.theta
  adm.output_state.phi   = state.phi

  # RLEE: Override x and y to match LLCEM.  X is north
  adm.output_state.x     = state.N
  adm.output_state.vx    = cos(state.psi) * cos(state.theta) * state.v
  adm.output_state.y     = state.E
  adm.output_state.vy    = sin(state.psi) * cos(state.theta) * state.v

  adm.update = nothing

  return adm.output_state
end

# Mirrors state_update() from scripted_dynamics.h
function state_update(C::LLADMConsts, state::LLADMState, ctrl::LLADMCommand, hddcmd::Float64, dt::Float64)
  # Inputs:
  # C = constants
  # state = current state (adm.state)
  # ctrl_update = current command (update::LLADMCommand)
  # hddcmd = vertical acceleration to apply
  # dt = time step

  const s_theta = sin(state.theta)
  const c_theta = cos(state.theta)
  const t_theta = tan(state.theta)
  const s_phi = sin(state.phi)
  const c_phi = cos(state.phi)
  const s_psi = sin(state.psi)
  const c_psi = cos(state.psi)

  # save v and saturated v
  const v = state.v
  const v_saturated = max(v,1.0)

  # compute and saturate q
  q = saturate(1.0 / (v_saturated * c_phi) * (hddcmd / c_theta + C._g * c_theta * s_phi * s_phi - ctrl.a * t_theta), -C._qmax, C._qmax)

  # compute and saturate r
  r = saturate(C._g * s_phi * c_theta / v_saturated, -C._rmax, C._rmax)

  # sometimes the aircraft bank angle is so large that it would prevent reaching desired h_dot,
  # so we dynamically saturate the bank angle so that h_dot is always acheived
  # compute phimax (in the Plopt/compute q block)
  hdd_cmd_phi = min(hddcmd, v_saturated * C._qmax * c_phi * c_theta)

  # calculate discriminant
  sqrt_arg = v_saturated * v_saturated * C._qmax * C._qmax - 4 * C._g * ctrl.a * s_theta + 4 * C._g * hdd_cmd_phi + 4 * C._g * C._g * c_theta * c_theta
  phi_max_2 = 0

  if sqrt_arg < 0
    phi_max_2 = 10000
  else
    # calculate cos(phi)
    cphi1 = (sqrt(sqrt_arg) - v_saturated * C._qmax) / (2 * C._g * c_theta)

    if abs(cphi1) < 1
      phi_max_2 = acos(cphi1) * 0.98 # add a small buffer to prevent jittering
    end
  end

  phimax = min(C._phi_max, phi_max_2)

  # compute and saturate p
  psidot_if_no_change = (q * s_phi + r * c_phi) / c_theta
  dpsidot = ctrl.psidot - psidot_if_no_change
  p = saturate(20.0 * dpsidot, -C._phidotmax, C._phidotmax) # 20 is user tuned parameter from CASSATT

  # limit max bank angle
  if (state.phi + p * dt > phimax)

    p = (phimax - state.phi) / dt
  end

  if (state.phi + p * dt < -phimax)

    p = (-phimax - state.phi) / dt
  end

  state.v += dt * (ctrl.a)
  state.N += dt * (v * c_theta * c_psi)
  state.E += dt * (v * c_theta * s_psi)
  state.h += dt * (v * s_theta)
  state.psi += dt * (q * s_phi / c_theta + r * c_phi / c_theta)
  state.psi = to_plusminus_pi(state.psi) #RLEE added: mod back to [-pi, pi]
  state.theta += dt * (q * c_phi - r * s_phi)
  state.theta = to_plusminus_pi(state.theta) #RLEE added: mod back to [-pi, pi]
  state.phi += dt * (p + q * s_phi * t_theta + r * c_phi * t_theta)
  state.phi = to_plusminus_pi(state.phi) #RLEE added: mode back to [-pi, pi]

  state.v = saturate(state.v, 1.7, C._v_high - 0.000001)
  state.t += dt

end

update(adm::LLADM, ctrl) = update(adm,convert(LLADMCommand, ctrl))

function update(adm::LLADM, ctrl::LLADMCommand)

  state = adm.state
  C = adm.params #constants

  t_start = state.t

  for i = 1:Int64(round(adm.timestep / C._dt))
    substep(state, ctrl, C)
  end

  #snap so small errors don't accumulate
  state.t = snap(state.t, t_start + adm.timestep)

  #make sure timing is in sync
  @test state.t == t_start + adm.timestep

  # calculate the new mitcas state vars
  # Note: x is E and y is N in scripted_dynamics.h
  adm.output_state.t     = state.t
  adm.output_state.h     = state.h
  adm.output_state.vh    = sin(state.theta) * state.v
  #adm.output_state.x     = state.E
  #adm.output_state.vx    = sin(state.psi) * cos(state.theta) * state.v
  #adm.output_state.y     = state.N
  #adm.output_state.vy    = cos(state.psi) * cos(state.theta) * state.v
  adm.output_state.psi   = state.psi
  adm.output_state.theta = state.theta
  adm.output_state.phi   = state.phi

  # RLEE: Override x and y to match LLCEM.  X is north
  adm.output_state.x     = state.N
  adm.output_state.vx    = cos(state.psi) * cos(state.theta) * state.v
  adm.output_state.y     = state.E
  adm.output_state.vy    = sin(state.psi) * cos(state.theta) * state.v

  return adm.output_state
end

function substep(state::LLADMState, ctrl::LLADMCommand, C::LLADMConsts)

  dh = state.v * sin(state.theta)

  dummy, ddh_script = response(dh, (ctrl.hdot - dh) / C._dt,
                                       ctrl.dh_min, ctrl.dh_max,
                                       C._min_vertical_rate, C._max_vertical_rate, C._dt)

  dummy, ddh_cmd = response(dh, ctrl.dh_min, ctrl.dh_max, 0.0, 0.0, ctrl.ddh,
                                 C._min_vertical_rate, C._max_vertical_rate, C._dt)

  ddh2apply = resolve_TCAS_and_script(ddh_cmd,ddh_script)

  #debug
  #println("ddh_script=", ddh_script, ", ddh_cmd=", ddh_cmd, ", ddh2apply=", ddh2apply)

  if C.LIMIT_DDH
    ddh2apply = saturate(ddh2apply, C._min_vertical_accel, C._max_vertical_accel)
  end

  state_update(C, state, ctrl, ddh2apply, C._dt)
end

function saturate(x::Float64, limit_min::Float64, limit_max::Float64)
  @test limit_min <= limit_max

  return min(max(x,limit_min),limit_max)
end

function response(dh::Float64, ddh_in::Float64, dh_cmd_min::Float64, dh_cmd_max::Float64,
                  min_dh::Float64, max_dh::Float64, dt::Float64)

  # default value for ddh_out
	ddh_out = ddh_in

	# new projected vertical rate
	dh_proj = dh + ddh_out * dt

	# whether or not inside green zone
	dh_inside_green = saturate(dh, dh_cmd_min, dh_cmd_max) == dh

	if (dh > dh_cmd_max && dh_proj < dh_cmd_min) || (dh_inside_green && dh_proj > dh_cmd_max)

		# OUTSIDE GREEN ZONE
		# cross max edge from above and go below min edge => hit max edge
		# OR
		# INSIDE GREEN ZONE
		# cross max edge from below and go above max edge => hit max edge

		dh_proj = dh_cmd_max

	elseif (dh < dh_cmd_min && dh_proj > dh_cmd_max) || (dh_inside_green && dh_proj < dh_cmd_min)

		# OUTSIDE GREEN ZONE
		# cross min edge from below and go above max edge => hit min edge
		# OR
		# INSIDE GREEN ZONE
		# cross min edge from above and go below min edge => hit min edge

		dh_proj = dh_cmd_min

	elseif dh_proj < min_dh

		# cross max descent rate => hit max descent rate
		dh_proj = min_dh

	elseif dh_proj > max_dh

		# cross max climb rate => hit max climb rate
		dh_proj = max_dh

  end

	# recompute ddh
	ddh_out = dt != 0.0 ? (dh_proj - dh) / dt : 0.0

	# compute next dh
	dh_next = snap(dh + ddh_out * dt, dh_cmd_min)
	dh_next = snap(dh_next, dh_cmd_max)
	dh_next = snap(dh_next, min_dh)
	dh_next = snap(dh_next, max_dh)

  return dh_next, ddh_out
end

function response(dh::Float64, dh_cmd_min::Float64, dh_cmd_max::Float64, ddh_no_response::Float64,
                  ddh_response::Float64, cmd_ddh::Float64, min_dh::Float64, max_dh::Float64, dt::Float64)

  ddh_in = acceleration_to_apply(dh, ddh_no_response, ddh_response, dh_cmd_min, dh_cmd_max, cmd_ddh)

  return response(dh, ddh_in, dh_cmd_min, dh_cmd_max, min_dh, max_dh, dt)
end

function acceleration_to_apply(dh::Float64, ddh_no_response::Float64, ddh_response::Float64,
                               dh_cmd_min::Float64, dh_cmd_max::Float64, cmd_ddh::Float64)

  ddh = ddh_no_response

  # resolution advisory handling
	if (saturate(dh, dh_cmd_min, dh_cmd_max) != dh)

		# commanded acceleration
		ddh = cmd_ddh

		# determine direction
		if (dh > dh_cmd_max)
      ddh = -ddh
    end

		# add response noise
		ddh += ddh_response;
  end

	return ddh
end

function resolve_TCAS_and_script(ddh_TCAS::Float64, ddh_script::Float64)

  opposite_senses = ddh_script * ddh_TCAS < 0
  same_senses = !opposite_senses

  cmd_greater = abs(ddh_TCAS) > abs(ddh_script)

  return (ddh_TCAS != 0.0 && ((same_senses && cmd_greater) || opposite_senses)) ? ddh_TCAS : ddh_script
end

snap(x, y) = approxeq(x, y, 1e-8) ? y : x

approxeq(x, y, eps) = (x == y) ? true : abs(x-y) < eps

#mods x to the range [-b, b]
function to_plusminus_b(x::AbstractFloat, b::AbstractFloat)
  z = mod(x, 2 * b)
  return (z > b) ? (z - 2 * b) : z
end

to_plusminus_pi(x::AbstractFloat) = to_plusminus_b(x, Float64(pi))

end


text:
using Revise

using DynamicBoundspODEsPILMS, LinearAlgebra, IntervalArithmetic, StaticArrays, TaylorSeries,
      TaylorIntegration, ForwardDiff, McCormick, BenchmarkTools, DocStringExtensions

using DiffResults: JacobianResult, MutableDiffResult
using ForwardDiff: Partials, JacobianConfig, vector_mode_dual_eval, value, vector_mode_jacobian!

import Base.copyto!

function f!(dx,x,p,t)
    dx[1] = x[1]^2 + p[2]
    dx[2] = x[2] + p[1]^2
    nothing
end
np = 2
nx = 2
k = 3
x = [0.1; 1.0]
p = [0.2; 0.1]
jtf! = DynamicBoundspODEsPILMS.JacTaylorFunctor!(f!, nx, np, k, Interval{Float64}(0.0), 0.0)

xIntv = Interval{Float64}.(x)
pIntv = Interval{Float64}.(p)
yIntv = [xIntv; pIntv]
DynamicBoundspODEsPILMS.jacobian_taylor_coeffs!(jtf!, yIntv)

jac = JacobianResult(jtf!.out, yIntv).derivs[1]
tjac = zeros(Interval{Float64}, 4, 8)
Jx = Matrix{Interval{Float64}}[zeros(Interval{Float64},2,2) for i in 1:4]
Jp = Matrix{Interval{Float64}}[zeros(Interval{Float64},2,2) for i in 1:4]
DynamicBoundspODEsPILMS.extract_JxJp!(Jx, Jp, jtf!.result, tjac, nx, np, k)

itf! = DynamicBoundspODEsPILMS.TaylorFunctor!(f!, nx, np, k, zero(Interval{Float64}), zero(Float64))
outIntv = zeros(Interval{Float64},8)
itf!(outIntv, yIntv)

y = [x; p]
rtf!  = DynamicBoundspODEsPILMS.TaylorFunctor!(f!, nx, np, k, zero(Float64), zero(Float64))
out = zeros(8)
rtf!(out, y)

coeff_out = zeros(Interval{Float64},2,4)
DynamicBoundspODEsPILMS.coeff_to_matrix!(coeff_out, jtf!.out, nx, k)

#=
hⱼ = 0.001
hmin = 0.00001
function euf!(out, x, p, t)
    out[1,1] = -x[1]
    nothing
end
eufY = [Interval{Float64}(0.5,1.5); Interval(0.0)]
itf_exist_unique! = DynamicBoundspODEsPILMS.TaylorFunctor!(euf!, 1, 1, k, zero(Interval{Float64}), zero(Float64))
jtf_exist_unique! = DynamicBoundspODEsPILMS.JacTaylorFunctor!(euf!, 1, 1, k, Interval{Float64}(0.0), 0.0)
DynamicBoundspODEsPILMS.jacobian_taylor_coeffs!(jtf_exist_unique!, eufY)
coeff_out = zeros(Interval{Float64},1,k)
DynamicBoundspODEsPILMS.coeff_to_matrix!(coeff_out, jtf!.out, 1, k)
Jx = Matrix{Interval{Float64}}[zeros(Interval{Float64},1,1) for i in 1:4]
Jp = Matrix{Interval{Float64}}[zeros(Interval{Float64},1,1) for i in 1:4]
tjac = zeros(Interval{Float64}, 2, 4)
outIntv_exist_unique! = zeros(Interval{Float64},4)
itf_exist_unique!(outIntv_exist_unique!, eufY)
coeff_out_exist_unique! = zeros(Interval{Float64},1,k+1)
DynamicBoundspODEsPILMS.coeff_to_matrix!(coeff_out_exist_unique!, outIntv_exist_unique!, 1, k)
DynamicBoundspODEsPILMS.extract_JxJp!(Jx, Jp, jtf_exist_unique!.result, tjac, 1, 1, k)
unique_result = DynamicBoundspODEsPILMS.UniquenessResult(1,1)
DynamicBoundspODEsPILMS.existence_uniqueness!(unique_result, itf_exist_unique!,
                                              eufY, hⱼ, hmin, coeff_out_exist_unique!, Jx, Jp)

bool1 = unique_result.step == 0.001
bool2 = unique_result.confirmed
bool3a = isapprox(unique_result.Ỹⱼ[1].lo, -1.50001E-6, atol=1E-10)
bool3b = isapprox(unique_result.Ỹⱼ[1].hi, 0.00150001, atol=1E-6)
bool4 = isapprox(unique_result.f̃k[1].lo, -7.50001E-16, atol = 1E-19)
bool5 = isapprox(unique_result.f̃k[1].hi, 7.50001E-13, atol = 1E-16)
=#

function fplohn!(out, x, p, t)
    out[1] = x[1]
    out[2] = -x[2]
    nothing
end
np = 1
nx = 2
k = 3
lf! = DynamicBoundspODEsPILMS.LohnersFunctor!(fplohn!, nx, np, k, zero(Interval{Float64}), zero(Float64))
hⱼ = 0.001
Ỹⱼ = [Interval(0.1, 5.1); Interval(0.1, 8.9); Interval(0.1, 8.9)]
Yⱼ = [Interval(0.1, 5.1); Interval(0.1, 8.9); Interval(0.1, 8.9)]
A = DynamicBoundspODEsPILMS.qr_stack(nx, 2)
yⱼ = mid.(Yⱼ)
Δⱼ = Yⱼ[1:2] - yⱼ[1:2]
lf!(hⱼ, Ỹⱼ, Yⱼ, A, yⱼ, Δⱼ)

#=
jetcoeffs!(zqwa, zqwb, zqwc, zqwd, zqwe, zqwr, p)
y = Interval{Float64}.([x; p])
out = g.out
cfg = ForwardDiff.JacobianConfig(nothing, out, y)

# extact is good... actual jacobians look odd...
#tv, xv = validated_integration(f!, Interval{Float64}.([3.0, 3.0]), 0.0, 0.3, 4, 1.0e-20, maxsteps=100 )
Q = [Yⱼ; P]
#@btime jacobianfunctor($outIntv, $yInterval)

d = g
zqwa = d.g!
zqwb = d.t
zqwc = d.xtaylor
zqwd = d.xout
zqwe = d.xaux
zqwr = d.taux
#@btime jetcoeffs!($zqwa, $zqwb, $zqwc, $zqwd, $zqwe, $zqwr, $s, $p)
#@code_warntype jetcoeffs!(zqwa, zqwb, zqwc, zqwd, zqwe, zqwr, p)

Jx = Matrix{Interval{Float64}}[zeros(Interval{Float64},2,2) for i in 1:4]
Jp = Matrix{Interval{Float64}}[zeros(Interval{Float64},2,2) for i in 1:4]
Jxsto = zeros(Interval{Float64},2,2)
Jpsto = zeros(Interval{Float64},2,2)

Yⱼ = [Interval{Float64}(-10.0, 20.0); Interval{Float64}(-10.0, 20.0)]
P = [Interval{Float64}(2.0, 3.0); Interval{Float64}(2.0, 3.0)]
Ycat = [Yⱼ; P]
yⱼ = mid.(Yⱼ)
Δⱼ = Yⱼ - yⱼ
At = zeros(2,2) + I
#Aⱼ =  DynamicBoundspODEsPILMS.QRDenseStorage(nx)
#Aⱼ₊₁ =  DynamicBoundspODEsPILMS.QRDenseStorage(nx)
A =  DynamicBoundspODEsPILMS.QRStack(nx, 2)
dtf = g
hⱼ = 0.001
yjcat = vcat(yⱼ,p)
# TODO: Remember rP is computed outside iteration and stored to JacTaylorFunctor
plohners = DynamicBoundspODEsPILMS.parametric_lohners!(itf!, rtf!, dtf, hⱼ, Ycat, Ycat,
                                                       A, yjcat, Δⱼ)

@btime DynamicBoundspODEsPILMS.parametric_lohners!($itf!, $rtf!, $dtf, $hⱼ, $Ycat, $Ycat,
                                                    $A, $yjcat, $Δⱼ)
=#


text:
# TODO document
reduce_axes(old_axes::Tuple{Vararg{Any,N}}, new_axes::Tuple, dims::Colon) where {N} = ()
function reduce_axes(old_axes::Tuple{Vararg{Any,N}}, new_axes::Tuple, dims) where {N}
    ntuple(Val(N)) do i
        if i in dims
            StaticRanges.shrink_last(getfield(old_axes, i), getfield(new_axes, i))
        else
            unsafe_reconstruct(getfield(old_axes, i), getfield(new_axes, i))
        end
    end
end

#=
We need to assign new indices to axes of `A` but `reshape` may have changed the
size of any axis
=#
@inline function reshape_axes(axs::Tuple, inds::Tuple{Vararg{Any,N}}) where {N}
    return map((a, i) -> resize_last(a, i), axs, inds)
end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:3] <= 1)

                     @variable(m, -8.232182974914569 <= q <= 6.070970235703473)

                     add_NL_constraint(m, :(1/(1 + exp(-(-0.17705543456567296 + 0.5634748900178774*1/(1 + exp(-(-0.9628319588335761 + 0.5826659544327617*$(x[1]) + -0.6099922924008303*$(x[2]) + 0.9101904633653426*$(x[3])))) + 0.6106029608527352*1/(1 + exp(-(0.06629760484384972 + -0.7340864295918403*$(x[1]) + 0.2309908240081593*$(x[2]) + 0.7867070707759196*$(x[3])))) + 0.5384388086098779*1/(1 + exp(-(-0.07370832879656497 + -0.32073040453400337*$(x[1]) + -0.07649911788086428*$(x[2]) + 0.03499506906851524*$(x[3]))))))) + 1/(1 + exp(-(0.08053767290922664 + 0.911477397877237*1/(1 + exp(-(-0.9628319588335761 + 0.5826659544327617*$(x[1]) + -0.6099922924008303*$(x[2]) + 0.9101904633653426*$(x[3])))) + 0.7134526118957916*1/(1 + exp(-(0.06629760484384972 + -0.7340864295918403*$(x[1]) + 0.2309908240081593*$(x[2]) + 0.7867070707759196*$(x[3])))) + -0.34097520379571256*1/(1 + exp(-(-0.07370832879656497 + -0.32073040453400337*$(x[1]) + -0.07649911788086428*$(x[2]) + 0.03499506906851524*$(x[3]))))))) + 1/(1 + exp(-(0.5992621446321471 + 0.15304126986204158*1/(1 + exp(-(-0.9628319588335761 + 0.5826659544327617*$(x[1]) + -0.6099922924008303*$(x[2]) + 0.9101904633653426*$(x[3])))) + -0.37134199071135*1/(1 + exp(-(0.06629760484384972 + -0.7340864295918403*$(x[1]) + 0.2309908240081593*$(x[2]) + 0.7867070707759196*$(x[3])))) + 0.8747239793986465*1/(1 + exp(-(-0.07370832879656497 + -0.32073040453400337*$(x[1]) + -0.07649911788086428*$(x[2]) + 0.03499506906851524*$(x[3]))))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
using DifferentialEquations, Test, Sundials

using DiffEqProblemLibrary.DAEProblemLibrary: importdaeproblems; importdaeproblems()
import DiffEqProblemLibrary.DAEProblemLibrary: prob_dae_resrob

prob = prob_dae_resrob
sol = solve(prob)

@test typeof(sol.alg) <: Sundials.IDA


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule Dex_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("Dex")
JLLWrappers.@generate_main_file("Dex", UUID("f1ef5e10-671a-599f-ac25-3c68827556ba"))
end  # module Dex_jll


text:
using JudyDicts

import Base.getindex

const NLOOPS = isempty(ARGS) ? 1000000 : int(ARGS[1])
const kk = [ bytestring(repeat(string(idx), 10)) for idx in 1:NLOOPS ]

macro nogc_elapsed(ex)
    quote
        gc_disable()
        local t0 = time_ns()
        local val = $(esc(ex))
        local ret = (time_ns()-t0)/1e9
        gc_enable(); gc()
        ret
    end
end

macro bmsetstrobjloop(nloops, arr)
    quote
        local i::Int
        @nogc_elapsed for i in 1:$(nloops)
            $(arr)[kk[i]] = kk[i]
        end
    end
end

macro bmgetstrobjloop(nloops, arr)
    quote
        local x::String = ""
        @nogc_elapsed for i in 1:$(nloops)
            x = $(arr)[kk[i]]
        end
    end
end

function compare_str_obj()
    println("items ", NLOOPS, " compare: JudyDict{String, ASCIIString} vs. Dict{String, ASCIIString}")
    local dict_ins, dict_access, ja_ins, ja_access

    let ja = JudyDict{String, ASCIIString}()
        ja_ins = @bmsetstrobjloop NLOOPS ja
        ja_access = @bmgetstrobjloop NLOOPS ja
    end
    gc()

    let d = Dict{String, ASCIIString}()
        dict_ins = @bmsetstrobjloop NLOOPS d
        dict_access = @bmgetstrobjloop NLOOPS d
    end
    gc()
    #println("dict ins:", dict_ins, " access:", dict_access)

    println("set => dict: ", dict_ins, ", judy: ", ja_ins);
    println("get => dict: ", dict_access, ", judy: ", ja_access);
end

macro bmgetstrloop(nloops, arr)
    quote
        local x::Int = 0
        @nogc_elapsed for i in 1:$(nloops)
            x += $(arr)[kk[i]]
        end
    end
end

macro bmsetstrloop(nloops, arr)
    quote
        @nogc_elapsed for i in 1:$(nloops)
            $(arr)[kk[i]] = i
        end
    end
end

function compare_str()
    println("items ", NLOOPS, " compare: JudyDict{String, Int} vs. Dict{String, Int64}")
    local dict_ins, dict_access, ja_ins, ja_access

    let ja = JudyDict{String, Int}()
        ja_ins = @bmsetstrloop NLOOPS ja
        ja_access = @bmgetstrloop NLOOPS ja
    end
    gc()

    let d = Dict{String, Int64}()
        dict_ins = @bmsetstrloop NLOOPS d
        dict_access = @bmgetstrloop NLOOPS d
    end
    gc()

    println("set => dict: ", dict_ins, ", judy: ", ja_ins);
    println("get => dict: ", dict_access, ", judy: ", ja_access);
end

macro bmsetintloop(nloops, arr)
    quote
        @nogc_elapsed for i in 1:$(nloops); $(arr)[i] = i; end
    end
end

macro bmgetintloop(nloops, arr)
    quote
        local x::Int = 0
        @nogc_elapsed for i in 1:$(nloops); x += $(arr)[i]; end
    end
end

function compare_int64()
    println("items ", NLOOPS, " compare: JudyDict{Int, Int} vs. Dict{Int64, Int64}")
    local dict_ins, dict_access, ja_ins, ja_access

    let d = Dict{Int64, Int64}()
        dict_ins = @bmsetintloop NLOOPS d
        dict_access = @bmgetintloop NLOOPS d
    end
    gc()

    let ja = JudyDict{Int, Int}()
        ja_ins = @bmsetintloop NLOOPS ja
        ja_access = @bmgetintloop NLOOPS ja
    end
    gc()

    println("set => dict: ", dict_ins, ", judy: ", ja_ins);
    println("get => dict: ", dict_access, ", judy: ", ja_access);
end

compare_int64()
println("")
compare_str()
println("")
compare_str_obj()


text:
T=3
tsname="joDWT"
@testset "$tsname" begin
for t=1:T # start test loop
    m=4^t
    v1=rand(m)
    v2=rand(m,m)
    vv2=vec(v2)

    wt=wavelet(WT.haar)
    A1=joDWT(m,wt)
    A2=joDWT(m,m,wt)
    verbose && println("$tsname ($m,$m) - planned")
    @testset "$m x $m" begin
        @test isadjoint(joDWT(m))[1]
        @test islinear(joDWT(m))[1]
        @test norm(A1*v1-dwt(v1,wt))<joTol
        @test norm(A1\v1-idwt(v1,wt))<joTol
        @test norm(adjoint(A1)*v1-idwt(v1,wt))<joTol
        @test norm((adjoint(A1)*A1)*v1-v1)<joTol
        @test isadjoint(joDWT(m,m))[1]
        @test islinear(joDWT(m,m))[1]
        @test norm(A2*vv2-vec(dwt(v2,wt)))<joTol
        @test norm(A2\vv2-vec(idwt(v2,wt)))<joTol
        @test norm(adjoint(A2)*vv2-vec(idwt(v2,wt)))<joTol
        @test norm((adjoint(A2)*A2)*vv2-vv2)<joTol
    end

end # end test loop
end


text:
module day19

using OffsetArrays  # rules ids are 0-indexed

include("inputs.jl")

struct Rule
    id::Int
    literal::Union{Char, Nothing}
    subrules::Union{Array{Array{Int,1},1}, Nothing}
end

ex = inputs.ex_d19();
input = inputs.input_d19();

function parse_input(input)
    rules_raw, msgs_raw = split(input, "\n\n", keepempty=false)
    rules = []

    for line in split(rules_raw, "\n")
        id, rest = split(line, ": ")
        id = parse(Int, id)
        literal, subrules = nothing, nothing

        if contains(rest, "\"")
            literal=rest[2]
        else
            parts = contains(rest, "|") ? split(rest, " | ") : parts = [rest]
            subrules = [[parse(Int,n) for n in split(part)] for part in parts]
        end
        push!(rules, Rule(id, literal, subrules))
    end

    rules = OffsetArray(sort!(rules, by=x->x.id), 0:length(rules)-1)
    msgs = split(msgs_raw, "\n")
    return rules, msgs
end

"""
Returns true if msg matches rule 0.
"""
function check(msg, rules)

    # Queue of pairs (remaining string, remaining rules)
    q = Vector([(msg, [0])])

    while ~isempty(q)
        msg, rule_ids = popfirst!(q)

        # consumed whole msg and all rules simultaneously, so that's a match
        if isempty(msg) && isempty(rule_ids)
            return true
        # consumed msg or rules but not both. Dead end.
        elseif isempty(msg) || isempty(rule_ids)
            continue
        end

        # have both msg and rule_ids. So try first rule in our ids
        rule = rules[rule_ids[1]]
        rule_ids = rule_ids[2:end]

        # first rule is literal, so if it matches the first character,
        # then add the rest of the msg and the rest of the rules to the queue
        if ~isnothing(rule.literal)
            if msg[1] == rule.literal
                push!(q, (msg[2:end], rule_ids))
            else
                continue
            end
        else
            for subrule_ids in rule.subrules
                push!(q, (msg, vcat(subrule_ids, rule_ids)))
            end
        end
    end

    return false
end

function part1(input)
    rules, msgs = parse_input(input)
    return count([check(msg, rules) for msg in msgs])
end

@assert part1(ex) == 2
@info "Day 19, Part 1 answer: $(part1(input))"

function part2(input)
    input = replace(input, "8: 42"=>"8: 42 | 42 8") |> x->
            replace(x, "11: 42 31"=>"11: 42 31 | 42 11 31")
    rules, msgs = parse_input(input)
    return count([check(msg, rules) for msg in msgs])
end

@info "Day 19, Part 2 answer: $(part2(input))"

end


text:
import Base: similar, ∈, ⊆, ∩, ∪, real, imag
import Random: rand
import IntervalArithmetic: ±, inf, sup, mid, diam, radius, midpoint_radius, hull

"""
    AbstractIntervalMatrix{IT} <: AbstractMatrix{IT}

Abstract supertype for interval matrix types.
"""
abstract type AbstractIntervalMatrix{IT} <: AbstractMatrix{IT} end

"""
    IntervalMatrix{T, IT, MT<:AbstractMatrix{IT}} <: AbstractIntervalMatrix{IT}

An interval matrix i.e. a matrix whose coefficients are intervals. This type is
parametrized in the number field, the interval type, and the matrix type.

### Fields

- `mat` -- matrix whose entries are intervals

### Examples

```jldoctest
julia> A = IntervalMatrix([-0.9±0.1 0±0; 0±0 -0.9±0.1])
2×2 IntervalMatrix{Float64, Interval{Float64}, Matrix{Interval{Float64}}}:
 [-1, -0.799999]   [0, 0]
  [0, 0]          [-1, -0.799999]
```

An interval matrix proportional to the identity matrix can be built using the
`UniformScaling` operator from the standard library `LinearAlgebra`. For example,

```jldoctest interval_uniform_scaling
julia> using LinearAlgebra

julia> IntervalMatrix(Interval(1)*I, 2)
2×2 IntervalMatrix{Float64, Interval{Float64}, Matrix{Interval{Float64}}}:
 [1, 1]  [0, 0]
 [0, 0]  [1, 1]
```
The number of columns can be specified as a third argument, creating a rectangular
``m × n`` matrix such that only the entries in the main diagonal,
``(1, 1), (2, 2), …,  (k, k)`` are specified, where ``k = \\min(m, n)``:

```jldoctest interval_uniform_scaling
julia> IntervalMatrix(Interval(-1, 1)*I, 2, 3)
2×3 IntervalMatrix{Float64, Interval{Float64}, Matrix{Interval{Float64}}}:
 [-1, 1]   [0, 0]  [0, 0]
  [0, 0]  [-1, 1]  [0, 0]

julia> IntervalMatrix(Interval(-1, 1)*I, 3, 2)
3×2 IntervalMatrix{Float64, Interval{Float64}, Matrix{Interval{Float64}}}:
 [-1, 1]   [0, 0]
  [0, 0]  [-1, 1]
  [0, 0]   [0, 0]
```

An uninitialized interval matrix can be constructed using `undef`:

```jldoctest undef_test
julia> m = IntervalMatrix{Float64}(undef, 2, 2);

julia> typeof(m)
IntervalMatrix{Float64, Interval{Float64}, Matrix{Interval{Float64}}}
```
Note that this constructor implicitly uses a dense matrix, `Matrix{Float64}`,
as the matrix (`mat`) field in the new interval matrix.
"""
struct IntervalMatrix{T, IT, MT<:AbstractMatrix{IT}} <: AbstractIntervalMatrix{IT}
    mat::MT
end

import Base:size, IndexStyle, getindex, setindex!, copy

IndexStyle(::Type{<:IntervalMatrix}) = IndexLinear()
size(M::IntervalMatrix) = size(M.mat)
getindex(M::IntervalMatrix, i::Int) = getindex(M.mat, i)
setindex!(M::IntervalMatrix, X, inds...) = setindex!(M.mat, X, inds...)
copy(M::IntervalMatrix) = IntervalMatrix(copy(M.mat))

function IntervalMatrix(A::MT) where {T, IT<:AbstractInterval{T}, MT<:AbstractMatrix{IT}}
    return IntervalMatrix{T, IT, MT}(A)
end

function IntervalMatrix(A::MT) where {T, IT<:Complex{Interval{T}}, MT <:AbstractMatrix{IT}}
    return IntervalMatrix{T, IT, MT}(A)
end

# constructor from uniform scaling
function IntervalMatrix(αI::UniformScaling{Interval{T}}, m::Integer, n::Integer=m) where {T}
    return IntervalMatrix(Matrix(αI, m, n))
end

# undef initializer, eg. IntervalMatrix{Float64}(undef, 2, 2)
function IntervalMatrix{T}(u::UndefInitializer, m::Integer, n::Integer=m) where {T}
    mat = Matrix{Interval{T}}(undef, m, n)
    return IntervalMatrix(mat)
end

# similar initializer
function similar(M::IntervalMatrix)
    return IntervalMatrix(similar(M.mat))
end

# constructor from a scalar matrix
function IntervalMatrix(A::AbstractMatrix{T}) where {T<:Number}
    return IntervalMatrix(map(Interval, A))
end

"""
    IntervalMatrix(A::MT, B::MT) where {T, MT<:AbstractMatrix{T}}

Return an interval matrix such that the lower and upper bounds of the intervals
are given by the matrices `A` and `B` respectively.

### Input

- `A` -- lower bound matrix
- `B` -- upper bound matrix

### Output

An interval matrix `M` such that `M[i, j]` corresponds to the interval whose lower bound
is `A[i, j]` and whose upper bound is `B[i, j]`, for each `i` and `j`. That is,
``M_{ij} = [A_{ij}, B_{ij}]``.

### Notes

The upper bound should be bigger or equal than the lower bound,
i.e. `B[i, j] ≥ A[i, j]` for each `i` and `j`.

### Examples

```jldoctest
julia> IntervalMatrix([1 2; 3 4], [1 2; 4 5])
2×2 IntervalMatrix{Float64, Interval{Float64}, Matrix{Interval{Float64}}}:
 [1, 1]  [2, 2]
 [3, 4]  [4, 5]
```
"""
function IntervalMatrix(A::MT, B::MT) where {T, MT<:AbstractMatrix{T}}
    size(A) == size(B) || throw(ArgumentError("the sizes of the lower and upper bound " *
                                "matrices should match, but they are $(size(A)) " *
                                "and $(size(B)) respectively"))

    return IntervalMatrix(map((x, y) -> Interval(x, y), A, B))
end

"""
    ±(C::MT, S::MT) where {T, MT<:AbstractMatrix{T}}

Return an interval matrix such that the center and radius of the intervals
is given by the matrices `C` and `S` respectively.

### Input

- `C` -- center matrix
- `S` -- radii matrix

### Output

An interval matrix `M` such that `M[i, j]` corresponds to the interval whose center
is `C[i, j]` and whose radius is `S[i, j]`, for each `i` and `j`. That is,
``M = C + [-S, S]``.

### Notes

The radii matrix should be nonnegative, i.e. `S[i, j] ≥ 0` for each `i` and `j`.

### Examples

```jldoctest
julia> [1 2; 3 4] ± [1 2; 4 5]
2×2 IntervalMatrix{Float64, Interval{Float64}, Matrix{Interval{Float64}}}:
  [0, 2]   [0, 4]
 [-1, 7]  [-1, 9]
```
"""
function ±(C::MT, S::MT) where {T, MT<:AbstractMatrix{T}}
    size(C) == size(S) || throw(ArgumentError("the sizes of the center matrix and the " *
                                "radii matrix should match, but they are $(size(C)) " *
                                "and $(size(S)) respectively"))

    return IntervalMatrix(map((x, y) -> x ± y, C, S))
end


text:
# Author: Vincent Arel-Bundock (w/ changes by JMW)
# Contact: varel@umich.edu
# License: BSD-3
# Original code from the python statsmodels project
# https://github.com/statsmodels/statsmodels

using DataFrames

module QuantileRegression
    export qreg

    using DataFrames, Distributions

    type QRegModel
        beta::Vector{Float64}
        vcov::Matrix{Float64}
        mf::ModelFrame
    end

    function bandwidth_hall_sheather(n::Integer, q::Real, alpha::Real)
        z = quantile(Normal(), q)
        num = 1.5 * pdf(Normal(), z)^2
        den = 2.0 * z^2 + 1
        h = n^(-1/3) * quantile(Normal(), (1. - alpha / 2.))^(2./3) * (num / den)^(1./3)
        return h
    end

    # operate in-place to save on memory allocation
    function kernel_epanechnikov!(u::Vector, x::Vector)
        for i in 1:length(x)
            if abs(x[i]) <= 1
                # The @inbounds macro tells the compiler it can remove
                # bounds-checking for array indexing
                @inbounds u[i] = 0.75 * (1 - x[i]^2)
            else
                @inbounds u[i] = 0
            end
        end
        return
    end

    function kernel_epanechnikov(x::Vector)
        u = similar(x)
        kernel_epanechnikov!(u, x)
        return u
    end

    function maxdiff(x::Array, y::Array)
        m = typemin(Float64)
        for i in length(x)
            t = abs(x[i] - y[i])
            if t > m
                m = t
            end
        end
        return m
    end

    function qreg_coef(y::Vector, X::Matrix, q::Real = 0.5;
                       tol::Real = 1e-12, max_iter::Integer = 1_000,
                       threshold::Real = 1e-5)
        n, p = size(X)
        xstar = copy(X)
        diff = Inf

        beta0 = Array(Float64, p)
        beta = Array(Float64, p)
        xtx = Array(Float64, p, p)
        xty = Array(Float64, p)
        xbeta = Array(Float64, n)
        resid = Array(Float64, n)

        for itr in 1:max_iter
            if diff > tol
                copy!(beta0, beta)

                At_mul_B(xtx, xstar, X)
                At_mul_B(xty, xstar, y)
                beta = xtx \ xty
                A_mul_B(xbeta, X, beta)

                for i in 1:n
                    @inbounds resid[i] = y[i] - xbeta[i]
                end

                for i in 1:n
                    if abs(resid[i]) < threshold
                        @inbounds resid[i] = sign(resid[i]) * threshold
                    end
                    if resid[i] < 0
                        @inbounds resid[i] = abs(q * resid[i])
                    else
                        @inbounds resid[i] = abs((1 - q) * resid[i])
                    end
                end

                broadcast!(/, xstar, X, resid)

                diff = maxdiff(beta0, beta)
            end
        end
        return beta
    end

    function qreg_vcov(y::Vector, X::Matrix, beta::Vector, q::Real)
        resid = y - X * beta
        n = length(resid)
        iqre = quantile(resid, .75) - quantile(resid, .25)
        h = bandwidth_hall_sheather(n, q, .05)
        h1 = min(std(y), iqre / 1.34)
        h2 = quantile(Normal(), q + h) - quantile(Normal(), q - h)
        h = h1 * h2
        # TODO: This line could be optimized, but is a hassle
        fhat0 = 1 / (n * h) * sum(kernel_epanechnikov(resid / h))
        d = copy(resid)
        for i in 1:n
            if d[i] > 0
                @inbounds d[i] = (q / fhat0)^2
            else
                @inbounds d[i] = ((1 - q) / fhat0)^2
            end
        end
        # This pinv() makes me uncomfortable
        xtxi = pinv(X' * X)
        xtdx = broadcast(*, X, d)' * X
        vcov = (xtxi * xtdx) * xtxi
        return vcov
    end

    function qreg(f::Expr, df::AbstractDataFrame, q::Real = 0.5)
        mf = ModelFrame(f, df)
        mm = ModelMatrix(mf)
        mr = vector(model_response(mf))
        coef = qreg_coef(mr, mm.m, q)
        vcov = qreg_vcov(mr, mm.m, coef, q)
        return QRegModel(coef, vcov, mf)
    end

    coef(x::QRegModel) = x.beta

    vcov(x::QRegModel) = x.vcov

    function DataFrames.describe(x::QRegModel)
        # TODO: Insert column names
        se = sqrt(diag(x.vcov))
        out = DataFrame()
        out["Estimate"] = x.beta
        out["Std.Error"] = se
        out["t value"] = x.beta ./ se
        out["2.5%"] = x.beta - 1.96 * se
        out["97.5%"] = x.beta + 1.96 * se
        return out
    end
end


text:
module LJsimulation

  include("Data.jl") #Data structures for the simulation
  export Data

  include("./initial-point.jl") #Initial points to do the simulations
  export initial_point

  include("./wrapcell.jl")
  include("./distpbc.jl") #Computes the distance between two points considering PBC
  include("./linkedlists.jl") #Linked lists of the particles
  export linkedlists

  include("./upair.jl") #Computes the Leonard-Jones potential energy between twoparticles
  include("./u_total_naive.jl") #Computes the total LJ potential energy
  export u_total_naive

  include("./u_total_ll.jl") #Computes the total LJ potential energy using linked lists
  export u_total_ll

  include("./u_total_parallel.jl") #Computes the total LJ potential energy using linked lists and parallelization
  export u_total_parallel

  include("./force_pair.jl")
  include("./force_energy.jl")
  export force_energy!

end


text:
module kiwiComposer
include("wave.jl")

export composer;
end  # module kiwiComposer

max_stars_repo_name: ZaneMuir/kiwi-composer

max_issues_repo_name: ZaneMuir/kiwi-composer

max_forks_repo_name: ZaneMuir/kiwi-composer


text:
#!/usr/bin/julia

function fib(n)
    f = BigInt[1 1; 1 0]
    for i = 1:n
        fn = f ^ i
        println(i, ": ", fn[2, 1])
    end
end

fib(parse(Int, ARGS[1]))

max_stars_repo_name: vfrico/sample-programs

max_issues_repo_name: vfrico/sample-programs

max_forks_repo_name: vfrico/sample-programs


text:
struct Uniform{T<:AbstractSet}
    set::T
end

"""
    Uniform(collection)

Create a uniform categorical distribution over a collection of objects.

The objects in the collection must be unique (this is tested on construction), and will be stored in a `Set`. To avoid this overhead, use `UnsafeUniform`.
"""
function Uniform(c)
    set = Set(c)
    if length(c) > length(set)
        error("""
              Error constructing Uniform($c).

              Objects must be unique (that is, length(Set(c)) == length(c)).
              """
             )
    end
    return Uniform(set)
end

support(d::Uniform) = d.set
sampletype(::Type{Uniform{T}}) where T = eltype(T)
Random.gentype(::Type{Uniform{T}}) where T = eltype(T)

function pdf(d::Uniform, s)
    if s in d.set
        return 1.0/length(d.set)
    else
        return 0.0
    end
end

Base.show(io::IO, m::MIME"text/plain", d::Uniform) = showdistribution(io, m, d, title="Uniform distribution")

"""
    UnsafeUniform(collection)

Create a uniform categorical distribution over a collection of objects.

No checks are performed to ensure uniqueness or check whether an object is actually in the set when evaluating the pdf.
"""
struct UnsafeUniform{T}
    collection::T
end

pdf(d::UnsafeUniform, s) = 1.0/length(d.collection)
support(d::UnsafeUniform) = d.collection
sampletype(::Type{UnsafeUniform{T}}) where T = eltype(T)
Random.gentype(::Type{UnsafeUniform{T}}) where T = eltype(T)

# Common Implementations

const Unif = Union{Uniform,UnsafeUniform}

rand(rng::AbstractRNG, d::Unif) = rand(rng, support(d))
mean(d::Unif) = mean(support(d))
mode(d::Unif) = mode(support(d))

function weighted_iterator(d::Unif)
    p = 1.0/length(support(d))
    return (x=>p for x in support(d))
end

Base.show(io::IO, m::MIME"text/plain", d::UnsafeUniform) = showdistribution(io, m, d, title="UnsafeUniform distribution")


text:
function calc_Vandermonde_xxxx_2D(x,y,ni,nj)

    N=1
    NN=2
    pows=reverse.(Iterators.product(fill(0:N-1,NN)...))[:]
    while length(pows)<nj
        N+=1
        pows=reverse.(Iterators.product(fill(0:N-1,NN)...))[:]
    end

    V=zeros(BigFloat, ni, nj)
    for i=1:ni
        for j=1:nj
            if pows[j][1]-4>=0
                # a1=BigInt(pows[j][1]-1)*BigInt(pows[j][1]-2)*BigInt(pows[j][1]-3)*BigInt(pows[j][1]-4)
                a1=BigInt(pows[j][1]-0)*BigInt(pows[j][1]-1)*BigInt(pows[j][1]-2)*BigInt(pows[j][1]-3)
                V[i,j] = a1*(x[i]^BigInt(pows[j][1]-4))*(y[i]^BigInt(pows[j][2]))
            end
        end
    end

    return V

end


text:
function get_provider(obj::PotentialCoefficientsReader, arg0::jboolean, arg1::jint, arg2::jint)
    return jcall(obj, "getProvider", RawSphericalHarmonicsProvider, (jboolean, jint, jint), arg0, arg1, arg2)
end

function load_data(obj::PotentialCoefficientsReader, arg0::InputStream, arg1::JString)
    return jcall(obj, "loadData", void, (InputStream, JString), arg0, arg1)
end


text:
gather(input::AbstractArray{T}, index::CuArray{Int}) where T = gather(cu(input), index)

function gather(input::CuMatrix{T}, index::CuArray{Int}) where T
    out = CuArrays.zeros(T, size(input,1), size(index)...)
    @inbounds for ind = CartesianIndices(index)
        out[:, ind] = input[:, index[ind]]
    end
    return out
end


text:
# syntax: proto3
using ProtoBuf
import ProtoBuf.meta

mutable struct Dummy <: ProtoType
    value::Float32
    Dummy(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct Dummy

mutable struct Humidity <: ProtoType
    humidity::Float32
    temperature::Float32
    temperature_raw::Float32
    temperature_is_calibrated::Bool
    Humidity(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct Humidity

mutable struct HumidityParams <: ProtoType
    current_temperature::Float32
    HumidityParams(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct HumidityParams

mutable struct UV <: ProtoType
    uv_index::Float32
    oms_risk::AbstractString
    UV(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct UV

mutable struct Pressure <: ProtoType
    pressure::Float32
    altitude::Float32
    temperature::Float32
    Pressure(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct Pressure

mutable struct Imu <: ProtoType
    yaw::Float32
    pitch::Float32
    roll::Float32
    accel_x::Float32
    accel_y::Float32
    accel_z::Float32
    gyro_x::Float32
    gyro_y::Float32
    gyro_z::Float32
    mag_x::Float32
    mag_y::Float32
    mag_z::Float32
    Imu(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct Imu

export Dummy, Humidity, HumidityParams, UV, Pressure, Imu


text:
# Sundials.jl algorithms

# Abstract Types
abstract type SundialsODEAlgorithm{Method, LinearSolver} <: DiffEqBase.AbstractODEAlgorithm end
abstract type SundialsDAEAlgorithm{LinearSolver} <: DiffEqBase.AbstractDAEAlgorithm end

# ODE Algorithms
struct CVODE_BDF{Method, LinearSolver, P, PS} <: SundialsODEAlgorithm{Method, LinearSolver}
    jac_upper::Int
    jac_lower::Int
    krylov_dim::Int
    stability_limit_detect::Bool
    max_hnil_warns::Int
    max_order::Int
    max_error_test_failures::Int
    max_nonlinear_iters::Int
    max_convergence_failures::Int
    prec::P
    psetup::PS
    prec_side::Int
end
Base.@pure function CVODE_BDF(;
    method = :Newton,
    linear_solver = :Dense,
    jac_upper = 0,
    jac_lower = 0,
    non_zero = 0,
    krylov_dim = 0,
    stability_limit_detect = false,
    max_hnil_warns = 10,
    max_order = 5,
    max_error_test_failures = 7,
    max_nonlinear_iters = 3,
    max_convergence_failures = 10,
    prec = nothing,
    psetup = nothing,
    prec_side = 0,
)
    if linear_solver == :Band && (jac_upper == 0 || jac_lower == 0)
        error("Banded solver must set the jac_upper and jac_lower")
    end
    if !(
        linear_solver in (
            :None,
            :Diagonal,
            :Dense,
            :LapackDense,
            :Band,
            :LapackBand,
            :BCG,
            :GMRES,
            :FGMRES,
            :PCG,
            :TFQMR,
            :KLU,
        )
    )
        error("Linear solver choice not accepted.")
    end
    CVODE_BDF{method, linear_solver, typeof(prec), typeof(psetup)}(
        jac_upper,
        jac_lower,
        krylov_dim,
        stability_limit_detect,
        max_hnil_warns,
        max_order,
        max_error_test_failures,
        max_nonlinear_iters,
        max_convergence_failures,
        prec,
        psetup,
        prec_side,
    )
end

struct CVODE_Adams{Method, LinearSolver, P, PS} <:
       SundialsODEAlgorithm{Method, LinearSolver}
    jac_upper::Int
    jac_lower::Int
    krylov_dim::Int
    stability_limit_detect::Bool
    max_hnil_warns::Int
    max_order::Int
    max_error_test_failures::Int
    max_nonlinear_iters::Int
    max_convergence_failures::Int
    prec::P
    psetup::PS
    prec_side::Int
end
Base.@pure function CVODE_Adams(;
    method = :Functional,
    linear_solver = :None,
    jac_upper = 0,
    jac_lower = 0,
    krylov_dim = 0,
    stability_limit_detect = false,
    max_hnil_warns = 10,
    max_order = 12,
    max_error_test_failures = 7,
    max_nonlinear_iters = 3,
    max_convergence_failures = 10,
    prec = nothing,
    psetup = nothing,
    prec_side = 0,
)
    if linear_solver == :Band && (jac_upper == 0 || jac_lower == 0)
        error("Banded solver must set the jac_upper and jac_lower")
    end
    if !(
        linear_solver in (
            :None,
            :Diagonal,
            :Dense,
            :LapackDense,
            :Band,
            :LapackBand,
            :BCG,
            :GMRES,
            :FGMRES,
            :PCG,
            :TFQMR,
            :KLU,
        )
    )
        error("Linear solver choice not accepted.")
    end
    CVODE_Adams{method, linear_solver, typeof(prec), typeof(psetup)}(
        jac_upper,
        jac_lower,
        krylov_dim,
        stability_limit_detect,
        max_hnil_warns,
        max_order,
        max_error_test_failures,
        max_nonlinear_iters,
        max_convergence_failures,
        prec,
        psetup,
        prec_side,
    )
end

struct ARKODE{Method, LinearSolver, MassLinearSolver, T, T1, T2, P, PS} <:
       SundialsODEAlgorithm{Method, LinearSolver}
    stiffness::T
    jac_upper::Int
    jac_lower::Int
    mass_upper::Int
    mass_lower::Int
    krylov_dim::Int
    mass_krylov_dim::Int
    max_hnil_warns::Int
    max_error_test_failures::Int
    max_nonlinear_iters::Int
    max_convergence_failures::Int
    predictor_method::Int
    nonlinear_convergence_coefficient::Float64
    dense_order::Int
    order::Int
    set_optimal_params::Bool
    crdown::Float64
    dgmax::Float64
    rdiv::Float64
    msbp::Int
    itable::T1
    etable::T2
    prec::P
    psetup::PS
    prec_side::Int
end

Base.@pure function ARKODE(
    stiffness = Implicit();
    method = :Newton,
    linear_solver = :Dense,
    mass_linear_solver = :Dense,
    jac_upper = 0,
    jac_lower = 0,
    mass_upper = 0,
    mass_lower = 0,
    non_zero = 0,
    krylov_dim = 0,
    mass_krylov_dim = 0,
    max_hnil_warns = 10,
    max_error_test_failures = 7,
    max_nonlinear_iters = 3,
    max_convergence_failures = 10,
    predictor_method = 0,
    nonlinear_convergence_coefficient = 0.1,
    dense_order = 3,
    order = 4,
    set_optimal_params = false,
    crdown = 0.3,
    dgmax = 0.2,
    rdiv = 2.3,
    msbp = 20,
    adaptivity_method = 0,
    itable = nothing,
    etable = nothing,
    prec = nothing,
    psetup = nothing,
    prec_side = 0,
)
    if linear_solver == :Band && (jac_upper == 0 || jac_lower == 0)
        error("Banded solver must set the jac_upper and jac_lower")
    end
    if !(
        linear_solver in (
            :None,
            :Diagonal,
            :Dense,
            :LapackDense,
            :Band,
            :LapackBand,
            :BCG,
            :GMRES,
            :FGMRES,
            :PCG,
            :TFQMR,
            :KLU,
        )
    )
        error("Linear solver choice not accepted.")
    end
    if !(
        mass_linear_solver in (
            :None,
            :Diagonal,
            :Dense,
            :LapackDense,
            :Band,
            :LapackBand,
            :BCG,
            :GMRES,
            :FGMRES,
            :PCG,
            :TFQMR,
            :KLU,
        )
    )
        error("Mass Matrix Linear solver choice not accepted.")
    end
    ARKODE{
        method,
        linear_solver,
        mass_linear_solver,
        typeof(stiffness),
        typeof(itable),
        typeof(etable),
        typeof(prec),
        typeof(psetup),
    }(
        stiffness,
        jac_upper,
        jac_lower,
        mass_upper,
        mass_lower,
        krylov_dim,
        mass_krylov_dim,
        max_hnil_warns,
        max_error_test_failures,
        max_nonlinear_iters,
        max_convergence_failures,
        predictor_method,
        nonlinear_convergence_coefficient,
        dense_order,
        order,
        set_optimal_params,
        crdown,
        dgmax,
        rdiv,
        msbp,
        itable,
        etable,
        prec,
        psetup,
        prec_side,
    )
end

# DAE Algorithms
struct IDA{LinearSolver, P, PS} <: SundialsDAEAlgorithm{LinearSolver}
    jac_upper::Int
    jac_lower::Int
    krylov_dim::Int
    max_order::Int
    max_error_test_failures::Int
    nonlinear_convergence_coefficient::Float64
    max_nonlinear_iters::Int
    max_convergence_failures::Int
    nonlinear_convergence_coefficient_ic::Float64
    max_num_steps_ic::Int
    max_num_jacs_ic::Int
    max_num_iters_ic::Int
    max_num_backs_ic::Int
    use_linesearch_ic::Bool
    init_all::Bool
    prec::P
    psetup::PS
    prec_side::Int
end
Base.@pure function IDA(;
    linear_solver = :Dense,
    jac_upper = 0,
    jac_lower = 0,
    krylov_dim = 0,
    max_order = 5,
    max_error_test_failures = 7,
    max_nonlinear_iters = 3,
    nonlinear_convergence_coefficient = 0.33,
    nonlinear_convergence_coefficient_ic = 0.0033,
    max_num_steps_ic = 5,
    max_num_jacs_ic = 4,
    max_num_iters_ic = 10,
    max_num_backs_ic = 100,
    use_linesearch_ic = true,
    init_all = false,
    max_convergence_failures = 10,
    prec = nothing,
    psetup = nothing,
    prec_side = 0,
)
    if linear_solver == :Band && (jac_upper == 0 || jac_lower == 0)
        error("Banded solver must set the jac_upper and jac_lower")
    end
    if !(
        linear_solver in (
            :None,
            :Diagonal,
            :Dense,
            :LapackDense,
            :Band,
            :LapackBand,
            :BCG,
            :GMRES,
            :FGMRES,
            :PCG,
            :TFQMR,
            :KLU,
        )
    )
        error("Linear solver choice not accepted.")
    end
    IDA{linear_solver, typeof(prec), typeof(psetup)}(
        jac_upper,
        jac_lower,
        krylov_dim,
        max_order,
        max_error_test_failures,
        nonlinear_convergence_coefficient,
        max_nonlinear_iters,
        max_convergence_failures,
        nonlinear_convergence_coefficient_ic,
        max_num_steps_ic,
        max_num_jacs_ic,
        max_num_iters_ic,
        max_num_backs_ic,
        use_linesearch_ic,
        init_all,
        prec,
        psetup,
        prec_side,
    )
end

method_choice(alg::SundialsODEAlgorithm{Method}) where {Method} = Method
method_choice(alg::SundialsDAEAlgorithm) = :Newton
linear_solver(
    alg::SundialsODEAlgorithm{Method, LinearSolver},
) where {Method, LinearSolver} = LinearSolver
linear_solver(alg::SundialsDAEAlgorithm{LinearSolver}) where {LinearSolver} = LinearSolver


text:
function optflow(first_img::AbstractArray{T, 2}, second_img::AbstractArray{T,2}, points::Array{SVector{2, Float64}, 1},  algorithm::LucasKanade) where T <: Gray
    displacement = Array{SVector{2, Float64}, 1}(undef, length(points))
    for i in eachindex(displacement)
            displacement[i] = SVector{2, Float64}(0.0, 0.0)
    end
    optflow!(first_img, second_img, points, displacement,  algorithm)
end

function optflow!(first_img::AbstractArray{T, 2}, second_img::AbstractArray{T,2}, points::Array{SVector{2, Float64}, 1}, displacement::Array{SVector{2, Float64}, 1}, algorithm::LucasKanade) where T <: Gray

    # Convert (row, col) convention to (x,y) coordinate convention.
    points = map(x -> SVector(last(x), first(x)), points)
    map!(x -> SVector(last(x), first(x)), displacement, displacement)

    # Replace NaN with zero in both images.
    first_img = map(x -> isnan(x) ? zero(x) : x, first_img)
    second_img = map(x -> isnan(x) ? zero(x) : x, second_img)

    # Construct Gaussian pyramids for both the images.
    first_pyramid, second_pyramid =  construct_pyramids(first_img, second_img,  algorithm.pyramid_levels)
    flow, status, reliability = initialise_arrays(length(displacement))

    for i = algorithm.pyramid_levels+1:-1:1
        # Extrapolate the second image so as to get intensity at non-integral points.
        itp = interpolate(second_pyramid[i], BSpline(Linear()))
        etp = extrapolate(itp, zero(eltype(second_pyramid[i])))

        Ix, Iy = imgradients(first_pyramid[i], KernelFactors.scharr, Fill(zero(eltype(first_pyramid[i])))) # TODO: Fix row,col vs x,y convention
        # Calculate Ixx, Iyy and Ixy taking Gaussian weights for pixels in the grid.
        # Return integral tables to facilitate fast sums over regions in Ixx, Iyy and Ixy.
        Ixx_integral_table, Ixy_integral_table, Iyy_integral_table =  compute_partial_derivatives(Ix, Iy)

        inner_bounds = map(i -> first(i) +  algorithm.window_size:last(i) -  algorithm.window_size, axes(first_pyramid[i]))
        # Running the algorithm for each input point.
        for j = 1:length(displacement)
            if status[j]
                # Position of the point in the current pyramid level.
                point = determine_pyramid_coordinates(points, i, j)
                # If the window falls inside the first image we pre-compute the
                # spatial gradient matrix because in many instance the window
                # will also be completely contained inside the second image.  If
                # it turns out that the  corresponding window is not totally
                # contained inside the second window then we will recompute the
                # spatial gradient matrix again later.
                if lies_in(inner_bounds, point)
                    square_grid = SVector{2}(point[1] -  algorithm.window_size:point[1] +  algorithm.window_size, point[2] -  algorithm.window_size:point[2] +  algorithm.window_size)
                    G = compute_spatial_gradient(square_grid, Ixx_integral_table, Ixy_integral_table, Iyy_integral_table)
                end

                pyramid_contribution = SVector{2}(0.0,0.0)
                for k = 1:algorithm.iterations
                    putative_flow = displacement[j] + pyramid_contribution
                    putative_correspondence = point + putative_flow
                    if lies_in(axes(second_img), putative_correspondence)
                        grid, offsets, is_truncated_window = get_grid(first_pyramid[i], point, putative_flow, algorithm.window_size)
                        # If the grid is not a square then we cannot reuse our precomputed spatial gradient.
                        if is_truncated_window
                            G = compute_spatial_gradient(grid, Ixx_integral_table, Ixy_integral_table, Iyy_integral_table)
                        end
                        A = @inbounds view(first_pyramid[i], grid[1], grid[2])
                        Ix_window = @inbounds view(Ix, grid[1], grid[2])
                        Iy_window = @inbounds view(Iy, grid[1], grid[2])
                        estimated_flow = compute_flow_vector(point + putative_flow, G, A, Ix_window, Iy_window, grid, offsets, etp)

                        pyramid_contribution += estimated_flow
                        # Minimum eigenvalue in the matrix is used as a measure of reliability.
                        F =  eigen(G)
                        min_eigenvalue = F.values[1] < F.values[2] ? F.values[1] : F.values[2]
                        grid_size = prod(length.(grid))
                        reliability[j] = min_eigenvalue / grid_size

                        # If the point is too unreliable then it should be declared lost.
                        if is_lost(first_pyramid[i], point + pyramid_contribution, min_eigenvalue, algorithm.eigenvalue_threshold, grid_size)
                            declare_lost!(status, flow, j)
                            break
                        end
                    else
                        # In this instance the corresponding point lies outside of the second image so we declare it as lost.
                        declare_lost!(status, flow, j)
                        break
                    end
                end

                # If the optical flow is too large then the point should be declared lost.
                if is_lost(SVector(displacement[j]), 2* algorithm.window_size)
                    declare_lost!(status, flow, j)
                end

                if status[j]
                    flow[j] = pyramid_contribution
                    # A guess for the next level of the pyramid.
                    displacement[j] = 2*(displacement[j] + flow[j])
                end
            end
        end
    end

    # Final optical flow.
    output_flow = 0.5*displacement
    return output_flow, status
end

function declare_lost!(status, flow, j)
    status[j] = false
    flow[j] = SVector{2}(0.0,0.0)
end

function construct_pyramids(first_img, second_img, pyramid_levels)
    first_pyramid = gaussian_pyramid(first_img, pyramid_levels, 2, 1.0)
    second_pyramid = gaussian_pyramid(second_img, pyramid_levels, 2, 1.0)
    return first_pyramid, second_pyramid
end

function initialise_arrays(N::Int)
    flow = Array{SVector{2, Float64}, 1}(undef, N)
    for i in eachindex(flow)
        flow[i] = SVector{2, Float64}(0.0, 0.0)
    end
    status = trues(N)
    reliability = zeros(Float64, N)
    return flow, status, reliability
end

function compute_partial_derivatives(Ix, Iy)
    Ixx = imfilter(Ix.*Ix, Kernel.gaussian(4))
    Iyy = imfilter(Iy.*Iy, Kernel.gaussian(4))
    Ixy = imfilter(Ix.*Iy, Kernel.gaussian(4))

    Ixx_integral_table = integral_image(Ixx)
    Ixy_integral_table = integral_image(Ixy)
    Iyy_integral_table = integral_image(Iyy)
    return Ixx_integral_table, Ixy_integral_table, Iyy_integral_table
end

function compute_flow_vector(corresponding_point, G, A, Ix, Iy, grid, offsets, etp)
    # Solve a linear system of equations in order to determine the flow vector.
    bk = prepare_linear_system(corresponding_point, A, Ix, Iy, grid, offsets, etp)
    G_inv = pinv2x2(G) #TODO: Add the pinv2x2 method to the StaticArrays package.
    ηk = G_inv*bk
end

function determine_pyramid_coordinates(points, i, j)
    px = floor(Int,(points[j][1]) / 2^i)
    py = floor(Int,(points[j][2]) / 2^i)
    point = SVector{2}(px ,py)
end

function compute_spatial_gradient(grid, Ixx_integral::AbstractArray, Ixy_integral::AbstractArray, Iyy_integral::AbstractArray)
    sum_Ixx = boxdiff(Ixx_integral, grid[1], grid[2])
    sum_Ixy = boxdiff(Ixy_integral, grid[1], grid[2])
    sum_Iyy = boxdiff(Iyy_integral, grid[1], grid[2])
    # Spatial gradient matrix.
    G = SMatrix{2,2,Float64}(sum_Ixx ,sum_Ixy, sum_Ixy, sum_Iyy)
end

function prepare_linear_system(corresponding_point, A, Ix, Iy, grid, offsets, etp)
    bx = 0.0
    by = 0.0
    P, Q = size(A)
    # Evaluates sum(δI.*I_x) and sum(δI.*I_y), where δI = A .- B
    @inbounds begin
        for q = 1:Q
           for p = 1:P
                Apq = A[p,q]
                r = corresponding_point[1] + offsets[1][p]
                c = corresponding_point[2] + offsets[2][q]
                Bpq = etp(r, c) # Subpixel precision
                bx = bx + (Apq - Bpq)*Ix[p,q]
                by = by + (Apq - Bpq)*Iy[p,q]
            end
        end
    end
    b = SVector{2}(bx, by)
end

function in_image(img::AbstractArray{T, 2}, point::SVector{2, U}, window::Int) where {T <: Gray, U <: Union{Int, Float64}}
    first_axis = first(axes(img))
    second_axis = last(axes(img))
    if point[1] < first(first_axis) || point[2] < first(second_axis) || point[1] > last(first_axis) || point[2] > last(second_axis)
        return false
    else
        return true
    end
end

function is_lost(img::AbstractArray{T, 2}, point::SVector{2, U}, min_eigenvalue::Float64, eigenvalue_threshold::Float64, window::Int) where {T <: Gray, U <: Union{Int, Float64}}
    if !(in_image(img, point, window))
        return true
    else
        val = min_eigenvalue/window

        if val < eigenvalue_threshold
            return true
        else
            return false
        end
    end
end

function is_lost(point::SVector{2, Float64}, window_size::Int)
    if point[1] > window_size || point[2] > window_size
        return true
    else
        return false
    end
end

function lies_in(area::Tuple{AbstractRange{Int},AbstractRange{Int}}, point::SVector{2, T}) where T <: Union{Int,Float64}
    if first(area[1]) <= point[1] && last(area[1]) >= point[1] && first(area[2]) <= point[2] && last(area[2]) >= point[2]
        return true
    else
        return false
    end
end

function get_grid(img::AbstractArray{T, 2}, point::SVector{2, U}, displacement::SVector{2, Float64}, window_size::Int) where {T <: Gray, U <: Union{Int, Float64}}
    first_axis, second_axis = axes(img)
    allowed_area = map(i -> first(i) + window_size:last(i) - window_size, axes(img))
    new_point = point + displacement

    if !lies_in(allowed_area, new_point) || !lies_in(allowed_area, point)
        is_truncated_window = true
    else
        is_truncated_window = false
    end

    # TODO: Convert these loops to explicit formulae.
    w_up = 0
    while point[1] - (w_up + 1) > first(first_axis)  && new_point[1] - (w_up + 1) >  first(first_axis)  &&  (w_up + 1) < window_size
        w_up += 1
    end

    w_down = 0
    while point[1] + (w_down + 1)  < last(first_axis)  && new_point[1] + (w_down + 1) <  last(first_axis)  && (w_down + 1) < window_size
        w_down += 1
    end

    w_left = 0
    while point[2] - (w_left + 1) > first(second_axis) && new_point[2] - (w_left + 1) >  first(second_axis) && (w_left + 1) < window_size
        w_left += 1
    end

    w_right = 0
    while point[2] + (w_right + 1) < last(second_axis) && new_point[2] + (w_right + 1) >  last(second_axis) && (w_right + 1) < window_size
        w_right += 1
    end

    grid_1 = (point[1] - w_up:point[1] + w_down, point[2] - w_left:point[2] + w_right)
    offsets = (UnitRange(-w_up, w_down), UnitRange(-w_left, w_right))

    return grid_1, offsets, is_truncated_window
end


text:
using Bridge
using StaticArrays, LinearAlgebra
import Bridge: b, σ, B, β, a, constdiff
const ℝ = SVector{N,T} where {N,T}

struct Lorenz{T} <: ContinuousTimeProcess{ℝ{3,T}}
    θ₁::T
    θ₂::T
    θ₃::T
    σ₁::T
    σ₂::T
    σ₃::T
    function Lorenz(θ₁::T, θ₂::T, θ₃::T, σ₁::T, σ₂::T, σ₃::T) where T
        new{T}(θ₁, θ₂, θ₃, σ₁, σ₂, σ₃)
    end
end

function b(t, x, P::Lorenz)
    ℝ{3}(P.θ₁*(x[2]-x[1]),
         P.θ₂*x[1] - x[2] - x[1]*x[3],
         x[1]*x[2] - P.θ₃*x[3])
end

function σ(t, x, P::Lorenz)
    @SMatrix[ P.σ₁ 0.0 0.0;
              0.0 P.σ₂ 0.0;
              0.0 0.0 P.σ₃]
end

constdiff(::Lorenz) = true
clone(P::Lorenz, θ) = Lorenz(θ...)
params(P::Lorenz) = [P.θ₁, P.θ₂, P.θ₃, P.σ₁, P.σ₂, P.σ₃]

struct LorenzAux{O,R,S1,S2} <: ContinuousTimeProcess{ℝ{3,R}}
    θ₁::R
    θ₂::R
    θ₃::R
    σ₁::R
    σ₂::R
    σ₃::R
    t::Float64
    u::S1
    T::Float64
    v::S2
    aux::Float64

    function LorenzAux(θ₁::R, θ₂::R, θ₃::R, σ₁::R, σ₂::R, σ₃::R, t, u::S1, T,
                       v::S2, ::O, aux=0.0) where {O,R,S1,S2}
        new{O,R,S1,S2}(θ₁, θ₂, θ₃, σ₁, σ₂, σ₃, t, u, T, v, aux)
    end
end

observables(::LorenzAux{O}) where O = O()

# Auxiliary diffusion when coordinates [1,2,3] are observed
# ---------------------------------------------------------
function B(t, P::LorenzAux{Val{(true,true,true)}})
    @SMatrix [ -P.θ₁  P.θ₁  0.0;
               P.θ₂-P.v[3]  -1.0  -P.v[1];
               P.v[2]  P.v[1]  -P.θ₃]
end

function β(t, P::LorenzAux{Val{(true,true,true)}})
    @SVector [0.0,
              P.v[1]*P.v[3],
              -P.v[1]*P.v[2]]
end

function σ(t, P::LorenzAux{Val{(true,true,true)}})
    @SMatrix[ P.σ₁ 0.0 0.0;
              0.0 P.σ₂ 0.0;
              0.0 0.0 P.σ₃]
end

a(t, P::LorenzAux) = σ(t,P) * σ(t, P)'
σ(t, x, P::LorenzAux) = σ(t, P)
b(t, x, P::LorenzAux) = B(t, P)*x + β(t, P)

constdiff(::LorenzAux) = true
clone(P::LorenzAux, θ) = LorenzAux(θ..., P.t, P.u, P.T, P.v, observables(P), P.aux)
clone(P::LorenzAux, θ, v) = LorenzAux(θ..., P.t, zero(v), P.T, v, observables(P), P.aux)
params(P::LorenzAux) = [P.θ₁, P.θ₂, P.θ₃, P.σ₁, P.σ₂, P.σ₃]
depends_on_params(::LorenzAux) = (1,2,3,4,5,6)

# Auxiliary diffusion when coordinates [1,2] are observed
# ---------------------------------------------------------
observables(::LorenzAux{O}) where O = O()

function B(t, P::LorenzAux{Val{(true,true,false)}})
    @SMatrix [ -P.θ₁  P.θ₁  0.0;
               P.θ₂-P.aux  -1.0  -P.v[1];
               P.v[2]  P.v[1]  -P.θ₃]
end

function β(t, P::LorenzAux{Val{(true,true,false)}})
    @SVector [0.0,
              P.v[1]*P.aux,
              -P.v[1]*P.v[2]]
end

function σ(t, P::LorenzAux{Val{(true,true,false)}})
    @SMatrix[ P.σ₁ 0.0 0.0;
              0.0 P.σ₂ 0.0;
              0.0 0.0 P.σ₃]
end


text:
struct BIOES
    tag2id::Dict{String,Int}
    id2tag::Vector{String}
end
BIOES() = BIOES(Dict{String,Int}(), String[])

function encode(tagset::BIOES, tags::Vector{String})
    map(tags) do tag
        get!(tagset.tag2id, tag) do
            push!(tagset.id2tag, tag)
            length(tagset.tag2id) + 1
        end
    end
end

function decode(tagset::BIOES, ids::Vector{Int})
    id2tag = tagset.id2tag
    spans = Tuple{Int,Int,String}[]
    bpos = 0
    for i = 1:length(ids)
        tag = id2tag[ids[i]]
        tag == "O" && continue
        startswith(tag,"B-") && (bpos = i)
        startswith(tag,"S-") && (bpos = i)
        nexttag = i == length(ids) ? "O" : id2tag[ids[i+1]]
        if (startswith(tag,"S-") || startswith(tag,"E-")) && bpos > 0
            basetag = id2tag[ids[bpos]][3:end]
            push!(spans, (bpos,i,basetag))
            bpos = 0
        end
    end
    spans
end

max_stars_repo_name: hshindo/Merlin-Examples

max_issues_repo_name: hshindo/Merlin-Examples

max_forks_repo_name: hshindo/Merlin-Examples