text:
#=
encoding.jl; This file implements functions to encode token-based data, be it alphabets, strings or full classes of text.
=#

# SEQUENCE EMBEDDING

# precomputing of the n-grams

ngrams_recursion(d::AbstractDict{K,V}, hdv::AbstractHDV) where {K,V<:AbstractHDV} =
    Dict(k => hdv * Π(v, 1) for (k, v) in d)

ngrams_recursion(d::AbstractDict{K,V}, hdv::AbstractHDV) where {K,V<:AbstractDict} =
    Dict(k => ngrams_recursion(v, hdv) for (k, v) in d)

function ngrams_recursion(d::AbstractDict{K,V}, hdvs) where {K,V}
    result = Dict{K,Dict{K,V}}()
    for k in keys(d)
        result[k] = ngrams_recursion(d, hdvs[k])
    end
    return result
end

# structure to store the N-grams
struct NGrams{V,D}
    d::D
    NGrams(V::Type, d::AbstractDict) = new{V,typeof(d)}(d)
end

Base.show(io::IO, ngrams::NGrams) = print("n-gram embedding of order $(order(ngrams)) for type $(vectortype(ngrams))")

# TODO: consider using an N-dim array?

# Generates nested dictionary with ngrams
# needs to be hard-coded for dispatch
#=
compute_1_grams(hdvs, alphabet=keys(hdvs)) = NGrams(eltype(hdvs), Dict(c=>hdvs[c] for c in alphabet))
compute_2_grams(hdvs, alphabet=1:length(hdvs)) = NGrams(eltype(hdvs), ngrams_recursion(compute_1_grams(hdvs, alphabet).d, hdvs))
compute_3_grams(hdvs, alphabet=1:length(hdvs)) = NGrams(eltype(hdvs), ngrams_recursion(compute_2_grams(hdvs, alphabet).d, hdvs))
compute_4_grams(hdvs, alphabet=1:length(hdvs)) = NGrams(eltype(hdvs), ngrams_recursion(compute_3_grams(hdvs, alphabet).d, hdvs))
compute_5_grams(hdvs, alphabet=1:length(hdvs)) = NGrams(eltype(hdvs), ngrams_recursion(compute_4_grams(hdvs, alphabet).d, hdvs))
compute_6_grams(hdvs, alphabet=1:length(hdvs)) = NGrams(eltype(hdvs), ngrams_recursion(compute_5_grams(hdvs, alphabet).d, hdvs))
compute_7_grams(hdvs, alphabet=1:length(hdvs)) = NGrams(eltype(hdvs), ngrams_recursion(compute_6_grams(hdvs, alphabet).d, hdvs))
compute_8_grams(hdvs, alphabet=1:length(hdvs)) = NGrams(eltype(hdvs), ngrams_recursion(compute_7_grams(hdvs, alphabet).d, hdvs))
=#
compute_1_grams(hdvs, alphabet=keys(hdvs)) = NGrams(eltype(hdvs), Dict(c=>hdvs[c] for c in alphabet))
compute_2_grams(hdvs, alphabet=keys(hdvs)) = NGrams(eltype(hdvs), ngrams_recursion(compute_1_grams(hdvs, alphabet).d, hdvs))
compute_3_grams(hdvs, alphabet=keys(hdvs)) = NGrams(eltype(hdvs), ngrams_recursion(compute_2_grams(hdvs, alphabet).d, hdvs))
compute_4_grams(hdvs, alphabet=keys(hdvs)) = NGrams(eltype(hdvs), ngrams_recursion(compute_3_grams(hdvs, alphabet).d, hdvs))
compute_5_grams(hdvs, alphabet=keys(hdvs)) = NGrams(eltype(hdvs), ngrams_recursion(compute_4_grams(hdvs, alphabet).d, hdvs))
compute_6_grams(hdvs, alphabet=keys(hdvs)) = NGrams(eltype(hdvs), ngrams_recursion(compute_5_grams(hdvs, alphabet).d, hdvs))
compute_7_grams(hdvs, alphabet=keys(hdvs)) = NGrams(eltype(hdvs), ngrams_recursion(compute_6_grams(hdvs, alphabet).d, hdvs))
compute_8_grams(hdvs, alphabet=keys(hdvs)) = NGrams(eltype(hdvs), ngrams_recursion(compute_7_grams(hdvs, alphabet).d, hdvs))

# get the embedding on the k-gram that starts at position i
get_gram_embedding(sequence, i, embeddings::Dict{K,V}) where {K,V<:Dict} =
            get_gram_embedding(sequence, i+1, embeddings[sequence[i]])
get_gram_embedding(sequence, i, embeddings::Dict{K,V}) where {K,V<:AbstractHDV} = embeddings[sequence[i]]
get_gram_embedding(sequence, i, embeddings::NGrams) = get_gram_embedding(sequence, i, embeddings.d)

order(embeddings::Dict{K,V}) where {K,V<:Dict} = order(first(embeddings::Dict{K,V})[2]) + 1
order(embeddings::Dict) = 1
order(embeddings::NGrams) = order(embeddings.d)

vectortype(embeddings::NGrams{V,D}) where {V, D} = V

get_leaf(d::AbstractDict{K,V}) where {K,V<:AbstractDict} = get_leaf(d[first(keys(d))])
get_leaf(d::AbstractDict) = d[first(keys(d))]

similar_vector(ngrams::NGrams) = similar(get_leaf(ngrams.d))

function sequence_embedding!(result::AbstractHDV, sequence, token_vectors, w=3)
    fill!(result.v, zero(eltype(result)))
    tmp = similar(result)
    n = length(sequence)
    for i in 1:n-w
        fill!(tmp.v, neutralbind(tmp))
        for k in (w-1):-1:0
            v = token_vectors[sequence[i+k]]
            offsetcombine!(tmp.v, bindfun(tmp), tmp.v, v.v, v.offset + k)
        end
        offsetcombine!(result.v, aggfun(result), result.v, tmp.v, 0)
    end
    result.m = length(sequence)-w
    normalize!(result)
    return result
end

function sequence_embedding!(result::AbstractHDV, sequence, ngrams_embedding::NGrams)
    fill!(result.v, zero(eltype(result)))
    n = length(sequence)
    w = order(ngrams_embedding)
    for i in 1:n-w
        hdv = get_gram_embedding(sequence, i, ngrams_embedding)
        offsetcombine!(result.v, aggfun(result), result.v, hdv.v, 0)
    end
    result.m = length(sequence)-w
    normalize!(result)
    return result
end

sequence_embedding(sequence, token_vectors::AbstractVector{V}, args...) where {V<:AbstractHDV} =
                sequence_embedding!(similar(first(token_vectors)), sequence, token_vectors, args...)

sequence_embedding(sequence, token_vectors::Dict{T,V}, args...) where {T,V<:AbstractHDV} =
                sequence_embedding!(similar(first(token_vectors)[2]), sequence, token_vectors, args...)

sequence_embedding(sequence, ngrams_embedding::NGrams, args...) =
                sequence_embedding!(similar_vector(ngrams_embedding), sequence, ngrams_embedding, args...)


text:
module JuliaSteps

export AbstractStep, GaussStep, PoissonStep, AppendSteps!, plot
using DataStructures

# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#    AbstractStep
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# Step must have functions:
# - `ErrorSum(this, ξ₀, ξ₁)` that calculate error in data[ξ₀:ξ₁].
# - `AppendSteps!(this)` that conduct step detection.
#

abstract type AbstractStep end

function InitStep!(this::AbstractStep, data, p)
    this.data = Vector{Float64}(data)
    this.len = length(data)
    if 0 < p < 1
        this.penalty = log(p/(1-p))
    else
        this.penalty = -0.5log(this.len)
    end
    this.steplist = Vector{Int}([1, this.len + 1])
    this.fit = zeros(Float64, this.len)
    return this
end

function Finalize(this::AbstractStep)
    this.μlist = zeros(Float64, length(this.steplist)-1)

    for i in 1:length(this.steplist)-1
        ξ₀ = this.steplist[i]
        ξ₁ = this.steplist[i+1]
        μ = sum(@views(this.data[ξ₀:ξ₁-1]))/(ξ₁-ξ₀)
        this.fit[ξ₀:ξ₁-1] .= μ
        this.μlist[i] = μ
    end
    this.Δμlist = abs.(diff(this.μlist))
    return
end

# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#    GaussStep
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

struct GaussMoment
    fw1::Vector{Float64}
    fw2::Vector{Float64}
    bw1::Vector{Float64}
    bw2::Vector{Float64}
    total1::Float64
    total2::Float64
    # empty struct
    GaussMoment() = new([], [], [], [], -1., -1.)
    # only forward
    GaussMoment(fw1::Vector{Float64}, fw2::Vector{Float64}, nothing, t1::Float64, t2::Float64) = new(fw1, fw2, Vector{Float64}(), Vector{Float64}(), t1, t2)
    # only backward
    GaussMoment(nothing, bw1::Vector{Float64}, bw2::Vector{Float64}, t1::Float64, t2::Float64) = new(Vector{Float64}(), Vector{Float64}(), bw1, bw2, t1, t2)
    function GaussMoment(data::Vector{Float64})
        fw1 = cumsum(data[1:end-1])
        fw2 = cumsum(data[1:end-1].^2)
        t1 = fw1[end] + data[end]
        t2 = fw2[end] + data[end]^2
        bw1 = t1 .- fw1
        bw2 = t2 .- fw2
        new(fw1, fw2, bw1, bw2, t1, t2)
    end
end

function GetOptimalSplitter(m::GaussMoment)
    N = length(m.fw1)
    χ²all = Getχ²(m)
    χ² = m.fw2 - m.fw1.^2 ./ vec(1:N) + m.bw2 - m.bw1.^2 ./ vec(N:-1:1)
    ξ = argmin(χ²)
    return χ²[ξ] - χ²all, ξ+1
end

Getχ²(m::GaussMoment) = m.total2 - m.total1^2/(length(m.fw1)+1)

function complement!(m::GaussMoment)
    if isempty(m.fw1)
        append!(m.fw1, m.total1 .- m.bw1)
        append!(m.fw2, m.total2 .- m.bw2)
    elseif isempty(m.bw1)
        append!(m.bw1, m.total1 .- m.fw1)
        append!(m.bw2, m.total2 .- m.fw2)
    end
end

function split(m::GaussMoment, i::Int)
    g₁ = GaussMoment(m.fw1[1:i-2], m.fw2[1:i-2], nothing, m.fw1[i-1], m.fw2[i-1])
    g₂ = GaussMoment(nothing, m.bw1[i:end], m.bw2[i:end], m.bw1[i-1], m.bw2[i-1])
    complement!(g₁)
    complement!(g₂)
    return g₁, g₂
end

mutable struct GaussStep <: AbstractStep
    data::Vector{Float64}
    len::Int
    penalty::Float64
    steplist::Vector{Int}
    fit::Vector{Float64}
    μlist::Vector{Float64}
    Δμlist::Vector{Float64}

    function GaussStep(data, p=-1.)
        this = new()
        InitStep!(this, data, p)
        return this
    end
end

function AppendSteps!(this::GaussStep)
    g = GaussMoment(this.data)
    χ² = Getχ²(g)
    heap = BinaryMinHeap{Tuple{Float64, Int, Int, GaussMoment}}()
    push!(heap, (GetOptimalSplitter(g)..., 1, g))

    while true
        Δχ², Δξ, ξ₀, g = pop!(heap)
        ΔlogL = this.penalty - this.len*log(1 + Δχ²/χ²) / 2
        if ΔlogL > 0
            ξ = ξ₀ + Δξ - 1
            g₁, g₂ = split(g, Δξ)
            length(g₁.fw1) > 1 && push!(heap, (GetOptimalSplitter(g₁)..., ξ₀, g₁))
            length(g₂.fw1) > 1 && push!(heap, (GetOptimalSplitter(g₂)..., ξ, g₂))
            push!(this.steplist, ξ)
            χ² += Δχ²
        else
            break
        end
    end
    sort!(this.steplist)
    Finalize(this)
    return
end

# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#    PoissonStep
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

struct PoissonMoment
    fw1::Vector{Float64}
    bw1::Vector{Float64}
    total1::Float64
    # empty struct
    PoissonMoment() = new([], [], -1.)
    # only forward
    PoissonMoment(fw1::Vector{Float64}, nothing, t1::Float64) = new(fw1, Vector{Float64}(), t1)
    # only backward
    PoissonMoment(nothing, bw1::Vector{Float64}, t1::Float64) = new(Vector{Float64}(), bw1, t1)
    function PoissonMoment(data::Vector{Float64})
        fw1 = cumsum(data[1:end-1])
        t1 = fw1[end] + data[end]
        bw1 = t1 .- fw1
        new(fw1, bw1, t1)
    end
end

GetSlogμ(m::PoissonMoment) = m.total1 * log(m.total1 / (length(m.fw1)+1))

function GetOptimalSplitter(m::PoissonMoment)
    N = length(m.fw1)
    Slogμall = GetSlogμ(m)
    Slogμ = m.fw1 .* log.((m.fw1 .+ 1e-12) ./ vec(1:N)) + m.bw1 .* log.((m.bw1 .+ 1e-12) ./ vec(N:-1:1))
    ξ = argmax(Slogμ)
    return Slogμ[ξ] - Slogμall, ξ+1
end

function complement!(m::PoissonMoment)
    if isempty(m.fw1)
        append!(m.fw1, m.total1 .- m.bw1)
    elseif isempty(m.bw1)
        append!(m.bw1, m.total1 .- m.fw1)
    end
end

function split(m::PoissonMoment, i::Int)
    p₁ = PoissonMoment(m.fw1[1:i-2], nothing, m.fw1[i-1])
    p₂ = PoissonMoment(nothing, m.bw1[i:end], m.bw1[i-1])
    complement!(p₁)
    complement!(p₂)
    return p₁, p₂
end

mutable struct PoissonStep <: AbstractStep
    data::Vector{Float64}
    len::Int
    penalty::Float64
    steplist::Vector{Int}
    fit::Vector{Float64}
    μlist::Vector{Float64}
    Δμlist::Vector{Float64}

    function PoissonStep(data, p=-1)
        this = new()
        InitStep!(this, data, p)
        return this
    end

end

function AppendOneStep!(this::PoissonStep, m::PoissonMoment, start=1)
    if length(m.fw1) < 2
        return
    end
    Δslogμ, Δξ = GetOptimalSplitter(m)
    ΔlogL = this.penalty + Δslogμ
    if ΔlogL > 0
        ξ = start + Δξ - 1
        m1, m2 = split(m, Δξ)
        append!(this.steplist, ξ)
        AppendOneStep!(this, m1, start)
        AppendOneStep!(this, m2, ξ)
    end
    return
end

function AppendSteps!(this::PoissonStep)
    m = PoissonMoment(this.data)
    AppendOneStep!(this, m)
    sort!(this.steplist)
    Finalize(this)
    return
end

end

ext: jl

max_stars_repo_name: hanjinliu/JuliaSteps

max_issues_repo_name: hanjinliu/JuliaSteps

max_forks_repo_name: hanjinliu/JuliaSteps


text:
using ScikitLearn
using RDatasets
using CSV
using DataFrames
using MLLabelUtils
using StatsBase: standardize, mean
@sk_import model_selection:train_test_split
@sk_import preprocessing:OneHotEncoder
iris = dataset("datasets", "iris")

# If you want to use ScikitLearn, you need to convert DataFrames to Arrays
# ScikitLearn only accepts data as arrays.
X = convert(Array, iris[!, 1:4])
y = convert(Array, iris[!, 5])

# train_test_split is a function from ScikitLearn, and splits the data
# according to the value passed to test_size. Here we are saying
# that we need to allocate 33% of the raw data as test data
X_train, X_test, y_train, y_test = train_test_split(X, y; test_size=0.33, random_state=42)

data = dataset("datasets", "airquality")
describe(data)
describe(dropmissing(data))

dropmissing!(data) # <- this will permananty change the data by dropping missing values

data = dataset("datasets", "airquality") # re-loading because we need a fresh copy of the data for next e.g.
# In some other cases, you just want to replace the missing values with the mean of the data or with 0.
# In those cases you can use coalesce function from julia's base library

# To replace the missing values with 0
data.Ozone = coalesce.(data.Ozone, 0)
# Let's check the data once again
describe(data)

# For Solar.R let's replace it with the mean
data."Solar.R" = coalesce.(data."Solar.R", mean(skipmissing(data."Solar.R")))
describe(data)

standardize(UnitRangeTransform, X; dims=1)
standardize(ZScoreTransform, X; dims=1, center=false, scale=true)

onehotencoded = convertlabel(LabelEnc.OneOfK{Float32}, y; obsdim=1)

using Discretizers
bin_edges = [0.1, 0.5, 0.7]
lindisc = LinearDiscretizer(bin_edges);

cat = [:yes, :no, :maybe]
catdisc = CategoricalDiscretizer(cat)
decode(catdisc, z)
z = encode(lindisc, a)

max_stars_repo_name: coinslab/ComputationalCognitiveModeling

max_issues_repo_name: coinslab/ComputationalCognitiveModeling

max_forks_repo_name: coinslab/ComputationalCognitiveModeling


text:
include("help.jl")
function stairs!(r) #Главная функция
    movements!(r,Sud)
    movements!(r,West)
    #Мы в левом нижнем углу

    side = Ost
    num_hor = moves!(r,side)
    movements!(r,West)
    mark_up!(r,num_hor)
end

function mark_up!(r::Robot,size::Int)
    i = 0
    while (1 <= size)
        for _ in 1:size
            move!(r,Ost)
            putmarker!(r)
        end
        if (isborder(r,West) == false)
            movements!(r,West)
            putmarker!(r)
        end
        if (isborder(r,Nord) == false)
        move!(r,Nord)
        size = size - 1
        else
            break
        end
    end
end

max_stars_repo_name: MCPavlo/mirea-progs

max_issues_repo_name: MCPavlo/mirea-progs

max_forks_repo_name: MCPavlo/mirea-progs


text:
function nlp_cvx_107_012(
    optimizer,
    objective_tol,
    primal_tol,
    dual_tol,
    termination_target = TERMINATION_TARGET_LOCAL,
    primal_target = PRIMAL_TARGET_LOCAL,
)
    model = Model(optimizer)

    @variable(model, x, start = 1.5)
    @variable(model, y, start = 0.5)

    @objective(model, Min, (x - 1.0)^2 + (y - 1.0)^2)
    @NLconstraint(model, x^2 + y^2 <= 1)

    optimize!(model)

    check_status(model, FEASIBLE_PROBLEM, termination_target, primal_target)
    check_objective(model, 0.17157287363083387, tol = objective_tol)
    return check_solution([x, y], [1 / sqrt(2), 1 / sqrt(2)], tol = primal_tol)
end


text:
include("runtests-prerequisite.jl")

abstract type IMyFantasticStruct <: PostgresORM.IEntity end
abstract type IMyGreatMovie <: PostgresORM.IEntity end

mutable struct MyGreatMovie <: IMyGreatMovie

  title::Union{Missing,String}
  year::Union{Missing,Int32}

  MyGreatMovie(args::NamedTuple) = MyGreatMovie(;args...)
  MyGreatMovie(;title = missing,
           year = missing,) = (
            x = new(missing, missing);
            x.title = title;
            x.year = year;
            return x
            )
end

mutable struct MyFantasticStruct <: IMyFantasticStruct

  id::Union{Missing,Int32}
  first_name::Union{Missing,String}
  favorite_movie::Union{Missing, IMyGreatMovie}

  MyFantasticStruct(args::NamedTuple) = MyFantasticStruct(;args...)
  MyFantasticStruct(;
        id = missing,
        first_name = missing,
        favorite_movie = missing,
 ) = (
            x = new(missing, missing, missing);
            x.id = id;
            x.first_name = first_name;
            x.favorite_movie = favorite_movie;
            return x
            )

end

module MyFantasticStructORM
    data_type = Main.MyFantasticStruct
    Main.PostgresORM.get_orm(x::Main.MyFantasticStruct) = return(MyFantasticStructORM)
    get_table_name() = "public.my_new_struct"
    const columns_selection_and_mapping =
      Dict(
         :id => "id",
         :first_name => "first_name",
         :favorite_movie => ["favorite_movie_title","favorite_movie_year"]
      )
    get_id_props() = return [:id]

    # A dictionnary of mapping between fields symbols and overriding types
    #   Left hanside is the field symbol ; right hand side is the type override
    const types_override = Dict(:favorite_movie => Main.MyGreatMovie)

    const track_changes = false
end

@testset "Test `util_get_property_real_type`" begin
  @test PostgresORM.Controller.util_get_property_real_type(MyFantasticStruct, :favorite_movie) == Main.MyGreatMovie
  @test PostgresORM.Controller.util_get_property_real_type(MyFantasticStruct, :first_name) == String
end


text:
Base.summary(annealing::AbstractAnnealing) = string(
    TYPE_COLOR,
    nameof(typeof(annealing)),
    NO_COLOR,
    " with ",
    TYPE_COLOR,
    typeof(annealing.H),
    NO_COLOR,
    " and u0 ",
    TYPE_COLOR,
    typeof(annealing.u0),
    NO_COLOR,
)

function Base.show(io::IO, A::AbstractAnnealing)
    println(io, summary(A))
    print(io, "u0 size: ")
    show(io, size(A.u0))
end

Base.summary(interaction::Interaction) = string(
    TYPE_COLOR,
    nameof(typeof(interaction)),
    NO_COLOR,
)

function Base.show(io::IO, I::Interaction)
    print(io, summary(I))
    print(io, " with ")
    show(io, I.coupling)
    println(io)
    print(io, "and bath: ")
    show(io, I.bath)
end

Base.summary(interactions::InteractionSet) = string(
    TYPE_COLOR,
    nameof(typeof(interactions)),
    NO_COLOR,
)

function Base.show(io::IO, I::InteractionSet)
    print(io, summary(I))
    print(io, " with ")
    show(io, length(I))
    print(io, " interactions")
end


text:
# This file is part of the Taylor1Series.jl Julia package, MIT license
#
# Luis Benet & David P. Sanders
# UNAM
#
# MIT Expat license
#

## Conversion
convert(::Type{Taylor1{T}}, a::Taylor1) where {T<:Number} =
    Taylor1(convert(Array{T,1}, a.coeffs), a.order)

convert(::Type{Taylor1{T}}, a::Taylor1{T}) where {T<:Number} = a

convert(::Type{Taylor1{Rational{T}}}, a::Taylor1{S}) where
        {T<:Integer, S<:AbstractFloat} = Taylor1(rationalize.(a[:]), a.order)

convert(::Type{Taylor1{T}}, b::Array{T,1}) where {T<:Number} =
    Taylor1(b, length(b)-1)

convert(::Type{Taylor1{T}}, b::Array{S,1}) where {T<:Number, S<:Number} =
    Taylor1(convert(Array{T,1},b), length(b)-1)

convert(::Type{Taylor1{T}}, b::S)  where {T<:Number, S<:Number} =
    Taylor1([convert(T,b)], 0)

convert(::Type{Taylor1{T}}, b::T)  where {T<:Number} = Taylor1([b], 0)

convert(::Type{Taylor1}, a::T) where {T<:Number} = Taylor1(a, 0)

convert(::Type{HomogeneousPolynomial{T}}, a::HomogeneousPolynomial) where {T<:Number} =
    HomogeneousPolynomial(convert(Array{T,1}, a.coeffs), a.order)

convert(::Type{HomogeneousPolynomial{T}}, a::HomogeneousPolynomial{T}) where {T<:Number} =
    a

function convert(::Type{HomogeneousPolynomial{Rational{T}}},
        a::HomogeneousPolynomial{S}) where {T<:Integer, S<:AbstractFloat}

    la = length(a.coeffs)
    v = Array{Rational{T}}(undef, la)
    v .= rationalize.(a[1:la], tol=eps(one(S)))
    return HomogeneousPolynomial(v, a.order)
end

convert(::Type{HomogeneousPolynomial{T}}, b::Array{S,1}) where {T<:Number, S<:Number} =
    HomogeneousPolynomial(convert(Array{T,1}, b), orderH(b))

convert(::Type{HomogeneousPolynomial{T}}, b::S) where {T<:Number, S<:Number}=
    HomogeneousPolynomial([convert(T,b)], 0)

convert(::Type{HomogeneousPolynomial{T}}, b::Array{T,1}) where {T<:Number} =
    HomogeneousPolynomial(b, orderH(b))

convert(::Type{HomogeneousPolynomial{T}}, b::T) where {T<:Number} =
    HomogeneousPolynomial([b], 0)

convert(::Type{HomogeneousPolynomial}, a::Number) = HomogeneousPolynomial([a],0)

convert(::Type{TaylorN{T}}, a::TaylorN) where {T<:Number} =
    TaylorN( convert(Array{HomogeneousPolynomial{T},1}, a.coeffs), a.order)

convert(::Type{TaylorN{T}}, a::TaylorN{T}) where {T<:Number} = a

convert(::Type{TaylorN{T}}, b::HomogeneousPolynomial{S}) where {T<:Number, S<:Number} =
    TaylorN( [convert(HomogeneousPolynomial{T}, b)], b.order)

convert(::Type{TaylorN{T}}, b::Array{HomogeneousPolynomial{S},1}) where {T<:Number, S<:Number} =
    TaylorN( convert(Array{HomogeneousPolynomial{T},1}, b), maxorderH(b))

convert(::Type{TaylorN{T}}, b::S)  where {T<:Number, S<:Number} =
    TaylorN( [HomogeneousPolynomial([convert(T, b)], 0)], 0)

convert(::Type{TaylorN{T}}, b::HomogeneousPolynomial{T}) where {T<:Number} =
    TaylorN( [b], b.order)

convert(::Type{TaylorN{T}}, b::Array{HomogeneousPolynomial{T},1}) where {T<:Number} =
    TaylorN( b, maxorderH(b))

convert(::Type{TaylorN{T}}, b::T) where {T<:Number} =
    TaylorN( [HomogeneousPolynomial([b], 0)], 0)

convert(::Type{TaylorN}, b::Number) = TaylorN( [HomogeneousPolynomial([b], 0)], 0)

function convert(::Type{TaylorN{Taylor1{T}}}, s::Taylor1{TaylorN{T}}) where {T<:NumberNotSeries}

    orderN = get_order()
    r = zeros(HomogeneousPolynomial{Taylor1{T}}, orderN)

    v = zeros(T, s.order+1)
    @inbounds for ordT in eachindex(s)
        v[ordT+1] = one(T)
        @inbounds for ordHP in 0:s[ordT].order
            @inbounds for ic in eachindex(s[ordT][ordHP].coeffs)
                coef = s[ordT][ordHP][ic]
                r[ordHP+1][ic] += Taylor1( coef.*v )
            end
        end
        v[ordT+1] = zero(T)
    end
    return TaylorN(r)
end

function convert(::Type{Taylor1{TaylorN{T}}}, s::TaylorN{Taylor1{T}}) where {T<:NumberNotSeries}

    ordert = 0
    for ordHP in eachindex(s)
        ordert = max(ordert, s[ordHP][1].order)
    end
    vT = Array{TaylorN{T}}(undef, ordert+1)
    @inbounds for ordT in eachindex(vT)
        vT[ordT] = TaylorN(zero(T), s.order)
    end

    @inbounds for ordN in eachindex(s)
        vHP = HomogeneousPolynomial(zeros(T, length(s[ordN])))
        @inbounds for ihp in eachindex(s[ordN].coeffs)
            @inbounds for ind in eachindex(s[ordN][ihp].coeffs)
                c = s[ordN][ihp][ind-1]
                vHP[ihp] = c
                vT[ind] += TaylorN(vHP, s.order)
                vHP[ihp] = zero(T)
            end
        end
    end
    return Taylor1(vT)
end

function convert(::Type{Array{TaylorN{Taylor1{T}},N}},
        s::Array{Taylor1{TaylorN{T}},N}) where {T<:NumberNotSeries, N}

    v = Array{TaylorN{Taylor1{T}}}(undef, size(s))
    @simd for ind in eachindex(s)
        @inbounds v[ind] = convert(TaylorN{Taylor1{T}}, s[ind])
    end
    return v
end

function convert(::Type{Array{Taylor1{TaylorN{T}},N}},
        s::Array{TaylorN{Taylor1{T}},N}) where {T<:NumberNotSeries, N}

    v = Array{Taylor1{TaylorN{T}}}(undef, size(s))
    @simd for ind in eachindex(s)
        @inbounds v[ind] = convert(Taylor1{TaylorN{T}}, s[ind])
    end
    return v
end

# Promotion
promote_rule(::Type{Taylor1{T}}, ::Type{Taylor1{T}}) where {T<:Number} = Taylor1{T}

promote_rule(::Type{Taylor1{T}}, ::Type{Taylor1{S}}) where {T<:Number, S<:Number} =
    Taylor1{promote_type(T,S)}

promote_rule(::Type{Taylor1{T}}, ::Type{Array{T,1}}) where {T<:Number} = Taylor1{T}

promote_rule(::Type{Taylor1{T}}, ::Type{Array{S,1}}) where {T<:Number, S<:Number} =
    Taylor1{promote_type(T,S)}

promote_rule(::Type{Taylor1{T}}, ::Type{T}) where {T<:Number} = Taylor1{T}

promote_rule(::Type{Taylor1{T}}, ::Type{S}) where {T<:Number, S<:Number} =
    Taylor1{promote_type(T,S)}

promote_rule(::Type{Taylor1{Taylor1{T}}}, ::Type{Taylor1{T}}) where {T<:Number} =
    Taylor1{Taylor1{T}}

promote_rule(::Type{HomogeneousPolynomial{T}},
    ::Type{HomogeneousPolynomial{S}}) where {T<:Number, S<:Number} =
        HomogeneousPolynomial{promote_type(T,S)}

promote_rule(::Type{HomogeneousPolynomial{T}},
    ::Type{HomogeneousPolynomial{T}}) where {T<:Number} = HomogeneousPolynomial{T}

promote_rule(::Type{HomogeneousPolynomial{T}},
    ::Type{Array{S,1}}) where {T<:Number, S<:Number} = HomogeneousPolynomial{promote_type(T,S)}

promote_rule(::Type{HomogeneousPolynomial{T}}, ::Type{S}) where
    {T<:Number, S<:NumberNotSeries} = HomogeneousPolynomial{promote_type(T,S)}

promote_rule(::Type{TaylorN{T}}, ::Type{TaylorN{S}}) where {T<:Number, S<:Number}=
    TaylorN{promote_type(T,S)}

promote_rule(::Type{TaylorN{T}}, ::Type{TaylorN{T}}) where {T<:Number} = TaylorN{T}

promote_rule(::Type{TaylorN{T}}, ::Type{HomogeneousPolynomial{S}}) where
    {T<:Number, S<:Number} = TaylorN{promote_type(T,S)}

promote_rule(::Type{TaylorN{T}}, ::Type{Array{HomogeneousPolynomial{S},1}}) where
    {T<:Number, S<:Number} = TaylorN{promote_type(T,S)}

promote_rule(::Type{TaylorN{T}}, ::Type{S}) where {T<:Number, S<:Number} =
    TaylorN{promote_type(T,S)}

# Order may matter
promote_rule(::Type{S}, ::Type{T}) where {S<:NumberNotSeries, T<:AbstractSeries} =
    promote_rule(T,S)

promote_rule(::Type{S}, ::Type{T}) where
    {S<:AbstractIrrational, T<:AbstractSeries} = promote_rule(T,S)

promote_rule(::Type{Taylor1{T}}, ::Type{TaylorN{S}}) where {T<:NumberNotSeries, S<:NumberNotSeries} =
    throw(ArgumentError("There is no reasonable promotion among `Taylor1{$T}` and `TaylorN{$S}` types"))

# Nested Taylor1's
function promote(a::Taylor1{Taylor1{T}}, b::Taylor1{T}) where {T<:NumberNotSeriesN}
    order_a = get_order(a)
    order_b = get_order(b)
    zb = zero(b)
    new_bcoeffs = similar(a.coeffs)
    new_bcoeffs[1] = b
    @inbounds for ind in 2:order_a+1
        new_bcoeffs[ind] = zb
    end
    return a, Taylor1(b, order_a)
end
promote(b::Taylor1{T}, a::Taylor1{Taylor1{T}}) where {T<:NumberNotSeriesN} =
    reverse(promote(a, b))


text:
function spawnrandompoint(randFloatX::Float64=rand(), randFloatY::Float64=rand())::Array{Float64,1}
    newxy = copy(SPAWN_TOP_LEFT_XY)
    newxy[1] += randFloatX * SPAWN_X_RANGE
    newxy[2] -= randFloatY * SPAWN_Y_RANGE
    return newxy
end

function insensor(xy::Array{Float64,1})::Bool
    y = xy[2]
    return y <= PPD_MIN_Y
end

function inwalls(x::Float64, y::Float64)::Bool
    withinx = SENSOR_HALF_WIDTH < abs(x) <= ESCAPE_X
    withiny = PPD_MIN_Y < y <= WALL_Y
    return withinx && withiny
end

function inppd(xy::Array{Float64,1})::Bool
    withinx = abs(xy[1]) <= SENSOR_HALF_WIDTH
    withiny = PPD_MIN_Y < xy[2] <= PPD_MAX_Y
    return withinx && withiny
end

function inwater(xy::Array{Float64,1})::Bool
    x, y = xy
    abovewallx = abs(x) <= ESCAPE_X
    abovewally = WALL_Y < y <= ESCAPE_Y
    inwellx = abs(x) <= SENSOR_HALF_WIDTH
    inwelly = PPD_MAX_Y < y <= WALL_Y
    return (abovewallx && abovewally) || (inwellx && inwelly)
end

function inenz(xy::Array{Float64,1})::Bool
    # No enzyme layer
    return false
end

function whereoutsidespawn(xy::Array{Float64,1})::String
    # No Enzyme Layer, only horizontal inteface between
    # Water and PPD Layer inside sensor well
    return "N"
end

function sensorcases(x::Float64)::String
    # only one sensor
    return "sensor"
end

function wallcases(initxy::Array{Float64,1},
                   dx::Float64,
                   dy::Float64,
                   ending_step_size::Float64)
    return approach_wall!(initxy, dx, dy, ending_step_size), "no collision"
end

max_stars_repo_name: s-y-wu/HMCResearchRandomWalks

max_issues_repo_name: s-y-wu/HMCResearchRandomWalks

max_forks_repo_name: s-y-wu/HMCResearchRandomWalks


text:
@testset "minimal_is_empty_test" begin

    @test isempty(emptyinterval()) === true

    @test isempty(interval(-Inf,+Inf)) === false

    @test isempty(interval(1.0,2.0)) === false

    @test isempty(interval(-1.0,2.0)) === false

    @test isempty(interval(-3.0,-2.0)) === false

    @test isempty(interval(-Inf,2.0)) === false

    @test isempty(interval(-Inf,0.0)) === false

    @test isempty(interval(-Inf,-0.0)) === false

    @test isempty(interval(0.0,Inf)) === false

    @test isempty(interval(-0.0,Inf)) === false

    @test isempty(interval(-0.0,0.0)) === false

    @test isempty(interval(0.0,-0.0)) === false

    @test isempty(interval(0.0,0.0)) === false

    @test isempty(interval(-0.0,-0.0)) === false

end

@testset "minimal_is_empty_dec_test" begin

    @test isempty(nai()) === false

    @test isempty(DecoratedInterval(emptyinterval(), trv)) === true

    @test isempty(DecoratedInterval(interval(-Inf,+Inf), def)) === false

    @test isempty(DecoratedInterval(interval(1.0,2.0), com)) === false

    @test isempty(DecoratedInterval(interval(-1.0,2.0), trv)) === false

    @test isempty(DecoratedInterval(interval(-3.0,-2.0), dac)) === false

    @test isempty(DecoratedInterval(interval(-Inf,2.0), trv)) === false

    @test isempty(DecoratedInterval(interval(-Inf,0.0), trv)) === false

    @test isempty(DecoratedInterval(interval(-Inf,-0.0), trv)) === false

    @test isempty(DecoratedInterval(interval(0.0,Inf), def)) === false

    @test isempty(DecoratedInterval(interval(-0.0,Inf), trv)) === false

    @test isempty(DecoratedInterval(interval(-0.0,0.0), com)) === false

    @test isempty(DecoratedInterval(interval(0.0,-0.0), trv)) === false

    @test isempty(DecoratedInterval(interval(0.0,0.0), trv)) === false

    @test isempty(DecoratedInterval(interval(-0.0,-0.0), trv)) === false

end

@testset "minimal_is_entire_test" begin

    @test isentire(emptyinterval()) === false

    @test isentire(interval(-Inf,+Inf)) === true

    @test isentire(interval(1.0,2.0)) === false

    @test isentire(interval(-1.0,2.0)) === false

    @test isentire(interval(-3.0,-2.0)) === false

    @test isentire(interval(-Inf,2.0)) === false

    @test isentire(interval(-Inf,0.0)) === false

    @test isentire(interval(-Inf,-0.0)) === false

    @test isentire(interval(0.0,Inf)) === false

    @test isentire(interval(-0.0,Inf)) === false

    @test isentire(interval(-0.0,0.0)) === false

    @test isentire(interval(0.0,-0.0)) === false

    @test isentire(interval(0.0,0.0)) === false

    @test isentire(interval(-0.0,-0.0)) === false

end

@testset "minimal_is_entire_dec_test" begin

    @test isentire(nai()) === false

    @test isentire(DecoratedInterval(emptyinterval(), trv)) === false

    @test isentire(DecoratedInterval(interval(-Inf,+Inf), trv)) === true

    @test isentire(DecoratedInterval(interval(-Inf,+Inf), def)) === true

    @test isentire(DecoratedInterval(interval(-Inf,+Inf), dac)) === true

    @test isentire(DecoratedInterval(interval(1.0,2.0), com)) === false

    @test isentire(DecoratedInterval(interval(-1.0,2.0), trv)) === false

    @test isentire(DecoratedInterval(interval(-3.0,-2.0), dac)) === false

    @test isentire(DecoratedInterval(interval(-Inf,2.0), trv)) === false

    @test isentire(DecoratedInterval(interval(-Inf,0.0), trv)) === false

    @test isentire(DecoratedInterval(interval(-Inf,-0.0), trv)) === false

    @test isentire(DecoratedInterval(interval(0.0,Inf), def)) === false

    @test isentire(DecoratedInterval(interval(-0.0,Inf), trv)) === false

    @test isentire(DecoratedInterval(interval(-0.0,0.0), com)) === false

    @test isentire(DecoratedInterval(interval(0.0,-0.0), trv)) === false

    @test isentire(DecoratedInterval(interval(0.0,0.0), trv)) === false

    @test isentire(DecoratedInterval(interval(-0.0,-0.0), trv)) === false

end

@testset "minimal_is_nai_dec_test" begin

    @test isnai(nai()) === true

    @test isnai(DecoratedInterval(interval(-Inf,+Inf), trv)) === false

    @test isnai(DecoratedInterval(interval(-Inf,+Inf), def)) === false

    @test isnai(DecoratedInterval(interval(-Inf,+Inf), dac)) === false

    @test isnai(DecoratedInterval(interval(1.0,2.0), com)) === false

    @test isnai(DecoratedInterval(interval(-1.0,2.0), trv)) === false

    @test isnai(DecoratedInterval(interval(-3.0,-2.0), dac)) === false

    @test isnai(DecoratedInterval(interval(-Inf,2.0), trv)) === false

    @test isnai(DecoratedInterval(interval(-Inf,0.0), trv)) === false

    @test isnai(DecoratedInterval(interval(-Inf,-0.0), trv)) === false

    @test isnai(DecoratedInterval(interval(0.0,Inf), def)) === false

    @test isnai(DecoratedInterval(interval(-0.0,Inf), trv)) === false

    @test isnai(DecoratedInterval(interval(-0.0,0.0), com)) === false

    @test isnai(DecoratedInterval(interval(0.0,-0.0), trv)) === false

    @test isnai(DecoratedInterval(interval(0.0,0.0), trv)) === false

    @test isnai(DecoratedInterval(interval(-0.0,-0.0), trv)) === false

end

@testset "minimal_equal_test" begin

    @test ==(interval(1.0,2.0), interval(1.0,2.0)) === true

    @test ==(interval(1.0,2.1), interval(1.0,2.0)) === false

    @test ==(emptyinterval(), emptyinterval()) === true

    @test ==(emptyinterval(), interval(1.0,2.0)) === false

    @test ==(interval(-Inf,+Inf), interval(-Inf,+Inf)) === true

    @test ==(interval(1.0,2.4), interval(-Inf,+Inf)) === false

    @test ==(interval(1.0,Inf), interval(1.0,Inf)) === true

    @test ==(interval(1.0,2.4), interval(1.0,Inf)) === false

    @test ==(interval(-Inf,2.0), interval(-Inf,2.0)) === true

    @test ==(interval(-Inf,2.4), interval(-Inf,2.0)) === false

    @test ==(interval(-2.0,0.0), interval(-2.0,0.0)) === true

    @test ==(interval(-0.0,2.0), interval(0.0,2.0)) === true

    @test ==(interval(-0.0,-0.0), interval(0.0,0.0)) === true

    @test ==(interval(-0.0,0.0), interval(0.0,0.0)) === true

    @test ==(interval(0.0,-0.0), interval(0.0,0.0)) === true

end

@testset "minimal_equal_dec_test" begin

    @test ==(DecoratedInterval(interval(1.0,2.0), def), DecoratedInterval(interval(1.0,2.0), trv)) === true

    @test ==(DecoratedInterval(interval(1.0,2.1), trv), DecoratedInterval(interval(1.0,2.0), trv)) === false

    @test ==(DecoratedInterval(emptyinterval(), trv), DecoratedInterval(emptyinterval(), trv)) === true

    @test ==(DecoratedInterval(emptyinterval(), trv), nai()) === false

    @test ==(nai(), DecoratedInterval(emptyinterval(), trv)) === false

    @test ==(nai(), nai()) === false

    @test ==(DecoratedInterval(emptyinterval(), trv), DecoratedInterval(interval(1.0,2.0), trv)) === false

    @test ==(nai(), DecoratedInterval(interval(1.0,2.0), trv)) === false

    @test ==(DecoratedInterval(interval(-Inf,+Inf), def), DecoratedInterval(interval(-Inf,+Inf), trv)) === true

    @test ==(DecoratedInterval(interval(1.0,2.4), trv), DecoratedInterval(interval(-Inf,+Inf), trv)) === false

    @test ==(DecoratedInterval(interval(1.0,Inf), trv), DecoratedInterval(interval(1.0,Inf), trv)) === true

    @test ==(DecoratedInterval(interval(1.0,2.4), def), DecoratedInterval(interval(1.0,Inf), trv)) === false

    @test ==(DecoratedInterval(interval(-Inf,2.0), trv), DecoratedInterval(interval(-Inf,2.0), trv)) === true

    @test ==(DecoratedInterval(interval(-Inf,2.4), def), DecoratedInterval(interval(-Inf,2.0), trv)) === false

    @test ==(DecoratedInterval(interval(-2.0,0.0), trv), DecoratedInterval(interval(-2.0,0.0), trv)) === true

    @test ==(DecoratedInterval(interval(-0.0,2.0), def), DecoratedInterval(interval(0.0,2.0), trv)) === true

    @test ==(DecoratedInterval(interval(-0.0,-0.0), trv), DecoratedInterval(interval(0.0,0.0), trv)) === true

    @test ==(DecoratedInterval(interval(-0.0,0.0), def), DecoratedInterval(interval(0.0,0.0), trv)) === true

    @test ==(DecoratedInterval(interval(0.0,-0.0), trv), DecoratedInterval(interval(0.0,0.0), trv)) === true

end

@testset "minimal_subset_test" begin

    @test issubset(emptyinterval(), emptyinterval()) === true

    @test issubset(emptyinterval(), interval(0.0,4.0)) === true

    @test issubset(emptyinterval(), interval(-0.0,4.0)) === true

    @test issubset(emptyinterval(), interval(-0.1,1.0)) === true

    @test issubset(emptyinterval(), interval(-0.1,0.0)) === true

    @test issubset(emptyinterval(), interval(-0.1,-0.0)) === true

    @test issubset(emptyinterval(), interval(-Inf,+Inf)) === true

    @test issubset(interval(0.0,4.0), emptyinterval()) === false

    @test issubset(interval(-0.0,4.0), emptyinterval()) === false

    @test issubset(interval(-0.1,1.0), emptyinterval()) === false

    @test issubset(interval(-Inf,+Inf), emptyinterval()) === false

    @test issubset(interval(0.0,4.0), interval(-Inf,+Inf)) === true

    @test issubset(interval(-0.0,4.0), interval(-Inf,+Inf)) === true

    @test issubset(interval(-0.1,1.0), interval(-Inf,+Inf)) === true

    @test issubset(interval(-Inf,+Inf), interval(-Inf,+Inf)) === true

    @test issubset(interval(1.0,2.0), interval(1.0,2.0)) === true

    @test issubset(interval(1.0,2.0), interval(0.0,4.0)) === true

    @test issubset(interval(1.0,2.0), interval(-0.0,4.0)) === true

    @test issubset(interval(0.1,0.2), interval(0.0,4.0)) === true

    @test issubset(interval(0.1,0.2), interval(-0.0,4.0)) === true

    @test issubset(interval(-0.1,-0.1), interval(-4.0, 3.4)) === true

    @test issubset(interval(0.0,0.0), interval(-0.0,-0.0)) === true

    @test issubset(interval(-0.0,-0.0), interval(0.0,0.0)) === true

    @test issubset(interval(-0.0,0.0), interval(0.0,0.0)) === true

    @test issubset(interval(-0.0,0.0), interval(0.0,-0.0)) === true

    @test issubset(interval(0.0,-0.0), interval(0.0,0.0)) === true

    @test issubset(interval(0.0,-0.0), interval(-0.0,0.0)) === true

end

@testset "minimal_subset_dec_test" begin

    @test issubset(nai(), nai()) === false

    @test issubset(nai(), DecoratedInterval(emptyinterval(), trv)) === false

    @test issubset(nai(), DecoratedInterval(interval(0.0,4.0), trv)) === false

    @test issubset(DecoratedInterval(emptyinterval(), trv), DecoratedInterval(interval(0.0,4.0), trv)) === true

    @test issubset(DecoratedInterval(emptyinterval(), trv), DecoratedInterval(interval(-0.0,4.0), def)) === true

    @test issubset(DecoratedInterval(emptyinterval(), trv), DecoratedInterval(interval(-0.1,1.0), trv)) === true

    @test issubset(DecoratedInterval(emptyinterval(), trv), DecoratedInterval(interval(-0.1,0.0), trv)) === true

    @test issubset(DecoratedInterval(emptyinterval(), trv), DecoratedInterval(interval(-0.1,-0.0), trv)) === true

    @test issubset(DecoratedInterval(emptyinterval(), trv), DecoratedInterval(interval(-Inf,+Inf), trv)) === true

    @test issubset(DecoratedInterval(interval(0.0,4.0), trv), DecoratedInterval(emptyinterval(), trv)) === false

    @test issubset(DecoratedInterval(interval(-0.0,4.0), def), DecoratedInterval(emptyinterval(), trv)) === false

    @test issubset(DecoratedInterval(interval(-0.1,1.0), trv), DecoratedInterval(emptyinterval(), trv)) === false

    @test issubset(DecoratedInterval(interval(-Inf,+Inf), trv), DecoratedInterval(emptyinterval(), trv)) === false

    @test issubset(DecoratedInterval(interval(0.0,4.0), trv), DecoratedInterval(interval(-Inf,+Inf), trv)) === true

    @test issubset(DecoratedInterval(interval(-0.0,4.0), trv), DecoratedInterval(interval(-Inf,+Inf), trv)) === true

    @test issubset(DecoratedInterval(interval(-0.1,1.0), trv), DecoratedInterval(interval(-Inf,+Inf), trv)) === true

    @test issubset(DecoratedInterval(interval(-Inf,+Inf), trv), DecoratedInterval(interval(-Inf,+Inf), trv)) === true

    @test issubset(DecoratedInterval(interval(1.0,2.0), trv), DecoratedInterval(interval(1.0,2.0), trv)) === true

    @test issubset(DecoratedInterval(interval(1.0,2.0), trv), DecoratedInterval(interval(0.0,4.0), trv)) === true

    @test issubset(DecoratedInterval(interval(1.0,2.0), def), DecoratedInterval(interval(-0.0,4.0), def)) === true

    @test issubset(DecoratedInterval(interval(0.1,0.2), trv), DecoratedInterval(interval(0.0,4.0), trv)) === true

    @test issubset(DecoratedInterval(interval(0.1,0.2), trv), DecoratedInterval(interval(-0.0,4.0), def)) === true

    @test issubset(DecoratedInterval(interval(-0.1,-0.1), trv), DecoratedInterval(interval(-4.0, 3.4), trv)) === true

    @test issubset(DecoratedInterval(interval(0.0,0.0), trv), DecoratedInterval(interval(-0.0,-0.0), trv)) === true

    @test issubset(DecoratedInterval(interval(-0.0,-0.0), trv), DecoratedInterval(interval(0.0,0.0), def)) === true

    @test issubset(DecoratedInterval(interval(-0.0,0.0), trv), DecoratedInterval(interval(0.0,0.0), trv)) === true

    @test issubset(DecoratedInterval(interval(-0.0,0.0), trv), DecoratedInterval(interval(0.0,-0.0), trv)) === true

    @test issubset(DecoratedInterval(interval(0.0,-0.0), def), DecoratedInterval(interval(0.0,0.0), trv)) === true

    @test issubset(DecoratedInterval(interval(0.0,-0.0), trv), DecoratedInterval(interval(-0.0,0.0), trv)) === true

end

@testset "minimal_less_test" begin

    @test ≤(emptyinterval(), emptyinterval()) === true

    @test ≤(interval(1.0,2.0), emptyinterval()) === false

    @test ≤(emptyinterval(), interval(1.0,2.0)) === false

    @test ≤(interval(-Inf,+Inf), interval(-Inf,+Inf)) === true

    @test ≤(interval(1.0,2.0), interval(-Inf,+Inf)) === false

    @test ≤(interval(0.0,2.0), interval(-Inf,+Inf)) === false

    @test ≤(interval(-0.0,2.0), interval(-Inf,+Inf)) === false

    @test ≤(interval(-Inf,+Inf), interval(1.0,2.0)) === false

    @test ≤(interval(-Inf,+Inf), interval(0.0,2.0)) === false

    @test ≤(interval(-Inf,+Inf), interval(-0.0,2.0)) === false

    @test ≤(interval(0.0,2.0), interval(0.0,2.0)) === true

    @test ≤(interval(0.0,2.0), interval(-0.0,2.0)) === true

    @test ≤(interval(0.0,2.0), interval(1.0,2.0)) === true

    @test ≤(interval(-0.0,2.0), interval(1.0,2.0)) === true

    @test ≤(interval(1.0,2.0), interval(1.0,2.0)) === true

    @test ≤(interval(1.0,2.0), interval(3.0,4.0)) === true

    @test ≤(interval(1.0,3.5), interval(3.0,4.0)) === true

    @test ≤(interval(1.0,4.0), interval(3.0,4.0)) === true

    @test ≤(interval(-2.0,-1.0), interval(-2.0,-1.0)) === true

    @test ≤(interval(-3.0,-1.5), interval(-2.0,-1.0)) === true

    @test ≤(interval(0.0,0.0), interval(-0.0,-0.0)) === true

    @test ≤(interval(-0.0,-0.0), interval(0.0,0.0)) === true

    @test ≤(interval(-0.0,0.0), interval(0.0,0.0)) === true

    @test ≤(interval(-0.0,0.0), interval(0.0,-0.0)) === true

    @test ≤(interval(0.0,-0.0), interval(0.0,0.0)) === true

    @test ≤(interval(0.0,-0.0), interval(-0.0,0.0)) === true

end

@testset "minimal_less_dec_test" begin

    @test ≤(nai(), nai()) === false

    @test ≤(DecoratedInterval(emptyinterval(), trv), nai()) === false

    @test ≤(DecoratedInterval(interval(1.0,2.0), trv), nai()) === false

    @test ≤(nai(), DecoratedInterval(interval(1.0,2.0), def)) === false

    @test ≤(DecoratedInterval(emptyinterval(), trv), DecoratedInterval(emptyinterval(), trv)) === true

    @test ≤(DecoratedInterval(interval(1.0,2.0), trv), DecoratedInterval(emptyinterval(), trv)) === false

    @test ≤(DecoratedInterval(emptyinterval(), trv), DecoratedInterval(interval(1.0,2.0), trv)) === false

    @test ≤(DecoratedInterval(interval(-Inf,+Inf), trv), DecoratedInterval(interval(-Inf,+Inf), trv)) === true

    @test ≤(DecoratedInterval(interval(1.0,2.0), def), DecoratedInterval(interval(-Inf,+Inf), trv)) === false

    @test ≤(DecoratedInterval(interval(0.0,2.0), trv), DecoratedInterval(interval(-Inf,+Inf), trv)) === false

    @test ≤(DecoratedInterval(interval(-0.0,2.0), trv), DecoratedInterval(interval(-Inf,+Inf), trv)) === false

    @test ≤(DecoratedInterval(interval(-Inf,+Inf), trv), DecoratedInterval(interval(1.0,2.0), trv)) === false

    @test ≤(DecoratedInterval(interval(-Inf,+Inf), trv), DecoratedInterval(interval(0.0,2.0), def)) === false

    @test ≤(DecoratedInterval(interval(-Inf,+Inf), trv), DecoratedInterval(interval(-0.0,2.0), trv)) === false

    @test ≤(DecoratedInterval(interval(0.0,2.0), trv), DecoratedInterval(interval(0.0,2.0), trv)) === true

    @test ≤(DecoratedInterval(interval(0.0,2.0), trv), DecoratedInterval(interval(-0.0,2.0), trv)) === true

    @test ≤(DecoratedInterval(interval(0.0,2.0), def), DecoratedInterval(interval(1.0,2.0), def)) === true

    @test ≤(DecoratedInterval(interval(-0.0,2.0), trv), DecoratedInterval(interval(1.0,2.0), trv)) === true

    @test ≤(DecoratedInterval(interval(1.0,2.0), trv), DecoratedInterval(interval(1.0,2.0), trv)) === true

    @test ≤(DecoratedInterval(interval(1.0,2.0), trv), DecoratedInterval(interval(3.0,4.0), def)) === true

    @test ≤(DecoratedInterval(interval(1.0,3.5), trv), DecoratedInterval(interval(3.0,4.0), trv)) === true

    @test ≤(DecoratedInterval(interval(1.0,4.0), trv), DecoratedInterval(interval(3.0,4.0), trv)) === true

    @test ≤(DecoratedInterval(interval(-2.0,-1.0), trv), DecoratedInterval(interval(-2.0,-1.0), trv)) === true

    @test ≤(DecoratedInterval(interval(-3.0,-1.5), trv), DecoratedInterval(interval(-2.0,-1.0), trv)) === true

    @test ≤(DecoratedInterval(interval(0.0,0.0), trv), DecoratedInterval(interval(-0.0,-0.0), trv)) === true

    @test ≤(DecoratedInterval(interval(-0.0,-0.0), trv), DecoratedInterval(interval(0.0,0.0), def)) === true

    @test ≤(DecoratedInterval(interval(-0.0,0.0), trv), DecoratedInterval(interval(0.0,0.0), trv)) === true

    @test ≤(DecoratedInterval(interval(-0.0,0.0), trv), DecoratedInterval(interval(0.0,-0.0), trv)) === true

    @test ≤(DecoratedInterval(interval(0.0,-0.0), def), DecoratedInterval(interval(0.0,0.0), trv)) === true

    @test ≤(DecoratedInterval(interval(0.0,-0.0), trv), DecoratedInterval(interval(-0.0,0.0), trv)) === true

end

@testset "minimal_precedes_test" begin

    @test precedes(emptyinterval(), interval(3.0,4.0)) === true

    @test precedes(interval(3.0,4.0), emptyinterval()) === true

    @test precedes(emptyinterval(), emptyinterval()) === true

    @test precedes(interval(1.0,2.0), interval(-Inf,+Inf)) === false

    @test precedes(interval(0.0,2.0), interval(-Inf,+Inf)) === false

    @test precedes(interval(-0.0,2.0), interval(-Inf,+Inf)) === false

    @test precedes(interval(-Inf,+Inf), interval(1.0,2.0)) === false

    @test precedes(interval(-Inf,+Inf), interval(-Inf,+Inf)) === false

    @test precedes(interval(1.0,2.0), interval(3.0,4.0)) === true

    @test precedes(interval(1.0,3.0), interval(3.0,4.0)) === true

    @test precedes(interval(-3.0, -1.0), interval(-1.0,0.0)) === true

    @test precedes(interval(-3.0, -1.0), interval(-1.0,-0.0)) === true

    @test precedes(interval(1.0,3.5), interval(3.0,4.0)) === false

    @test precedes(interval(1.0,4.0), interval(3.0,4.0)) === false

    @test precedes(interval(-3.0, -0.1), interval(-1.0,0.0)) === false

    @test precedes(interval(0.0,0.0), interval(-0.0,-0.0)) === true

    @test precedes(interval(-0.0,-0.0), interval(0.0,0.0)) === true

    @test precedes(interval(-0.0,0.0), interval(0.0,0.0)) === true

    @test precedes(interval(-0.0,0.0), interval(0.0,-0.0)) === true

    @test precedes(interval(0.0,-0.0), interval(0.0,0.0)) === true

    @test precedes(interval(0.0,-0.0), interval(-0.0,0.0)) === true

end

@testset "minimal_precedes_dec_test" begin

    @test precedes(nai(), DecoratedInterval(interval(3.0,4.0), def)) === false

    @test precedes(DecoratedInterval(interval(3.0,4.0), trv), nai()) === false

    @test precedes(nai(), DecoratedInterval(emptyinterval(), trv)) === false

    @test precedes(nai(), nai()) === false

    @test precedes(DecoratedInterval(emptyinterval(), trv), DecoratedInterval(interval(3.0,4.0), def)) === true

    @test precedes(DecoratedInterval(interval(3.0,4.0), trv), DecoratedInterval(emptyinterval(), trv)) === true

    @test precedes(DecoratedInterval(emptyinterval(), trv), DecoratedInterval(emptyinterval(), trv)) === true

    @test precedes(DecoratedInterval(interval(1.0,2.0), trv), DecoratedInterval(interval(-Inf,+Inf), trv)) === false

    @test precedes(DecoratedInterval(interval(0.0,2.0), trv), DecoratedInterval(interval(-Inf,+Inf), trv)) === false

    @test precedes(DecoratedInterval(interval(-0.0,2.0), trv), DecoratedInterval(interval(-Inf,+Inf), trv)) === false

    @test precedes(DecoratedInterval(interval(-Inf,+Inf), trv), DecoratedInterval(interval(1.0,2.0), trv)) === false

    @test precedes(DecoratedInterval(interval(-Inf,+Inf), trv), DecoratedInterval(interval(-Inf,+Inf), trv)) === false

    @test precedes(DecoratedInterval(interval(1.0,2.0), trv), DecoratedInterval(interval(3.0,4.0), trv)) === true

    @test precedes(DecoratedInterval(interval(1.0,3.0), trv), DecoratedInterval(interval(3.0,4.0), def)) === true

    @test precedes(DecoratedInterval(interval(-3.0, -1.0), def), DecoratedInterval(interval(-1.0,0.0), trv)) === true

    @test precedes(DecoratedInterval(interval(-3.0, -1.0), trv), DecoratedInterval(interval(-1.0,-0.0), trv)) === true

    @test precedes(DecoratedInterval(interval(1.0,3.5), trv), DecoratedInterval(interval(3.0,4.0), trv)) === false

    @test precedes(DecoratedInterval(interval(1.0,4.0), trv), DecoratedInterval(interval(3.0,4.0), trv)) === false

    @test precedes(DecoratedInterval(interval(-3.0, -0.1), trv), DecoratedInterval(interval(-1.0,0.0), trv)) === false

    @test precedes(DecoratedInterval(interval(0.0,0.0), trv), DecoratedInterval(interval(-0.0,-0.0), trv)) === true

    @test precedes(DecoratedInterval(interval(-0.0,-0.0), trv), DecoratedInterval(interval(0.0,0.0), def)) === true

    @test precedes(DecoratedInterval(interval(-0.0,0.0), trv), DecoratedInterval(interval(0.0,0.0), trv)) === true

    @test precedes(DecoratedInterval(interval(-0.0,0.0), def), DecoratedInterval(interval(0.0,-0.0), trv)) === true

    @test precedes(DecoratedInterval(interval(0.0,-0.0), trv), DecoratedInterval(interval(0.0,0.0), trv)) === true

    @test precedes(DecoratedInterval(interval(0.0,-0.0), trv), DecoratedInterval(interval(-0.0,0.0), trv)) === true

end

@testset "minimal_interior_test" begin

    @test isinterior(emptyinterval(), emptyinterval()) === true

    @test isinterior(emptyinterval(), interval(0.0,4.0)) === true

    @test isinterior(interval(0.0,4.0), emptyinterval()) === false

    @test isinterior(interval(-Inf,+Inf), interval(-Inf,+Inf)) === true

    @test isinterior(interval(0.0,4.0), interval(-Inf,+Inf)) === true

    @test isinterior(emptyinterval(), interval(-Inf,+Inf)) === true

    @test isinterior(interval(-Inf,+Inf), interval(0.0,4.0)) === false

    @test isinterior(interval(0.0,4.0), interval(0.0,4.0)) === false

    @test isinterior(interval(1.0,2.0), interval(0.0,4.0)) === true

    @test isinterior(interval(-2.0,2.0), interval(-2.0,4.0)) === false

    @test isinterior(interval(-0.0,-0.0), interval(-2.0,4.0)) === true

    @test isinterior(interval(0.0,0.0), interval(-2.0,4.0)) === true

    @test isinterior(interval(0.0,0.0), interval(-0.0,-0.0)) === false

    @test isinterior(interval(0.0,4.4), interval(0.0,4.0)) === false

    @test isinterior(interval(-1.0,-1.0), interval(0.0,4.0)) === false

    @test isinterior(interval(2.0,2.0), interval(-2.0,-1.0)) === false

end

@testset "minimal_interior_dec_test" begin

    @test isinterior(nai(), nai()) === false

    @test isinterior(nai(), DecoratedInterval(emptyinterval(), trv)) === false

    @test isinterior(nai(), DecoratedInterval(interval(0.0,4.0), trv)) === false

    @test isinterior(DecoratedInterval(interval(0.0,4.0), def), nai()) === false

    @test isinterior(DecoratedInterval(emptyinterval(), trv), DecoratedInterval(emptyinterval(), trv)) === true

    @test isinterior(DecoratedInterval(emptyinterval(), trv), DecoratedInterval(interval(0.0,4.0), trv)) === true

    @test isinterior(DecoratedInterval(interval(0.0,4.0), def), DecoratedInterval(emptyinterval(), trv)) === false

    @test isinterior(DecoratedInterval(interval(-Inf,+Inf), trv), DecoratedInterval(interval(-Inf,+Inf), trv)) === true

    @test isinterior(DecoratedInterval(interval(0.0,4.0), trv), DecoratedInterval(interval(-Inf,+Inf), trv)) === true

    @test isinterior(DecoratedInterval(emptyinterval(), trv), DecoratedInterval(interval(-Inf,+Inf), trv)) === true

    @test isinterior(DecoratedInterval(interval(-Inf,+Inf), trv), DecoratedInterval(interval(0.0,4.0), trv)) === false

    @test isinterior(DecoratedInterval(interval(0.0,4.0), trv), DecoratedInterval(interval(0.0,4.0), trv)) === false

    @test isinterior(DecoratedInterval(interval(1.0,2.0), def), DecoratedInterval(interval(0.0,4.0), trv)) === true

    @test isinterior(DecoratedInterval(interval(-2.0,2.0), trv), DecoratedInterval(interval(-2.0,4.0), def)) === false

    @test isinterior(DecoratedInterval(interval(-0.0,-0.0), trv), DecoratedInterval(interval(-2.0,4.0), trv)) === true

    @test isinterior(DecoratedInterval(interval(0.0,0.0), def), DecoratedInterval(interval(-2.0,4.0), trv)) === true

    @test isinterior(DecoratedInterval(interval(0.0,0.0), trv), DecoratedInterval(interval(-0.0,-0.0), trv)) === false

    @test isinterior(DecoratedInterval(interval(0.0,4.4), trv), DecoratedInterval(interval(0.0,4.0), trv)) === false

    @test isinterior(DecoratedInterval(interval(-1.0,-1.0), trv), DecoratedInterval(interval(0.0,4.0), def)) === false

    @test isinterior(DecoratedInterval(interval(2.0,2.0), def), DecoratedInterval(interval(-2.0,-1.0), trv)) === false

end

@testset "minimal_strictly_less_test" begin

    @test <(emptyinterval(), emptyinterval()) === true

    @test <(interval(1.0,2.0), emptyinterval()) === false

    @test <(emptyinterval(), interval(1.0,2.0)) === false

    @test <(interval(-Inf,+Inf), interval(-Inf,+Inf)) === true

    @test <(interval(1.0,2.0), interval(-Inf,+Inf)) === false

    @test <(interval(-Inf,+Inf), interval(1.0,2.0)) === false

    @test <(interval(1.0,2.0), interval(1.0,2.0)) === false

    @test <(interval(1.0,2.0), interval(3.0,4.0)) === true

    @test <(interval(1.0,3.5), interval(3.0,4.0)) === true

    @test <(interval(1.0,4.0), interval(3.0,4.0)) === false

    @test <(interval(0.0,4.0), interval(0.0,4.0)) === false

    @test <(interval(-0.0,4.0), interval(0.0,4.0)) === false

    @test <(interval(-2.0,-1.0), interval(-2.0,-1.0)) === false

    @test <(interval(-3.0,-1.5), interval(-2.0,-1.0)) === true

end

@testset "minimal_strictly_less_dec_test" begin

    @test <(nai(), nai()) === false

    @test <(DecoratedInterval(emptyinterval(), trv), nai()) === false

    @test <(DecoratedInterval(interval(1.0,2.0), trv), nai()) === false

    @test <(nai(), DecoratedInterval(interval(1.0,2.0), def)) === false

    @test <(DecoratedInterval(emptyinterval(), trv), DecoratedInterval(emptyinterval(), trv)) === true

    @test <(DecoratedInterval(interval(1.0,2.0), trv), DecoratedInterval(emptyinterval(), trv)) === false

    @test <(DecoratedInterval(emptyinterval(), trv), DecoratedInterval(interval(1.0,2.0), def)) === false

    @test <(DecoratedInterval(interval(-Inf,+Inf), trv), DecoratedInterval(interval(-Inf,+Inf), trv)) === true

    @test <(DecoratedInterval(interval(1.0,2.0), trv), DecoratedInterval(interval(-Inf,+Inf), trv)) === false

    @test <(DecoratedInterval(interval(-Inf,+Inf), trv), DecoratedInterval(interval(1.0,2.0), trv)) === false

    @test <(DecoratedInterval(interval(1.0,2.0), trv), DecoratedInterval(interval(1.0,2.0), trv)) === false

    @test <(DecoratedInterval(interval(1.0,2.0), trv), DecoratedInterval(interval(3.0,4.0), trv)) === true

    @test <(DecoratedInterval(interval(1.0,3.5), def), DecoratedInterval(interval(3.0,4.0), trv)) === true

    @test <(DecoratedInterval(interval(1.0,4.0), trv), DecoratedInterval(interval(3.0,4.0), def)) === false

    @test <(DecoratedInterval(interval(0.0,4.0), trv), DecoratedInterval(interval(0.0,4.0), def)) === false

    @test <(DecoratedInterval(interval(-0.0,4.0), def), DecoratedInterval(interval(0.0,4.0), trv)) === false

    @test <(DecoratedInterval(interval(-2.0,-1.0), def), DecoratedInterval(interval(-2.0,-1.0), def)) === false

    @test <(DecoratedInterval(interval(-3.0,-1.5), trv), DecoratedInterval(interval(-2.0,-1.0), trv)) === true

end

@testset "minimal_strictly_precedes_test" begin

    @test strictprecedes(emptyinterval(), interval(3.0,4.0)) === true

    @test strictprecedes(interval(3.0,4.0), emptyinterval()) === true

    @test strictprecedes(emptyinterval(), emptyinterval()) === true

    @test strictprecedes(interval(1.0,2.0), interval(-Inf,+Inf)) === false

    @test strictprecedes(interval(-Inf,+Inf), interval(1.0,2.0)) === false

    @test strictprecedes(interval(-Inf,+Inf), interval(-Inf,+Inf)) === false

    @test strictprecedes(interval(1.0,2.0), interval(3.0,4.0)) === true

    @test strictprecedes(interval(1.0,3.0), interval(3.0,4.0)) === false

    @test strictprecedes(interval(-3.0,-1.0), interval(-1.0,0.0)) === false

    @test strictprecedes(interval(-3.0,-0.0), interval(0.0,1.0)) === false

    @test strictprecedes(interval(-3.0,0.0), interval(-0.0,1.0)) === false

    @test strictprecedes(interval(1.0,3.5), interval(3.0,4.0)) === false

    @test strictprecedes(interval(1.0,4.0), interval(3.0,4.0)) === false

    @test strictprecedes(interval(-3.0,-0.1), interval(-1.0,0.0)) === false

end

@testset "minimal_strictly_precedes_dec_test" begin

    @test strictprecedes(nai(), DecoratedInterval(interval(3.0,4.0), trv)) === false

    @test strictprecedes(DecoratedInterval(interval(3.0,4.0), def), nai()) === false

    @test strictprecedes(nai(), DecoratedInterval(emptyinterval(), trv)) === false

    @test strictprecedes(nai(), nai()) === false

    @test strictprecedes(DecoratedInterval(emptyinterval(), trv), DecoratedInterval(interval(3.0,4.0), trv)) === true

    @test strictprecedes(DecoratedInterval(interval(3.0,4.0), def), DecoratedInterval(emptyinterval(), trv)) === true

    @test strictprecedes(DecoratedInterval(emptyinterval(), trv), DecoratedInterval(emptyinterval(), trv)) === true

    @test strictprecedes(DecoratedInterval(interval(1.0,2.0), trv), DecoratedInterval(interval(-Inf,+Inf), trv)) === false

    @test strictprecedes(DecoratedInterval(interval(-Inf,+Inf), trv), DecoratedInterval(interval(1.0,2.0), trv)) === false

    @test strictprecedes(DecoratedInterval(interval(-Inf,+Inf), trv), DecoratedInterval(interval(-Inf,+Inf), trv)) === false

    @test strictprecedes(DecoratedInterval(interval(1.0,2.0), trv), DecoratedInterval(interval(3.0,4.0), trv)) === true

    @test strictprecedes(DecoratedInterval(interval(1.0,3.0), def), DecoratedInterval(interval(3.0,4.0), trv)) === false

    @test strictprecedes(DecoratedInterval(interval(-3.0,-1.0), trv), DecoratedInterval(interval(-1.0,0.0), def)) === false

    @test strictprecedes(DecoratedInterval(interval(-3.0,-0.0), def), DecoratedInterval(interval(0.0,1.0), trv)) === false

    @test strictprecedes(DecoratedInterval(interval(-3.0,0.0), trv), DecoratedInterval(interval(-0.0,1.0), trv)) === false

    @test strictprecedes(DecoratedInterval(interval(1.0,3.5), trv), DecoratedInterval(interval(3.0,4.0), trv)) === false

    @test strictprecedes(DecoratedInterval(interval(1.0,4.0), trv), DecoratedInterval(interval(3.0,4.0), def)) === false

    @test strictprecedes(DecoratedInterval(interval(-3.0,-0.1), trv), DecoratedInterval(interval(-1.0,0.0), trv)) === false

end

@testset "minimal_disjoint_test" begin

    @test isdisjoint(emptyinterval(), interval(3.0,4.0)) === true

    @test isdisjoint(interval(3.0,4.0), emptyinterval()) === true

    @test isdisjoint(emptyinterval(), emptyinterval()) === true

    @test isdisjoint(interval(3.0,4.0), interval(1.0,2.0)) === true

    @test isdisjoint(interval(0.0,0.0), interval(-0.0,-0.0)) === false

    @test isdisjoint(interval(0.0,-0.0), interval(-0.0,0.0)) === false

    @test isdisjoint(interval(3.0,4.0), interval(1.0,7.0)) === false

    @test isdisjoint(interval(3.0,4.0), interval(-Inf,+Inf)) === false

    @test isdisjoint(interval(-Inf,+Inf), interval(1.0,7.0)) === false

    @test isdisjoint(interval(-Inf,+Inf), interval(-Inf,+Inf)) === false

end

@testset "minimal_disjoint_dec_test" begin

    @test isdisjoint(nai(), DecoratedInterval(interval(3.0,4.0), def)) === false

    @test isdisjoint(DecoratedInterval(interval(3.0,4.0), trv), nai()) === false

    @test isdisjoint(DecoratedInterval(emptyinterval(), trv), nai()) === false

    @test isdisjoint(nai(), nai()) === false

    @test isdisjoint(DecoratedInterval(emptyinterval(), trv), DecoratedInterval(interval(3.0,4.0), def)) === true

    @test isdisjoint(DecoratedInterval(interval(3.0,4.0), trv), DecoratedInterval(emptyinterval(), trv)) === true

    @test isdisjoint(DecoratedInterval(emptyinterval(), trv), DecoratedInterval(emptyinterval(), trv)) === true

    @test isdisjoint(DecoratedInterval(interval(3.0,4.0), trv), DecoratedInterval(interval(1.0,2.0), def)) === true

    @test isdisjoint(DecoratedInterval(interval(0.0,0.0), trv), DecoratedInterval(interval(-0.0,-0.0), trv)) === false

    @test isdisjoint(DecoratedInterval(interval(0.0,-0.0), trv), DecoratedInterval(interval(-0.0,0.0), trv)) === false

    @test isdisjoint(DecoratedInterval(interval(3.0,4.0), def), DecoratedInterval(interval(1.0,7.0), def)) === false

    @test isdisjoint(DecoratedInterval(interval(3.0,4.0), trv), DecoratedInterval(interval(-Inf,+Inf), trv)) === false

    @test isdisjoint(DecoratedInterval(interval(-Inf,+Inf), trv), DecoratedInterval(interval(1.0,7.0), trv)) === false

    @test isdisjoint(DecoratedInterval(interval(-Inf,+Inf), trv), DecoratedInterval(interval(-Inf,+Inf), trv)) === false

end


text:
include("Args.jl") # Parsing the command-line
##################################################
# General parameters
##################################################
"""
    PARALLEL

Determining if the code is run in parallel.
"""
const PARALLEL = parsed_args["parallel"]
if ((PARALLEL != "yes") && (PARALLEL != "no"))
    error("ERROR: UNKNOWN PARALLEL") # Unknown parallel procedure
end

"""
    nbK_default

Default number of sampling points for the 3D integrals.
"""
const nbK_default = parsed_args["nbK"]

"""
    nbAvr_default

Default number of sampling points for the orbit-averaging integrals.
"""
const nbAvr_default = parsed_args["nbAvr"]

"""
    qCalc

Anisotropy parameter q for the Plummer model.
"""
const qCalc = parsed_args["q"]

##################################################
# Options to be removed in GitHub
##################################################

"""
    gAni_default

Anisotropy parameter gAni for the Isochrone model.
"""
const gAni_default = parsed_args["gAni"]

"""
    logJrcutoff

log10(Jr) lower cut-off used in the computation of dRc/dt
"""
const logJrcutoff = parsed_args["logJrcutoff"]

"""
    logLcutoff

log10(L) lower cut-off used in the computation of dRc/dt
"""
const logLcutoff = parsed_args["logLcutoff"]

##################################################
##################################################

include("Constants.jl")
include("Mean.jl")
include("Bath.jl")
include("LocalDeflection.jl")
include("OrbitAverage.jl")
include("ActionSpace.jl")
include("Flux.jl")

max_stars_repo_name: KerwannTEP/CAT

max_issues_repo_name: KerwannTEP/CAT

max_forks_repo_name: KerwannTEP/CAT


text:
# using RoMEPlotting
using Caesar, Caesar.ZmqCaesar
using IncrementalInference, RoME
using Test

@testset "ZMQ Interface" begin
    # Unit tests
    include("multilangzmq/callbackCompatibilityTest.jl")
    # Integration tests
    include("multilangzmq/zmqInternal.jl")
    # End to end tests
    # Do this manually...
    # include("multilangzmq/serverTest.jl")
end

# println("[TEST] CloudGraphs API calls...")
# if false
#   println("[TEST] with CloudGraphs with local DB data layer (multicore)...")
#   include("fourdoortestcloudgraph.jl")
#   println("[SUCCESS]")
# else
#   @warn "[NOT TESTING] CloudGraphs interface, auth required -- you can enable it here at $(@__FILE__)."
# end

# if false
#   # using VictoriaParkTypes
#
# else
#   @warn "not running full Victoria park dataset now"
# end


text:
module VarianceComponentModels

using LinearAlgebra, Statistics, MathProgBase, Ipopt#, KNITRO#, Mosek#, Gurobi
import Base: eltype, length, size, +
export VarianceComponentModel, VarianceComponentVariate,
  TwoVarCompModelRotate, TwoVarCompVariateRotate, VarianceComponentAuxData,
  residual, residual!, nvarcomps, nmeanparams, nvarparams, nparams,
  mean!, mean, cov!, cov

"""
`VarianceComponentModel` stores the model parameters of a variance component
model.

# Fields
- `B`: `p x d` mean parameters
- `Σ`: tuple of `d x d` variance component parameters
- `A`: constraint matrix for `vec(B)`
- `sense`: vector of characters `'='`, `'<'` or `'>'`
- `b`: constraint vector for `vec(B)`
- `lb`: lower bounds for `vec(B)`
- `ub`: upper bounds for `vec(B)`
"""
mutable struct VarianceComponentModel{T <: AbstractFloat, M,
  BT <: AbstractVecOrMat, ΣT <: AbstractMatrix}
  # model parameters
  B::BT
  Σ::NTuple{M, ΣT}
  # constraints: A * vec(B) sense b, lb ≤ vec(B) ≤ ub
  A::Matrix{T}
  sense::Union{Char, Vector{Char}}
  b::Union{T, Vector{T}}
  lb::Union{T, Vector{T}}
  ub::Union{T, Vector{T}}
end

"""
    VarianceComponentModel(B, Σ)
    VarianceComponentModel(B, Σ, A, sense, b, lb, ub)

Default constructor of [`VarianceComponentModel`](@ref) type.
"""
function VarianceComponentModel(
  B::AbstractVecOrMat,
  Σ::NTuple{M, AbstractMatrix},
  A::Matrix = zeros(eltype(B), 0, length(B)),
  sense::Union{Char, Vector{Char}} = Array{Char}(undef, 0),
  b::Union{Number, Vector} = zeros(eltype(B), 0),
  lb::Union{Number, Vector} = convert(eltype(B), -Inf),
  ub::Union{Number, Vector} = convert(eltype(B), Inf),
  ) where {M}

  VarianceComponentModel{eltype(B), M, typeof(B), eltype(Σ)}(B, Σ, A, sense, b,
    lb, ub)
end

"""
    VarianceComponentModel(Σ)

Construct [`VarianceComponentModel`](@ref) from `Σ` alone. `B` is treated empty.
"""
function VarianceComponentModel(Σ::NTuple{M, AbstractMatrix}) where {M}
  B = zeros(eltype(Σ[1]), 0, size(Σ[1], 1))
  VarianceComponentModel(B, Σ)
end

"""
`TwoVarCompModelRotate` stores the rotated two variance component model.

# Fields
- `Brot`: rotated mean parameters `B * eigvec`
- `eigval`: eigenvalues of `eig(Σ[1], Σ[2])`
- `eigvec`: eigenvectors of `eig(Σ[1], Σ[2])`
- `logdetΣ2`: log-determinant of `Σ[2]`
"""
struct TwoVarCompModelRotate{T <: AbstractFloat, BT <: AbstractVecOrMat}
  Brot::BT
  eigval::Vector{T}
  eigvec::Matrix{T}
  logdetΣ2::T
end

"""
    TwoVarCompModelRotate(Brot, eigval, eigvec, logdetΣ2)

Default constructor of [`TwoVarCompModelRotate`](@ref) type.
"""
function TwoVarCompModelRotate(
  Brot::AbstractVecOrMat,
  eigval::Vector,
  eigvec::Matrix,
  logdetΣ2
  )

  TwoVarCompModelRotate{eltype(Brot), typeof(Brot)}(Brot, eigval, eigvec, logdetΣ2)
end

"""
    TwoVarCompModelRotate(twovarcomp)

Constructor of [`TwoVarCompModelRotate`](@ref) instance from a
[`TwoVarCompModelRotate`](@ref) instance itself.
"""
function TwoVarCompModelRotate(
  twovarcomp::TwoVarCompModelRotate{T}) where T <: AbstractFloat

  ## JK ###
  twovarcomp
end

"""
    TwoVarCompModelRotate(vcmodel)

Constructor of a [`TwoVarCompModelRotate`](@ref) instance from a
[`VarianceComponentModel`](@ref) instance.
"""
function TwoVarCompModelRotate(vcm::VarianceComponentModel{T, 2, BT, ΣT}) where {T, BT, ΣT}

  # generalized eigenvalue decomposition of (Σ1, Σ2)
  F = eigen(Symmetric(vcm.Σ[1]), Symmetric(vcm.Σ[2])) # output is of BlasFloat type
  λ = convert(Vector{T}, F.values)
  Φ = convert(Matrix{T}, F.vectors)
  # correct negative eigenvalues due to roundoff
  map!(x -> max(x, zero(T)), λ, λ)
  Brot = isempty(vcm.B) ? Array{T}(undef, size(vcm.B)) : vcm.B * Φ
  logdetΣ2 = convert(T, logdet(vcm.Σ[2]))
  TwoVarCompModelRotate{T, BT}(Brot, λ, Φ, logdetΣ2)
  #TwoVarCompModelRotate(Brot, λ, Φ, logdetΣ2) where {T, BT}
end

"""
`VarianceComponentVariate` stores the data of a variance component
model.

# Feilds
- `Y`: `n x d` responses
- `X`: `n x p` predictors
- `V`: tuple of `n x n` covariance matrices
"""
mutable struct VarianceComponentVariate{T <: AbstractFloat, M,
  YT <: AbstractVecOrMat, XT <: AbstractVecOrMat, VT <: AbstractMatrix}
  # data
  Y::YT
  X::XT
  V::NTuple{M, VT}
end

"""
    VarianceComponentVariate(Y, X, V)

Default constructor of [`VarianceComponentVariate`](@ref) type.
"""
function VarianceComponentVariate(
  Y::AbstractVecOrMat,
  X::AbstractVecOrMat,
  V::NTuple{M, AbstractMatrix}
  ) where {M}

  VarianceComponentVariate{eltype(Y), M, typeof(Y), typeof(X), eltype(V)}(Y, X, V)
end

"""
   VarianceComponentVariate(Y, V)

Constructor of a [`VarianceComponentVariate`](@ref) instance from `Y` and `V`
alone. `X` is created empty.
"""
function VarianceComponentVariate(
  Y::AbstractVecOrMat,
  V::NTuple{M, AbstractMatrix}
  ) where {M}

  X = zeros(eltype(Y), size(Y, 1), 0)
  VarianceComponentVariate{eltype(Y), M, typeof(Y), typeof(X), eltype(V)}(Y, X, V)
end

"""
`TwoVarCompVariateRotate` stores the rotated two variance component data.

# Fields
- `Yrot`: rotated responses `eigvec * Y`
- `Xrot`: rotated covariates `eigvec * X`
- `eigval`: eigenvalues of `eig(V[1], V[2])`
- `eigvec`: eigenvectors of `eig(V[1], V[2])`
- `logdetV2`: log-determinant of `V[2]`
"""
struct TwoVarCompVariateRotate{T <: AbstractFloat, YT <: AbstractVecOrMat,
  XT <: AbstractVecOrMat}
  # data
  Yrot::YT
  Xrot::XT
  eigval::Vector{T}
  eigvec::Matrix{T}
  logdetV2::T
end

"""
    TwoVarCompVariateRotate(Yrot, Xrot, eigval, eigvec,logdetV2)

Default constructor of a [`TwoVarCompVariateRotate`](@ref) instance.
"""
function TwoVarCompVariateRotate(
  Yrot::AbstractVecOrMat,
  Xrot::AbstractVecOrMat,
  eigval::Vector,
  eigvec::Matrix,
  logdetV2::Real)

  ## JZ ###
  TwoVarCompVariateRotate{eltype(Yrot), typeof(Yrot), typeof(Xrot)}(Yrot, Xrot,
    eigval, eigvec, logdetV2)
end

"""
    TwoVarCompVariateRotate(twovarcomp)

Constructor of a [`TwoVarCompVariateRotate`](@ref) instance from a
[`TwoVarCompVariateRotate`](@ref) instance itself.
"""
function TwoVarCompVariateRotate(
  twovarcomp::TwoVarCompVariateRotate{T}) where T <: AbstractFloat

  ## JK ###
  twovarcomp
end

"""
    TwoVarCompVariateRotate(vcobs)

Constructor of a [`TwoVarCompVariateRotate`](@ref) instance from a
[`VarianceComponentVariate`](@ref) instance.
"""
function TwoVarCompVariateRotate(vcobs::VarianceComponentVariate{T, 2}) where T <: AbstractFloat

  zeroT = zero(T)
  # (generalized)-eigendecomposition of (V1, V2)
  if isa(vcobs.V[2], UniformScaling) ||
    (isdiag(vcobs.V[2]) && norm(diag(vcobs.V[2]) .- one(T)) < 1.0e-8)
    F = eigen(Symmetric(vcobs.V[1]))
    deval = convert(Vector{T}, F.values)
    U = convert(Matrix{T}, F.vectors)
    logdetV2 = zero(T)
  else
    F = eigen(Symmetric(vcobs.V[1]), Symmetric(vcobs.V[2]))
    deval = convert(Vector{T}, F.values)
    U = convert(Matrix{T}, F.vectors)
    logdetV2 = convert(T, logdet(vcobs.V[2]))
  end
  # corect negative eigenvalues due to roundoff error
  @inbounds @simd for i in eachindex(deval)
    deval[i] = deval[i] > zeroT ? deval[i] : zeroT
  end
  # rotate responses
  Yrot = transpose(U) * vcobs.Y
  Xrot = isempty(vcobs.X) ? Array{T}(undef, size(Yrot, 1), 0) : transpose(U) * vcobs.X
  # output
  TwoVarCompVariateRotate(Yrot, Xrot, deval, U, logdetV2)
end

function TwoVarCompVariateRotate(vcdata::Array{T}) where T <: VarianceComponentVariate

  map(x -> TwoVarCompVariateRotate(x), vcdata)
end

"""
    VarianceComponentModel(vcobs)

Construct a [`VarianceComponentModel`](@ref) instance from a [`VarianceComponentVariate`](@ref)
instance. `B` is initialized to zero; `Σ` is initialized to a tupe of identity
matrices.
"""
function VarianceComponentModel(vcobs::VarianceComponentVariate{T, M}) where {T, M}

  p, d, m = size(vcobs.X, 2), size(vcobs.Y, 2), length(vcobs.V)
  B = zeros(T, p, d)
  Σ = ntuple(x -> Matrix{T}(I, d, d), m)::NTuple{M, Matrix{T}}
  VarianceComponentModel(B, Σ)
end

"""
    VarianceComponentModel(vcobsrot)

Construct [`VarianceComponentModel`](@ref) instance from a [`TwoVarCompVariateRotate`](@ref)
instance.
"""
function VarianceComponentModel(vcobsrot::TwoVarCompVariateRotate)

  p, d = size(vcobsrot.Xrot, 2), size(vcobsrot.Yrot, 2)
  T = eltype(vcobsrot)
  B = zeros(T, p, d)
  Σ = (Matrix{T}(I, d, d), Matrix{T}(I, d, d))
  VarianceComponentModel(B, Σ)
end

"""
Auxillary data for variance component model. This can be used as pre-allocated
intermediate variable in iterative algorithm.
"""
struct VarianceComponentAuxData{
  T1 <: AbstractMatrix,
  T2 <: AbstractVector
  }

  res::T1     # same shape as response, p-by-d
  Xwork::T1   # nd-by-pd
  ywork::T2   # nd
  obswt::T2   # nd
end

function VarianceComponentAuxData(vcobs::VarianceComponentVariate)
  T     = eltype(vcobs)
  res   = zeros(T, size(vcobs.Y, 1), size(vcobs.Y, 2))
  Xwork = zeros(T, length(vcobs.Y), nmeanparams(vcobs))
  ywork = zeros(T, length(vcobs.Y))
  obswt = zeros(T, length(vcobs.Y))
  VarianceComponentAuxData{typeof(Xwork), typeof(ywork)}(res,
    Xwork, ywork, obswt)
end

function VarianceComponentAuxData(vcobsrot::TwoVarCompVariateRotate)
  T     = eltype(vcobsrot)
  res   = zeros(T, size(vcobsrot.Yrot, 1), size(vcobsrot.Yrot, 2))
  Xwork = zeros(T, length(vcobsrot.Yrot), nmeanparams(vcobsrot))
  ywork = zeros(T, length(vcobsrot.Yrot))
  obswt = zeros(T, length(vcobsrot.Yrot))
  VarianceComponentAuxData{typeof(Xwork), typeof(ywork)}(res,
    Xwork, ywork, obswt)
end

function VarianceComponentAuxData(vcdata::Array{T}) where
  T <: Union{VarianceComponentVariate, TwoVarCompVariateRotate}

  map(x -> VarianceComponentAuxData(x), vcdata)
end

"""
    eltype(vcm)
    eltype(vcobs)
    eltype(vcmrot)
    eltype(vcobsrot)

Element type of variance component model or data.
"""
Base.eltype(vcm::VarianceComponentModel) = Base.eltype(vcm.B)
Base.eltype(vcobs::VarianceComponentVariate) = Base.eltype(vcobs.Y)
Base.eltype(vcmrot::TwoVarCompModelRotate) = Base.eltype(vcmrot.Brot)
Base.eltype(vcobsrot::TwoVarCompVariateRotate) = Base.eltype(vcobsrot.Yrot)
# Dimension of response, d
"""
    length(vcm)
    length(vcobs)
    length(vcmrot)
    length(vcobsrot)

Length `d` of response.
"""
length(vcm::VarianceComponentModel) = size(vcm.Σ[1], 1)
length(vcobs::VarianceComponentVariate) = size(vcobs.Y, 2)
length(vcmrot::TwoVarCompModelRotate) = length(vcmrot.eigval)
length(vcobsrot::TwoVarCompVariateRotate) = size(vcobsrot.Yrot, 2)
# Size of response, (n, d)
"""
    size(vcobs)
    size(vcobsrot)

Size `(n, d)` of the response matrix of a [`VarianceComponentVariate`](@ref) or
[`TwoVarCompVariateRotate`](@ref).
"""
size(vcobs::VarianceComponentVariate) = size(vcobs.Y)
size(vcobsrot::TwoVarCompVariateRotate) = size(vcobsrot.Yrot)
# Number of variance components, m
"""
    nvarcomps(vcm)
    nvarcomps(vcobs)
    nvarcomps(vcmrot)
    nvarcomps(vcobsrot)

Number of varianece components.
"""
nvarcomps(vcm::VarianceComponentModel) = length(vcm.Σ)
nvarcomps(vcobs::VarianceComponentVariate) = length(vcobs.V)
nvarcomps(vcmrot::TwoVarCompModelRotate) = 2
nvarcomps(vcobsrot::TwoVarCompVariateRotate) = 2
# Number of mean parameters, p * d
"""
    nmeanparams(vcm)
    nmeanparams(vcmrot)
    nmeanparams(vcobs)
    nmeanparams(vcobsrot)

Number of mean parameters `p * d` of [`VarianceComponentModel`](@ref).
"""
nmeanparams(vcm::VarianceComponentModel) = length(vcm.B)
nmeanparams(vcmrot::TwoVarCompModelRotate) = length(vcmrot.Brot)
nmeanparams(vcobs::VarianceComponentVariate) = size(vcobs.X, 2) * size(vcobs.Y, 2)
nmeanparams(vcobsrot::TwoVarCompVariateRotate) =
  size(vcobsrot.Xrot, 2) * size(vcobsrot.Yrot, 2)
# Number of free parameters in Cholesky factors, m * d * (d + 1) / 2
"""
    nvarparams(vcm)

Number of free parameters in Cholesky factors `m * d * (d + 1) / 2` of
[`VarianceComponentModel`](@ref).
"""
nvarparams(vcm::Union{VarianceComponentModel, TwoVarCompModelRotate}) =
  nvarcomps(vcm) * binomial(length(vcm) + 1, 2)
# Number of model parameters, p * d + m * d * (d + 1) / 2
"""
    nparams(vcm)

Total number of model parameters `p * d + m * d * (d + 1) / 2` of
[`VarianceComponentModel`](@ref).
"""
nparams(vcm::Union{VarianceComponentModel, TwoVarCompModelRotate}) =
  nmeanparams(vcm) + nvarparams(vcm)

"""
    cov!(C, vcm, vcobs)

Calculate the `nd x nd` covariance matrix of a [`VarianceComponentModel`](@ref) at
a [`VarianceComponentVariate`](@ref) observation and over-write `C`.
"""
function cov!(
  C::AbstractMatrix,
  vcm::VarianceComponentModel,
  vcobs::VarianceComponentVariate
  )

  fill!(C, zero(eltype(vcm)))
  for i in 1:length(vcm.Σ)
    kronaxpy!(vcm.Σ[i], vcobs.V[i], C)
  end
  return C
end

"""
    cov(vcm, vcobs)

Calculate the `nd x nd` covariance matrix of a [`VarianceComponentModel`](@ref) at
a [`VarianceComponentVariate`](@ref) observation.
"""
function cov(
  vcm::VarianceComponentModel,
  vcobs::VarianceComponentVariate
  )

  n, d = size(vcobs)
  C = zeros(eltype(vcm), n * d, n * d)
  cov!(C, vcm, vcobs)
end

"""
    mean!(μ, vcm, vcobs)

Calculate the `n x d` mean matrix of a a [`VarianceComponentModel`](@ref) at
a [`VarianceComponentVariate`](@ref) observation and over-write `μ`.
"""
function mean!(
  μ::AbstractMatrix,
  vcm::VarianceComponentModel,
  vcobs::VarianceComponentVariate
  )

  mul!(μ, vcobs.X, vcm.B)
end

"""
    mean(vcm, vcobs)

Calculate the `n x d` mean matrix of a [`VarianceComponentModel`](@ref) at
a [`VarianceComponentVariate`](@ref) observation.
"""
function mean(vcm::VarianceComponentModel, vcobs::VarianceComponentVariate)
  μ = zeros(eltype(vcobs), size(vcobs.Y))
  mean!(μ, vcm, vcobs)
end

"""
    residual(vcm, vcobs)

Residual of [`VarianceComponentVariate`](@ref) under [`VarianceComponentModel`](@ref).
"""
function residual(vcm::VarianceComponentModel, vcobs::VarianceComponentVariate)
  return isempty(vcobs.X) ? vcobs.Y : vcobs.Y - vcobs.X * vcm.B
end

function residual(
  vcm::T1,
  vcdata::Array{T2}
  ) where {T1 <: VarianceComponentModel, T2 <: VarianceComponentVariate}

  map(x -> residual(vcm, x), vcdata)
end

function residual!(
  resid::AbstractVecOrMat,
  vcm::VarianceComponentModel,
  vcobs::VarianceComponentVariate
  )

  if isempty(vcobs.X)
    copyto!(resid, vcobs.Y)
  else
    oneT = one(eltype(vcm))
    copyto!(resid, vcobs.Y)
    LinearAlgebra.BLAS.gemm!('N', 'N', -oneT, vcobs.X, vcm.B, oneT, resid)
  end
  resid
end

function residual!(
  resid::AbstractArray,
  vcm::T1,
  vcdata::Array{T2}
  ) where {T1 <: VarianceComponentModel, T2 <: VarianceComponentVariate}

  map(x -> residual!(x[1], vcm, x[2]), zip(resid, vcdata))
end

"""
    residual(vcmrot, vcobsrot)

Residual of [`TwoVarCompVariateRotate`](@ref) under [`TwoVarCompModelRotate`](@ref).
"""
function residual(vcmrot::TwoVarCompModelRotate, vcobsrot::TwoVarCompVariateRotate)
  if isempty(vcmrot.Brot)
    return vcobsrot.Yrot * vcmrot.eigvec
  else
    return vcobsrot.Yrot * vcmrot.eigvec - vcobsrot.Xrot * vcmrot.Brot
  end
end

function residual(
  vcmrot::T1,
  vcdatarot::Array{T2}
  ) where {T1 <: TwoVarCompModelRotate, T2 <: TwoVarCompVariateRotate}

  map(x -> residual(vcmrot, x), vcdatarot)
end

function residual!(
  resid::AbstractMatrix,
  vcmrot::TwoVarCompModelRotate,
  vcobsrot::TwoVarCompVariateRotate
  )

  if isempty(vcmrot.Brot)
    mul!(resid, vcobsrot.Yrot, vcmrot.eigvec)
  else
    #copy!(resid, vcobsrot.Yrot * vcmrot.eigvec - vcobsrot.Xrot * vcmrot.Brot)
    oneT = one(eltype(vcmrot))
    mul!(resid, vcobsrot.Yrot, vcmrot.eigvec)
    LinearAlgebra.BLAS.gemm!('N', 'N', -oneT, vcobsrot.Xrot, vcmrot.Brot, oneT, resid)
  end
  resid
end

function residual!(
  resid::Array{AbstractMatrix},
  vcmrot::T1,
  vcdatarot::Array{T2}) where {T1 <: TwoVarCompModelRotate, T2 <: TwoVarCompVariateRotate}

  map!(x -> residual(vcmrot, x), resid, vcdatarot)
end

# utilities for multivariate calculus
include("multivariate_calculus.jl")
# source for fitting models with 2 variance components
include("two_variance_component.jl")

end # VarianceComponentModels

max_stars_repo_name: Hua-Zhou/VarianceComponentModels

max_issues_repo_name: Hua-Zhou/VarianceComponentModels

max_forks_repo_name: Hua-Zhou/VarianceComponentModels


text:
grid = ones(Float64, 4, 4)
grid[1:2, :] = 4.0
grid[3:4, :] = 0.2
cut = MinCut(size(grid)...)
reset!(cut, grid, 1.0, 1.0)
out = ones(UInt8, 4, 4)
out[3:4, :] = zero(UInt8)
@test segment!(cut) == out


text:
# depends on scene.jl

##############################################
# rapidly exploring random tree for 2D point #
##############################################

mutable struct RRTTree
    confs::Matrix{Float64} # 2 x max_num_nodes
    control_xs::Vector{Float64}
    control_ys::Vector{Float64}
    costs_from_start::Vector{Float64}
    parents::Vector{UInt}
    num_nodes::UInt
    max_num_nodes::UInt
end

function RRTTree(root_conf::Point, iters::UInt)
    max_num_nodes = iters + 2

    # pre-allocate
    confs = Matrix{Float64}(undef, 2, max_num_nodes)
    control_xs = Vector{Float64}(undef, max_num_nodes)
    control_ys = Vector{Float64}(undef, max_num_nodes)
    costs_from_start = Vector{Float64}(undef, max_num_nodes)
    parents = Vector{UInt}(undef, max_num_nodes)

    # initialize
    confs[1,1] = root_conf.x
    confs[2,1] = root_conf.y
    control_xs[1] = NaN
    control_ys[1] = NaN
    costs_from_start[1] = 0.0
    parents[1] = 0

    return RRTTree(confs, control_xs, control_ys, costs_from_start, parents, 1, max_num_nodes)
end

@inline function add_node!(tree::RRTTree, parent::UInt, new_conf::Point, control::Point, cost_from_start::Float64)
    node = tree.num_nodes + 1
    tree.confs[1,node] = new_conf.x
    tree.confs[2,node] = new_conf.y
    tree.control_xs[node] = control.x # TODO fix make dense..
    tree.control_ys[node] = control.y
    tree.costs_from_start[node] = cost_from_start
    tree.parents[node] = parent
    tree.num_nodes = node
    return node
end

# TODO modify this so it stochastically sometimes does not pick the nearest neighbor?

function nearest_neighbor(conf::Point, confs::Matrix{Float64}, num_nodes::UInt)
    nearest::UInt = 1
    @inbounds dx::Float64 = confs[1,1] - conf.x
    @inbounds dy::Float64 = confs[2,1] - conf.y
    min_so_far::Float64 = dx * dx + dy * dy
    local d::Float64
    for node in 2:num_nodes
        @inbounds dx = confs[1,node] - conf.x
        @inbounds dy = confs[2,node] - conf.y
        d = dx * dx + dy * dy
        if d < min_so_far
            nearest = node
            min_so_far = d
        end
    end
    return nearest
end

struct SelectControlResult
    start_conf::Point
    new_conf::Point
    control::Point
    failed::Bool # new_conf is undefined in this case
    cost::Float64 # cost of this control action (e.g. distance)
end

@inline function random_config(bounds)
    x = rand()* (bounds.xmax - bounds.xmin) + bounds.xmin
    y = rand() * (bounds.ymax - bounds.ymin) + bounds.ymin
    return Point(x, y)
end

function select_control(walls, target_conf::Point, start_conf::Point, dt::Float64)

    dist_to_target = dist(start_conf, target_conf)
    diff = Point(target_conf.x - start_conf.x, target_conf.y - start_conf.y)
    distance_to_move = min(dt, dist_to_target)
    scale = distance_to_move / dist_to_target
    control = Point(scale * diff.x, scale * diff.y)

    # go in the direction of target_conf from start_conf
    new_conf = Point(start_conf.x + control.x, start_conf.y + control.y)

    # test the obstacles
    failed = false
    for wall in walls
        if line_intersects_line(wall.a, wall.b, start_conf, new_conf)
            failed = true
            break
        end
    end
    cost = distance_to_move
    return SelectControlResult(start_conf, new_conf, control, failed, cost)
end

function rrt(bounds::Bounds, walls::Vector, init::Point, iters::UInt, dt::Float64)
    tree = RRTTree(init, iters)
    local near_node::UInt
    for iter=1:iters
        rand_conf::Point = random_config(bounds)
        near_node = nearest_neighbor(rand_conf, tree.confs, tree.num_nodes)
        result = select_control(walls, rand_conf, Point(tree.confs[1,near_node], tree.confs[2,near_node]), dt)
        if !result.failed
            cost_from_start = tree.costs_from_start[near_node] + result.cost
            add_node!(tree, near_node, result.new_conf, result.control, cost_from_start)
        end
    end
    return tree
end

##############################
# path planner that uses RRT #
##############################

struct PlannerParams
    rrt_iters::UInt
    rrt_dt::Float64 # the maximum proposal distance
    refine_spacing::Float64
    refine_iters::UInt
    refine_std::Float64
    refine_gap::Float64
end

struct Path
    start::Point
    goal::Point
    points::Array{Point,1}
end

function concatenate(a::Path, b::Path)
    if a.goal.x != b.start.x || a.goal.y != b.start.y
        error("goal of first path muts be start of second path")
    end
    points = Array{Point,1}()
    for point in a.points
        push!(points, point)
    end
    for point in b.points[2:end]
        push!(points, point)
    end
    @assert points[1].x == a.start.x
    @assert points[1].y == a.start.y
    @assert points[end].x == b.goal.x
    @assert points[end].y == b.goal.y
    Path(a.start, b.goal, points)
end

function add_intermediate_points!(new_points, next_point, spacing)
    # add intermediate points, evently spaced
    d = dist(next_point, new_points[end])
    while d > spacing # add points
        intermediate_point = ((d-spacing)/d) * new_points[end] + (spacing/d) * next_point
        @assert isapprox(dist(intermediate_point, new_points[end]), spacing)
        push!(new_points, intermediate_point)
        d = dist(next_point, new_points[end])
    end
end

"""
Remove intermediate nodes that are not necessary, so that refinement
optimization is a lower-dimensional optimization problem.
"""
function simplify_path(scene::Scene, original::Path; spacing=0.05)
    new_points = Array{Point,1}()
    push!(new_points, original.start)
    for i=2:length(original.points) - 1

        # try skipping this point
        if !line_of_site(scene, new_points[end], original.points[i + 1])
            add_intermediate_points!(new_points, original.points[i], spacing)
            push!(new_points, original.points[i])
        end

        # then we can skip it?
    end
    @assert line_of_site(scene, new_points[end], original.goal)
    add_intermediate_points!(new_points, original.goal, spacing)
    push!(new_points, original.goal)
    path = Path(original.start, original.goal, new_points)
    #@assert check_path_clear(path, scene)
    return path
end

function local_score(scene, a, b, num_tests, perturb_width)
    score = 0.
    for i in 1:num_tests
        perturbed_a = Point(a.x + uniform(-perturb_width, perturb_width), a.y + uniform(-perturb_width, perturb_width))
        perturbed_b = Point(b.x + uniform(-perturb_width, perturb_width), b.y + uniform(-perturb_width, perturb_width))
        score += line_of_site(scene, perturbed_a, perturbed_b)
    end
    return score / num_tests
end

"""
Optimize the path to minimize its length while avoiding obstacles in the scene.
"""
function refine_path(
        scene::Scene, original::Path, iters::UInt, std::Float64, gap::Float64) # 0.02

    # do stochastic optimization
    new_points = deepcopy(original.points)
    num_interior_points = length(original.points) -2
    if num_interior_points == 0
        return original
    end
    for i=1:iters
        point_idx = 2 + (i % num_interior_points)
        @assert point_idx > 1 # not start
        @assert point_idx < length(original.points) # not goal
        prev_point = new_points[point_idx-1]
        point = new_points[point_idx]
        next_point = new_points[point_idx+1]
        adjusted = Point(point.x + randn() * std, point.y + randn() * std)
        old_dist_sq = dist_sq(prev_point, point) + dist_sq(point, next_point)
        new_dist_sq = dist_sq(prev_point, adjusted) + dist_sq(adjusted, next_point)

        cur_1 = Point(point.x + gap, point.y)
        cur_2 = Point(point.x - gap, point.y)
        cur_3 = Point(point.x, point.y + gap)
        cur_4 = Point(point.x, point.y - gap)

        adjusted_1 = Point(adjusted.x + gap, adjusted.y)
        adjusted_2 = Point(adjusted.x - gap, adjusted.y)
        adjusted_3 = Point(adjusted.x, adjusted.y + gap)
        adjusted_4 = Point(adjusted.x, adjusted.y - gap)

        cur_cost = 0.001 * old_dist_sq
        adjusted_cost = 0.001 * new_dist_sq
        for wall in scene.walls

            cur_cost += (
                intersects(wall, prev_point, cur_1) +
                intersects(wall, cur_1, next_point) +
                intersects(wall, prev_point, cur_2) +
                intersects(wall, cur_2, next_point) +
                intersects(wall, prev_point, cur_3) +
                intersects(wall, cur_3, next_point) +
                intersects(wall, prev_point, cur_4) +
                intersects(wall, cur_4, next_point))
                #Inf * (
                    #intersects(wall, prev_point, point) +
                    #intersects(wall, point, next_point)
                #))

            adjusted_cost += (
                intersects(wall, prev_point, adjusted_1) +
                intersects(wall, adjusted_1, next_point) +
                intersects(wall, prev_point, adjusted_2) +
                intersects(wall, adjusted_2, next_point) +
                intersects(wall, prev_point, adjusted_3) +
                intersects(wall, adjusted_3, next_point) +
                intersects(wall, prev_point, adjusted_4) +
                intersects(wall, adjusted_4, next_point))
                #Inf * (
                    #intersects(wall, prev_point, adjusted) +
                    #intersects(wall, adjusted, next_point)
                #))
        end
        if adjusted_cost < cur_cost
            # buffer scores between 0 and 4
            new_points[point_idx] = adjusted
        end
    end
    Path(original.start, original.goal, new_points)
end

function optimize_path(
        scene::Scene, original::Path,
        refine_iters::UInt, refine_std::Float64, spacing::Float64, refine_gap::Float64)
    simplified = simplify_path(scene, original; spacing=spacing)
    #@assert check_path_clear(simplified, scene)
    refined = refine_path(scene, simplified, refine_iters, refine_std, refine_gap)
    #@assert check_path_clear(refined, scene)
    return refined
end

@inline function get_conf_point(tree, node::UInt)
    return Point(tree.confs[1,node], tree.confs[2,node])
end

function check_path_clear(path::Path, scene::Scene)
    prev = path.start
    @assert path.points[1] == path.start
    @assert path.points[end] == path.goal
    for i in 2:length(path.points)
        next = path.points[i]
        if !line_of_site(scene, prev, next)
            return false
        end
        prev = next
    end
    return true
end

"""
Plan path from start to goal that avoids obstacles in the scene.
"""
function plan_path(
        start::Point, goal::Point, scene::Scene,
        params::PlannerParams=PlannerParams(2000, 3.0, 0.05, 10000, 1.0, 0.02))

    tree = rrt(scene.bounds, scene.walls, start, params.rrt_iters, params.rrt_dt)

    # find the best path along the tree to the goal, if one exists
    best_node = 0
    min_cost = Inf
    path_found = false
    best_conf = Point(NaN,NaN)
    all_valid_nodes = UInt[]
    for node in 1:tree.num_nodes
        # check for line-of-site to the goal
        conf = get_conf_point(tree, node)
        clear_path = line_of_site(scene, conf, goal)
        cost = tree.costs_from_start[node] + (clear_path ? dist(conf, goal) : Inf)
        if clear_path
            push!(all_valid_nodes, node)
        end
        if cost < min_cost
            path_found = true
            best_node = node
            min_cost = cost
            best_conf = conf
        end
    end
    if path_found
        best_node = rand(all_valid_nodes)
        best_conf = get_conf_point(tree, best_node)
    end

    local path::Union{Nothing,Path}
    if path_found
        # extend the tree to the goal configuration
        control = Point(goal.x - best_conf.x, goal.y - best_conf.y)
        goal_node = add_node!(tree, best_node, goal, control, NaN)
        points = Array{Point,1}()
        node = goal_node
        push!(points, get_conf_point(tree, goal_node))
        # the path will contain the start and goal
        while node != 1
            node = tree.parents[node]
            push!(points, get_conf_point(tree, node))
        end
        @assert points[end] == start # the start point
        @assert points[1] == goal
        path = Path(start, goal, reverse(points))
        #@assert check_path_clear(path, scene)
    else
        path = nothing
    end

    return path, tree
end

function plan_and_optimize_path(scene, prev_loc, loc, params)
    (path, tree) = plan_path(prev_loc, loc, scene, params)
    if isnothing(path)
        return Point[], true, tree
    else
        path = optimize_path(scene, path, params.refine_iters, params.refine_std, params.refine_spacing, params.refine_gap)
        #@assert check_path_clear(path, scene)
        points = path.points
        @assert points[1] == prev_loc
        @assert points[end] == loc
        return points[2:end], false, tree
    end
end

export PlannerParams, plan_and_optimize_path


text:
## user values (<: llvm::User)

abstract type User <: Value end

# operand iteration

export operands

struct UserOperandSet
    user::User
end

operands(user::User) = UserOperandSet(user)

Base.eltype(::UserOperandSet) = Value

function Base.getindex(iter::UserOperandSet, i::Integer)   # TODO: otherwise unitrange indexing errors
    @boundscheck 1 <= i <= length(iter) || throw(BoundsError(iter, i))
    return Value(API.LLVMGetOperand(iter.user, i-1))
end

Base.setindex!(iter::UserOperandSet, val::Value, i) =
    API.LLVMSetOperand(iter.user, i-1, val)

function Base.iterate(iter::UserOperandSet, i=1)
    i >= length(iter) + 1 ? nothing : (iter[i], i+1)
end

Base.length(iter::UserOperandSet) = API.LLVMGetNumOperands(iter.user)

Base.lastindex(iter::UserOperandSet) = length(iter)