text:
# using Languages
using SimilaritySearch, TextSearch
using Test, SparseArrays, LinearAlgebra, StatsBase
const fit = TextSearch.fit

const text0 = "@user;) #jello.world"
const text1 = "hello world!! @user;) #jello.world :)"
const text2 = "a b c d e f g h i j k l m n o p q"
const corpus = ["hello world :)", "@user;) excellent!!", "#jello world."]

@testset "Character q-grams" begin
    config = TextConfig()
    config.group_usr = false
    config.nlist = []
    config.qlist = [3]
    config.slist = []
    a = [(p.first, p.second) for p in compute_bow(tokenize(config, text0))[1]]
    b = [(:rld, 1), (Symbol("@us"), 1), (:orl, 1), (:wor, 1), (Symbol("lo."), 1), (:use, 1), (:ser, 1), (Symbol("o.w"), 1), (:llo, 1), (Symbol("r;)"), 1), (:jel, 1), (Symbol(" #j"), 1), (Symbol("er;"), 1), (Symbol("ld "), 1), (Symbol(";) "), 1), (Symbol(" @u"), 1), (:ell, 1), (Symbol(") #"), 1), (Symbol(".wo"), 1), (Symbol("#je"), 1)]
    @test sort(a) == sort(b)
end

@testset "Word n-grams" begin
    config = TextConfig()
    config.group_usr = false
    config.nlist = [1, 2]
    config.qlist = []
    config.slist = []
    a = [(p.first, p.second) for p in compute_bow(tokenize(config, text0))[1]]
    b = [(Symbol("#jello"), 1), (Symbol("#jello ."), 1), (:., 1), (Symbol(". world"), 1), (Symbol(";)"), 1), (Symbol(";) #jello"), 1), (Symbol("@user"), 1), (Symbol("@user ;)"), 1), (:world, 1)]
    @test sort(a) == sort(b)
 end

@testset "Skip-grams" begin
    config = TextConfig()
    config.nlist = []
    config.qlist = []
    config.del_punc = true
    config.slist = [(2,1), (2, 2), (3, 1), (3, 2)]
    #L = collect(compute_bow(text2, config))
    #sort!(L)
    a = [(p.first, p.second) for p in compute_bow(tokenize(config, text2))[1]]
    b = [(Symbol("a c"), 1), (Symbol("a c e"), 1), (Symbol("a d"), 1), (Symbol("a d g"), 1), (Symbol("b d"), 1), (Symbol("b d f"), 1), (Symbol("b e"), 1), (Symbol("b e h"), 1), (Symbol("c e"), 1), (Symbol("c e g"), 1), (Symbol("c f"), 1), (Symbol("c f i"), 1), (Symbol("d f"), 1), (Symbol("d f h"), 1), (Symbol("d g"), 1), (Symbol("d g j"), 1), (Symbol("e g"), 1), (Symbol("e g i"), 1), (Symbol("e h"), 1), (Symbol("e h k"), 1), (Symbol("f h"), 1), (Symbol("f h j"), 1), (Symbol("f i"), 1), (Symbol("f i l"), 1), (Symbol("g i"), 1), (Symbol("g i k"), 1), (Symbol("g j"), 1), (Symbol("g j m"), 1), (Symbol("h j"), 1), (Symbol("h j l"), 1), (Symbol("h k"), 1), (Symbol("h k n"), 1), (Symbol("i k"), 1), (Symbol("i k m"), 1), (Symbol("i l"), 1), (Symbol("i l o"), 1), (Symbol("j l"), 1), (Symbol("j l n"), 1), (Symbol("j m"), 1), (Symbol("j m p"), 1), (Symbol("k m"), 1), (Symbol("k m o"), 1), (Symbol("k n"), 1), (Symbol("k n q"), 1), (Symbol("l n"), 1), (Symbol("l n p"), 1), (Symbol("l o"), 1), (Symbol("m o"), 1), (Symbol("m o q"), 1), (Symbol("m p"), 1), (Symbol("n p"), 1), (Symbol("n q"), 1), (Symbol("o q"), 1)]
    @test sort(a) == sort(b)
end

@testset "Tokenizer, DVEC, and vectorize" begin # test_vmodel
    config = TextConfig()
    config.nlist = [1]
    config.qlist = []
    config.slist = []
    config.group_usr = false

    @test tokenize(config, text1) == [Symbol(h) for h in ["hello", "world", "!!",  "@user", ";)", "#jello", ".", "world", ":)"]]
    model = fit(VectorModel, config, corpus)
    x = vectorize(model, TfModel, text1)
    @test nnz(x) == 8
    x = vectorize(model, TfModel, text2)
    @test nnz(x) == 0
end

const labeled_corpus = [("me gusta", 1), ("me encanta", 1), ("lo odio", 2), ("odio esto", 2), ("me encanta esto LOL!", 1)]
const sentiment_text = "lol, esto me encanta"

@testset "DistModel tests" begin
    config = TextConfig()
    config.nlist = [1]
    dmodel = fit(DistModel, config, [x[1] for x in labeled_corpus], [x[2] for x in labeled_corpus])
    @show sentiment_text
    @show dmodel
    #a = vectorize(dmodel, TfIdf, sentiment_text)
    #b = [(:me1,1.0),(:me2,0.0),(:encanta1,1.0),(:encanta2,0.0),(:esto1,0.4),(:esto2,0.6),(:lol1,1.0),(:lol2,0.0)]
    #@test string(a) == string(b)

end

@testset "EntModel tests" begin
    config = TextConfig()
    config.nlist = [1]
    X = [x[1] for x in labeled_corpus]
    y = [x[2] for x in labeled_corpus]
    dmodel = fit(DistModel, config, X, y, weights=:balance, smooth=1)
    emodel = fit(EntModel, dmodel)
    emodel_ = fit(EntModel, config, X, y, weights=:balance, smooth=1)
    a = vectorize(emodel, X)
    b = vectorize(emodel_, X)
    @test 0.999 < dot(a, b)
 end

@testset "distances" begin
    u = Dict(:el => 0.9, :hola => 0.1, :mundo => 0.2) |> normalize!
    v = Dict(:el => 0.4, :hola => 0.2, :mundo => 0.4) |> normalize!
    w = Dict(:xel => 0.4, :xhola => 0.2, :xmundo => 0.4) |> normalize!

    dist = angle_distance
    @test dist(u, v) ≈ 0.5975474808029686
    @test dist(u, u) <= eps(Float32)
    @test dist(w, u) ≈ 1.5707963267948966
end

@testset "operations" begin
    u = Dict(:el => 0.1, :hola => 0.2, :mundo => 0.4)
    v = Dict(:el => 0.2, :hola => 0.4, :mundo => 0.8)
    w = Dict(:el => 0.1^2, :hola => 0.2^2, :mundo => 0.4^2)
    y = Dict(:el => 0.1/9, :hola => 0.2/9, :mundo => 0.4/9)
    @test u == u
    @test u != v
    @test u + u == v
    @test u * u == w
    @test u * (1/9) == y
    @test (1/9) * u == y
    @test dot(normalize!(u + v - v), normalize!(u)) > 0.99
end

@testset "io" begin
    buff = IOBuffer("""{"key1": "value1a", "key2c": "value2a"}
{"key1": "value1b", "key2c": "value2b"}
{"key1": "value1c", "key2b": "value2c"}
{"key1": "value1d", "key2a": "value2d"}""")
    itertweets(buff) do x
        @info x
    end
end

_corpus = [
    "la casa roja",
    "la casa verde",
    "la casa azul",
    "la manzana roja",
    "la pera verde esta rica",
    "la manzana verde esta rica",
    "la hoja verde",
]

@testset "transpose bow" begin
    config = TextConfig()
    config.nlist = [1]
    config.qlist = []
    config.slist = []
    model = fit(VectorModel, config, _corpus)
    X = [vectorize(model, FreqModel, x) for x in _corpus]
    dX = transpose(X)
end

@testset "invindex" begin
    config = TextConfig()
    #config.qlist = [3, 4]
    #config.nlist = [1, 2, 3]
    config.nlist = [1]

    model = fit(VectorModel, config, _corpus)
    invindex = fit(InvIndex, [vectorize(model, TfidfModel, text) for text in _corpus])
    q = vectorize(model, TfidfModel, "la casa roja")
    res = search(invindex, cosine_distance, q, KnnResult(4))
    ires = [r.objID for r in res]
    @test sort(ires) == [1, 2, 3, 4]
    shortindex = prune(invindex, 3)
    @show shortindex
    res = search(shortindex, cosine_distance, q, KnnResult(4))
    ires = [r.objID for r in res]
    @test sort(ires) == [1, 2, 3, 4]
end

@testset "centroid computing" begin
    config = TextConfig()
    config.nlist = [1]
    config.qlist = []
    config.slist = []
    model = fit(VectorModel, config, _corpus)
    X = [vectorize(model, FreqModel, x) for x in _corpus]
    x = sum(X) |> normalize!
    vec = bow(model, x)
    expected = Dict(:la => 0.7366651330405098,:verde => 0.39921969741172364,:azul => 0.11248181187626208,:pera => 0.08712803682959973,:esta => 0.17425607365919946,:roja => 0.22496362375252416,:hoja => 0.11248181187626208,:casa => 0.33744543562878626,:rica => 0.17425607365919946,:manzana => 0.19960984870586182)
    @test 0.999 < dot(vec, expected)
end

@testset "rocchio" begin
    config = TextConfig()
    config.nlist = [1,2]
    config.qlist = [3]
    config.slist = []
    corpus = [x[1] for x in labeled_corpus]

    model = fit(VectorModel, config, corpus)
    X = [vectorize(model, TfidfModel, x) for x in corpus]
    y = [x[2] for x in labeled_corpus]
    rocchio = fit(Rocchio, X, y)
    @test sum(predict.(rocchio, X) .== y)/length(y) == 1.0

    for p in transform.(rocchio, X)
        println(p)
    end
end

@testset "rocchio bagging" begin
    config = TextConfig()
    config.nlist = [1, 2]
    config.qlist = [3]
    config.slist = []
    corpus = [x[1] for x in labeled_corpus]
    model = fit(VectorModel, config, corpus)
    X = [vectorize(model, TfidfModel, x) for x in corpus]
    y = [x[2] for x in labeled_corpus]
    rocchio = fit(RocchioBagging, X, y)
    for p in transform.(rocchio, X)
        println("===>", p)
    end

    @test mean(predict.(rocchio, X) .== y) > 0.5
end

@testset "neardup" begin
    config = TextConfig()
    config.nlist = [1]
    config.qlist = []
    config.slist = []

    function create_corpus()
        alphabet = Char.(97:100)
        corpus = String[]
        L = Char[]
        for i in 1:10000
            resize!(L, 0)
            for j in 1:7
                append!(L, rand(alphabet, 2))
                push!(L, ' ')
            end

            push!(corpus, join(L))
        end

        corpus
    end

    corpus = create_corpus()
    model = fit(VectorModel, config, corpus)
    @show corpus[1:10]
    X = [vectorize(model, TfModel, x) for x in corpus]
    L, D  = neardup(X, 0.2)
    @test length(X) > length(unique(L))
end


text:
export ChebyshevT

"""
    ChebyshevT{T<:Number, X}(coeffs::AbstractVector)

Chebyshev polynomial of the first kind.

Construct a polynomial from its coefficients `coeffs`, lowest order first, optionally in
terms of the given variable `var`, which can be a character, symbol, or string.

!!! note
    `ChebyshevT` is not axis-aware, and it treats `coeffs` simply as a list of coefficients with the first
    index always corresponding to the coefficient of `T_0(x)`.

# Examples

```jldoctest ChebyshevT
julia> using Polynomials

julia> p = ChebyshevT([1, 0, 3, 4])
ChebyshevT(1⋅T_0(x) + 3⋅T_2(x) + 4⋅T_3(x))

julia> ChebyshevT([1, 2, 3, 0], :s)
ChebyshevT(1⋅T_0(s) + 2⋅T_1(s) + 3⋅T_2(s))

julia> one(ChebyshevT)
ChebyshevT(1.0⋅T_0(x))

julia> p(0.5)
-4.5

julia> Polynomials.evalpoly(5.0, p, false) # bypasses the domain check done in p(5.0)
2088.0
```

The latter shows how to evaluate a `ChebyshevT` polynomial outside of its domain, which is `[-1,1]`. (For newer versions of `Julia`, `evalpoly` is an exported function from Base with methods extended in this package, so the module qualification is unnecessary.

!!! Note:
    The Chebyshev polynomials are also implemented in `ApproxFun`, `ClassicalOrthogonalPolynomials.jl`, `FastTransforms.jl`, and `SpecialPolynomials.jl`.

"""
struct ChebyshevT{T <: Number, X} <: AbstractPolynomial{T, X}
    coeffs::Vector{T}
    function ChebyshevT{T, X}(coeffs::AbstractVector{S}) where {T <: Number,X, S}

        if Base.has_offset_axes(coeffs)
            @warn "ignoring the axis offset of the coefficient vector"
        end

        N = findlast(!iszero, coeffs)
        N == nothing && return new{T,X}(zeros(T,1))
        cs = T[coeffs[i] for i ∈ firstindex(coeffs):N]
        new{T,X}(cs)
    end
end

@register ChebyshevT

function Base.convert(P::Type{<:Polynomial}, ch::ChebyshevT)
    T = eltype(P)
    X = indeterminate(P,ch)
    Q = ⟒(P){T,X}

    if length(ch) < 3
        return Q(ch.coeffs)
    end

    c0 = Q(ch[end - 1])
    c1 = Q(ch[end])
    x = variable(Q)
    @inbounds for i in degree(ch):-1:2
        tmp = c0
        c0 = Q(ch[i - 2]) - c1
        c1 = tmp + c1 * x * 2
    end
    return c0 + c1 * x
end
Base.convert(C::Type{<:ChebyshevT}, p::Polynomial) = p(variable(C))

domain(::Type{<:ChebyshevT}) = Interval(-1, 1)
function Base.one(::Type{P}) where {P<:ChebyshevT}
    T,X = eltype(P), indeterminate(P)
    ChebyshevT{T,X}(ones(T,1))
end
function variable(::Type{P}) where {P<:ChebyshevT}
    T,X = eltype(P), indeterminate(P)
    ChebyshevT{T,X}([zero(T), one(T)])
end
constantterm(p::ChebyshevT) = p(0)
"""
    (::ChebyshevT)(x)

Evaluate the Chebyshev polynomial at `x`. If `x` is outside of the domain of [-1, 1], an error will be thrown. The evaluation uses Clenshaw Recursion.

# Examples
```jldoctest ChebyshevT
julia> using Polynomials

julia> c = ChebyshevT([2.5, 1.5, 1.0])
ChebyshevT(2.5⋅T_0(x) + 1.5⋅T_1(x) + 1.0⋅T_2(x))

julia> c(0)
1.5

julia> c.(-1:0.5:1)
5-element Vector{Float64}:
 2.0
 1.25
 1.5
 2.75
 5.0
```
"""
function evalpoly(x::S, ch::ChebyshevT{T}) where {T,S}
    x ∉ domain(ch) && throw(ArgumentError("$x outside of domain"))
    evalpoly(x, ch, false)
end

# no checking, so can be called directly through any third argument
function evalpoly(x::S, ch::ChebyshevT{T}, checked) where {T,S}
    R = promote_type(T, S)
    length(ch) == 0 && return zero(R)
    length(ch) == 1 && return R(ch[0])
    c0 = ch[end - 1]
    c1 = ch[end]
    @inbounds for i in lastindex(ch) - 2:-1:0
        c0, c1 = ch[i] - c1, c0 + c1 * 2x
    end
    return R(c0 + c1 * x)
end

function vander(P::Type{<:ChebyshevT}, x::AbstractVector{T}, n::Integer) where {T <: Number}
    A = Matrix{T}(undef, length(x), n + 1)
    A[:, 1] .= one(T)
    if n > 0
        A[:, 2] .= x
        @inbounds for i in 3:n + 1
            A[:, i] .= A[:, i - 1] .* 2x .- A[:, i - 2]
        end
    end
    return A
end

function integrate(p::ChebyshevT{T,X}) where {T,X}
    R = eltype(one(T) / 1)
    Q = ChebyshevT{R,X}
    if hasnan(p)
        return Q([NaN])
    end
    n = length(p)
    if n == 1
        return Q([zero(R), p[0]])
    end
    a2 = Vector{R}(undef, n + 1)
    a2[1] = zero(R)
    a2[2] = p[0]
    a2[3] = p[1] / 4
    @inbounds for i in 2:n - 1
        a2[i + 2] = p[i] / (2 * (i + 1))
        a2[i] -= p[i] / (2 * (i - 1))
    end

    return Q(a2)
end

function derivative(p::ChebyshevT{T}, order::Integer = 1) where {T}
    order < 0 && throw(ArgumentError("Order of derivative must be non-negative"))
    R  = eltype(one(T)/1)
    order == 0 && return convert(ChebyshevT{R}, p)
    hasnan(p) && return ChebyshevT(R[NaN], indeterminate(p))
    order > length(p) && return zero(ChebyshevT{R})

    q =  convert(ChebyshevT{R}, copy(p))
    n = length(p)
    der = Vector{R}(undef, n)

    for j in n:-1:2
        der[j] = 2j * q[j]
        q[j - 2] += j * q[j] / (j - 2)
    end
    if n > 1
        der[2] = 4q[2]
    end
    der[1] = q[1]

    pp = ChebyshevT(der, indeterminate(p))
    return order > 1 ?  derivative(pp, order - 1) : pp

end

function companion(p::ChebyshevT{T}) where T
    d = length(p) - 1
    d < 1 && throw(ArgumentError("Series must have degree greater than 1"))
    d == 1 && return diagm(0 => [-p[0] / p[1]])
    R = eltype(one(T) / one(T))

    scl = vcat(1.0, fill(R(√0.5), d - 1))

    diag = vcat(√0.5, fill(R(0.5), d - 2))
    comp = diagm(1 => diag,
                 -1 => diag)
    monics = p.coeffs ./ p.coeffs[end]
    comp[:, end] .-= monics[1:d] .* scl ./ scl[end] ./ 2
    return R.(comp)
end

# scalar +
function Base.:+(p::ChebyshevT{T,X}, c::S) where {T,X, S<:Number}
    R = promote_type(T,S)
    cs = collect(R, values(p))
    cs[1] += c
    ChebyshevT{R,X}(cs)
end
function Base.:+(p::P, c::T) where {T,X,P<:ChebyshevT{T,X}}
    cs = collect(T, values(p))
    cs[1] += c
    P(cs)
end

function Base.:+(p1::ChebyshevT{T,X}, p2::ChebyshevT{T,X}) where {T,X}
    n = max(length(p1), length(p2))
    c = T[p1[i] + p2[i] for i = 0:n]
    return ChebyshevT{T,X}(c)
end

function Base.:*(p1::ChebyshevT{T,X}, p2::ChebyshevT{T,X}) where {T,X}
    z1 = _c_to_z(p1.coeffs)
    z2 = _c_to_z(p2.coeffs)
    prod = fastconv(z1, z2)
    cs = _z_to_c(prod)
    ret = ChebyshevT(cs,X)
    return truncate!(ret)
end

function Base.divrem(num::ChebyshevT{T,X}, den::ChebyshevT{S,Y}) where {T,X,S,Y}
    assert_same_variable(num, den)
    n = length(num) - 1
    m = length(den) - 1

    R = typeof(one(T) / one(S))
    P = ChebyshevT{R,X}

    if n < m
        return zero(P), convert(P, num)
    elseif m == 0
        den[0] ≈ 0 && throw(DivideError())
        return num ./ den[end], zero(P)
    end

    znum = _c_to_z(num.coeffs)
    zden = _c_to_z(den.coeffs)
    quo, rem = _z_division(znum, zden)
    q_coeff = _z_to_c(quo)
    r_coeff = _z_to_c(rem)
    return P(q_coeff), P(r_coeff)
end

function showterm(io::IO, ::Type{ChebyshevT{T,X}}, pj::T, var, j, first::Bool, mimetype) where {N, T,X}
    iszero(pj) && return false
    !first &&  print(io, " ")
    print(io, hasneg(T)  && isneg(pj) ? "- " :  (!first ? "+ " : ""))
    print(io, "$(abs(pj))⋅T_$j($var)")
    return true
end

#=
zseries =#

function _c_to_z(cs::AbstractVector{T}) where {T}
    n = length(cs)
    U = typeof(one(T) / 2)
    zs = zeros(U, 2n - 1)
    zs[n:end] = cs ./ 2
    return zs .+ reverse(zs)
end

function _z_to_c(z::AbstractVector{T}) where {T}
    n = (length(z) + 1) ÷ 2
    cs = z[n:end]
    cs[2:n] *= 2
    return cs
end

function _z_division(z1::AbstractVector{T}, z2::AbstractVector{S}) where {T,S}
    R = eltype(one(T) / one(S))
    length(z1)
    length(z2)
    if length(z2) == 1
        z1 ./= z2
        return z1, zero(R)
    elseif length(z1) < length(z2)
        return zero(R), R.(z1)
    end
    dlen = length(z1) - length(z2)
    scl = z2[1]
    z2 ./= scl
    quo = Vector{R}(undef, dlen + 1)
    i = 1
    j = dlen + 1
    while i < j
        r = z1[i]
        quo[i] = z1[i]
        quo[end - i + 1] = r
        tmp = r .* z2
        z1[i:i + length(z2) - 1] .-= tmp
        z1[j:j + length(z2) - 1] .-= tmp
        i += 1
        j -= 1
    end

    r = z1[i]
    quo[i] = r
    tmp = r * z2
    z1[i:i + length(z2) - 1] .-= tmp
    quo ./= scl
    rem = z1[i + 1:i - 2 + length(z2)]
    return quo, rem
end

ext: jl


text:
"""
The `Downloads` module exports a function [`download`](@ref), which provides cross-platform, multi-protocol,
in-process download functionality implemented with [libcurl](https://curl.haxx.se/libcurl/).   It is used
for the `Base.download` function in Julia 1.6 or later.

More generally, the module exports functions and types that provide lower-level control and diagnostic information
for file downloading:
- [`download`](@ref) — download a file from a URL, erroring if it can't be downloaded
- [`request`](@ref) — request a URL, returning a `Response` object indicating success
- [`Response`](@ref) — a type capturing the status and other metadata about a request
- [`RequestError`](@ref) — an error type thrown by `download` and `request` on error
- [`Downloader`](@ref) — an object encapsulating shared resources for downloading
"""
module Downloads

using Base.Experimental: @sync
using NetworkOptions
using ArgTools

include("Curl/Curl.jl")
using .Curl

export download, request, Downloader, Response, RequestError

## public API types ##

"""
    Downloader(; [ grace::Real = 30 ])

`Downloader` objects are used to perform individual `download` operations.
Connections, name lookups and other resources are shared within a `Downloader`.
These connections and resources are cleaned up after a configurable grace period
(default: 30 seconds) since anything was downloaded with it, or when it is
garbage collected, whichever comes first. If the grace period is set to zero,
all resources will be cleaned up immediately as soon as there are no more
ongoing downloads in progress. If the grace period is set to `Inf` then
resources are not cleaned up until `Downloader` is garbage collected.
"""
mutable struct Downloader
    multi::Multi
    ca_roots::Union{String, Nothing}
    easy_hook::Union{Function, Nothing}

    Downloader(multi::Multi) = new(multi, get_ca_roots(), EASY_HOOK[])
end
Downloader(; grace::Real=30) = Downloader(Multi(grace_ms(grace)))

function grace_ms(grace::Real)
    grace < 0 && throw(ArgumentError("grace period cannot be negative: $grace"))
    grace <= typemax(UInt64) ÷ 1000 ? round(UInt64, 1000*grace) : typemax(UInt64)
end

function easy_hook(downloader::Downloader, easy::Easy, info::NamedTuple)
    downloader.easy_hook !== nothing && downloader.easy_hook(easy, info)
end

get_ca_roots() = Curl.SYSTEM_SSL ? ca_roots() : ca_roots_path()

function set_ca_roots(downloader::Downloader, easy::Easy)
    ca_roots = downloader.ca_roots
    ca_roots !== nothing && set_ca_roots_path(easy, ca_roots)
end

const DOWNLOAD_LOCK = ReentrantLock()
const DOWNLOADER = Ref{Union{Downloader, Nothing}}(nothing)
const EASY_HOOK = Ref{Union{Function, Nothing}}(nothing)

"""
    struct Response
        proto   :: String
        url     :: String
        status  :: Int
        message :: String
        headers :: Vector{Pair{String,String}}
    end

`Response` is a type capturing the properties of a successful response to a
request as an object. It has the following fields:

- `proto`: the protocol that was used to get the response
- `url`: the URL that was ultimately requested after following redirects
- `status`: the status code of the response, indicating success, failure, etc.
- `message`: a textual message describing the nature of the response
- `headers`: any headers that were returned with the response

The meaning and availability of some of these responses depends on the protocol
used for the request. For many protocols, including HTTP/S and S/FTP, a 2xx
status code indicates a successful response. For responses in protocols that do
not support headers, the headers vector will be empty. HTTP/2 does not include a
status message, only a status code, so the message will be empty.
"""
struct Response
    proto   :: Union{String, Nothing}
    url     :: String # redirected URL
    status  :: Int
    message :: String
    headers :: Vector{Pair{String,String}}
end

"""
    struct RequestError <: ErrorException
        url      :: String
        code     :: Int
        message  :: String
        response :: Response
    end

`RequestError` is a type capturing the properties of a failed response to a
request as an exception object:

- `url`: the original URL that was requested without any redirects
- `code`: the libcurl error code; `0` if a protocol-only error occurred
- `message`: the libcurl error message indicating what went wrong
- `response`: response object capturing what response info is available

The same `RequestError` type is thrown by `download` if the request was
successful but there was a protocol-level error indicated by a status code that
is not in the 2xx range, in which case `code` will be zero and the `message`
field will be the empty string. The `request` API only throws a `RequestError`
if the libcurl error `code` is non-zero, in which case the included `response`
object is likely to have a `status` of zero and an empty message. There are,
however, situations where a curl-level error is thrown due to a protocol error,
in which case both the inner and outer code and message may be of interest.
"""
struct RequestError <: Exception
    url      :: String # original URL
    code     :: Int
    message  :: String
    response :: Response
end

function Base.showerror(io::IO, err::RequestError)
    print(io, "$(error_message(err)) while requesting $(err.url)")
end

function error_message(err::RequestError)
    errstr = err.message
    status = err.response.status
    message = err.response.message
    status_re = Regex(status == 0 ? "" : "\\b$status\\b")

    err.code == Curl.CURLE_OK &&
        return isempty(message) ? "Error status $status" :
            contains(message, status_re) ? message :
                "$message (status $status)"

    isempty(message) && !isempty(errstr) &&
        return status == 0 ? errstr : "$errstr (status $status)"

    isempty(message) && (message = "Status $status")
    isempty(errstr)  && (errstr = "curl error $(err.code)")

    !contains(message, status_re) && !contains(errstr, status_re) &&
        (errstr = "status $status; $errstr")

    return "$message ($errstr)"
end

## download API ##

"""
    download(url, [ output = tempfile() ];
        [ method = "GET", ]
        [ headers = <none>, ]
        [ timeout = <none>, ]
        [ progress = <none>, ]
        [ verbose = false, ]
        [ downloader = <default>, ]
    ) -> output

        url        :: AbstractString
        output     :: Union{AbstractString, AbstractCmd, IO}
        method     :: AbstractString
        headers    :: Union{AbstractVector, AbstractDict}
        timeout    :: Real
        progress   :: (total::Integer, now::Integer) --> Any
        verbose    :: Bool
        downloader :: Downloader

Download a file from the given url, saving it to `output` or if not specified, a
temporary path. The `output` can also be an `IO` handle, in which case the body
of the response is streamed to that handle and the handle is returned. If
`output` is a command, the command is run and output is sent to it on stdin.

If the `downloader` keyword argument is provided, it must be a `Downloader`
object. Resources and connections will be shared between downloads performed by
the same `Downloader` and cleaned up automatically when the object is garbage
collected or there have been no downloads performed with it for a grace period.
See `Downloader` for more info about configuration and usage.

If the `headers` keyword argument is provided, it must be a vector or dictionary
whose elements are all pairs of strings. These pairs are passed as headers when
downloading URLs with protocols that supports them, such as HTTP/S.

The `timeout` keyword argument specifies a timeout for the download in seconds,
with a resolution of milliseconds. By default no timeout is set, but this can
also be explicitly requested by passing a timeout value of `Inf`.

If the `progress` keyword argument is provided, it must be a callback funtion
which will be called whenever there are updates about the size and status of the
ongoing download. The callback must take two integer arguments: `total` and
`now` which are the total size of the download in bytes, and the number of bytes
which have been downloaded so far. Note that `total` starts out as zero and
remains zero until the server gives an indication of the total size of the
download (e.g. with a `Content-Length` header), which may never happen. So a
well-behaved progress callback should handle a total size of zero gracefully.

If the `verbose` option is set to true, `libcurl`, which is used to implement
the download functionality will print debugging information to `stderr`.
"""
function download(
    url        :: AbstractString,
    output     :: Union{ArgWrite, Nothing} = nothing;
    method     :: Union{AbstractString, Nothing} = nothing,
    headers    :: Union{AbstractVector, AbstractDict} = Pair{String,String}[],
    timeout    :: Real = Inf,
    progress   :: Union{Function, Nothing} = nothing,
    verbose    :: Bool = false,
    downloader :: Union{Downloader, Nothing} = nothing,
) :: ArgWrite
    arg_write(output) do output
        response = request(
            url,
            output = output,
            method = method,
            headers = headers,
            timeout = timeout,
            progress = progress,
            verbose = verbose,
            downloader = downloader,
        )::Response
        status_ok(response.proto, response.status) && return output
        throw(RequestError(url, Curl.CURLE_OK, "", response))
    end
end

## request API ##

"""
    request(url;
        [ input = <none>, ]
        [ output = <none>, ]
        [ method = input ? "PUT" : output ? "GET" : "HEAD", ]
        [ headers = <none>, ]
        [ timeout = <none>, ]
        [ progress = <none>, ]
        [ verbose = false, ]
        [ throw = true, ]
        [ downloader = <default>, ]
    ) -> Union{Response, RequestError}

        url        :: AbstractString
        input      :: Union{AbstractString, AbstractCmd, IO}
        output     :: Union{AbstractString, AbstractCmd, IO}
        method     :: AbstractString
        headers    :: Union{AbstractVector, AbstractDict}
        timeout    :: Real
        progress   :: (dl_total, dl_now, ul_total, ul_now) --> Any
        verbose    :: Bool
        throw      :: Bool
        downloader :: Downloader

Make a request to the given url, returning a `Response` object capturing the
status, headers and other information about the response. The body of the
reponse is written to `output` if specified and discarded otherwise. For HTTP/S
requests, if an `input` stream is given, a `PUT` request is made; otherwise if
an `output` stream is given, a `GET` request is made; if neither is given a
`HEAD` request is made. For other protocols, appropriate default methods are
used based on what combination of input and output are requested. The following
options differ from the `download` function:

- `input` allows providing a request body; if provided default to `PUT` request
- `progress` is a callback taking four integers for upload and download progress
- `throw` controls whether to throw or return a `RequestError` on request error

Note that unlike `download` which throws an error if the requested URL could not
be downloaded (indicated by non-2xx status code), `request` returns a `Response`
object no matter what the status code of the response is. If there is an error
with getting a response at all, then a `RequestError` is thrown or returned.
"""
function request(
    url        :: AbstractString;
    input      :: Union{ArgRead, Nothing} = nothing,
    output     :: Union{ArgWrite, Nothing} = nothing,
    method     :: Union{AbstractString, Nothing} = nothing,
    headers    :: Union{AbstractVector, AbstractDict} = Pair{String,String}[],
    timeout    :: Real = Inf,
    progress   :: Union{Function, Nothing} = nothing,
    verbose    :: Bool = false,
    throw      :: Bool = true,
    downloader :: Union{Downloader, Nothing} = nothing,
) :: Union{Response, RequestError}
    lock(DOWNLOAD_LOCK) do
        yield() # let other downloads finish
        downloader isa Downloader && return
        while true
            downloader = DOWNLOADER[]
            downloader isa Downloader && return
            DOWNLOADER[] = Downloader()
        end
    end
    local response
    have_input = input !== nothing
    have_output = output !== nothing
    input = something(input, devnull)
    output = something(output, devnull)
    input_size = arg_read_size(input)
    progress = p_func(progress, input, output)
    arg_read(input) do input
        arg_write(output) do output
            with_handle(Easy()) do easy
                # setup the request
                set_url(easy, url)
                set_timeout(easy, timeout)
                set_verbose(easy, verbose)
                add_headers(easy, headers)

                # libcurl does not set the default header reliably so set it
                # explicitly unless user has specified it, xref
                # https://github.com/JuliaLang/Pkg.jl/pull/2357
                if !any(kv -> lowercase(kv[1]) == "user-agent", headers)
                    Curl.add_header(easy, "User-Agent", Curl.USER_AGENT)
                end

                if have_input
                    enable_upload(easy)
                    if input_size !== nothing
                        set_upload_size(easy, input_size)
                    end
                    if applicable(seek, input, 0)
                        set_seeker(easy) do offset
                            seek(input, Int(offset))
                        end
                    end
                else
                    set_body(easy, have_output && method != "HEAD")
                end
                method !== nothing && set_method(easy, method)
                progress !== nothing && enable_progress(easy)
                set_ca_roots(downloader, easy)
                info = (url = url, method = method, headers = headers)
                easy_hook(downloader, easy, info)

                # do the request
                add_handle(downloader.multi, easy)
                try # ensure handle is removed
                    @sync begin
                        @async for buf in easy.output
                            write(output, buf)
                        end
                        if progress !== nothing
                            @async for prog in easy.progress
                                progress(prog...)
                            end
                        end
                        if have_input
                            @async upload_data(easy, input)
                        end
                    end
                finally
                    remove_handle(downloader.multi, easy)
                end

                # return the response or throw an error
                response = Response(get_response_info(easy)...)
                easy.code == Curl.CURLE_OK && return response
                message = get_curl_errstr(easy)
                response = RequestError(url, easy.code, message, response)
                throw && Base.throw(response)
            end
        end
    end
    return response
end

## helper functions ##

function p_func(progress::Function, input::ArgRead, output::ArgWrite)
    hasmethod(progress, NTuple{4,Int}) && return progress
    hasmethod(progress, NTuple{2,Int}) ||
        throw(ArgumentError("invalid progress callback"))

    input === devnull && output !== devnull &&
        return (total, now, _, _) -> progress(total, now)
    input !== devnull && output === devnull &&
        return (_, _, total, now) -> progress(total, now)

    (dl_total, dl_now, ul_total, ul_now) ->
        progress(dl_total + ul_total, dl_now + ul_now)
end
p_func(progress::Nothing, input::ArgRead, output::ArgWrite) = nothing

arg_read_size(path::AbstractString) = filesize(path)
arg_read_size(io::Base.GenericIOBuffer) = bytesavailable(io)
arg_read_size(::Base.DevNull) = 0
arg_read_size(::Any) = nothing

end # module


text:
struct PendingInvocation
    type::CLRObject
    target::CLRObject
    name::Union{Symbol,Nothing}
end

struct GenericModifier{S}
    subject::S
    params
end

function Base.getindex(pi::PendingInvocation, params...)
    GenericModifier(pi, params)
end

function (call::PendingInvocation)(args...)
    try
        if !isnothing(call.name)
            invokemember(call.type, call.target, call.name, args...)
        else
            invokemember(BindingFlags.CreateInstance, call.type, CLRObject(0), "", args...)
        end
    catch ex
        if ex isa CLRException && isclrtype(ex.object, T"System.Reflection.AmbiguousMatchException")
            rethrow(ErrorException("""Multiple overloads match the given binding criteria.
            Use '$(call.name)[T]()' to call a generic method with a type argument T."""))
        end
        rethrow()
    end
end

function (gencall::GenericModifier{PendingInvocation})(args...)
    type = gencall.subject.type
    name = gencall.subject.name
    ngenarg = length(gencall.params)
    if isnothing(name)
        genty = makegenerictype(type, gencall.params...)
        return invokemember(BindingFlags.CreateInstance, genty, CLRObject(0), "", args...)
    end
    candidates = getmember(type, name)
    for m in candidates
        genargs = invokemember(m, :GetGenericArguments)
        len = invokemember(genargs, :Length)
        if len == ngenarg
            resolved = if ngenarg != 0
                makegenericmethod(m, gencall.params...)
            else
                m
            end
            arr = arrayof(T"System.Object", length(args))
            for (i, x) in enumerate(args)
                arraystore(arr, i - 1, x)
            end
            return invokemember(resolved, :Invoke, gencall.subject.target, arr)
        end
    end
    error("No method named '$name' that takes $ngenarg generic arguments in type $type")
end

function Base.getproperty(obj::CLRObject, sym::Symbol)
    ty = clrtypeof(obj)
    members = nothing
    pi = nothing
    if isassignable(T"System.Type", ty)
        if sym == :new
            return PendingInvocation(obj, CLRObject(0), nothing)
        end
        members = getmember(obj, sym)
        if isempty(members)
            members = getmember(ty, sym)
            if isempty(members)
                error("No such member '$sym' in type $obj")
            end
            pi = PendingInvocation(ty, obj, sym)

        else
            pi = PendingInvocation(obj, CLRObject(0), sym)
        end
    else
        members = getmember(ty, sym)
        if isempty(members)
            error("No such member '$sym' in type $ty")
        end
        pi = PendingInvocation(ty, obj, sym)
    end
    return if any(x->isassignable(T"System.Reflection.MethodInfo", clrtypeof(x)), members)
        pi
    else
        invokemember(pi.type, pi.target, pi.name)
    end
end

function Base.setproperty!(obj::CLRObject, sym::Symbol, val)
    ty = clrtypeof(obj)
    flags = BindingFlags.Static | BindingFlags.Instance | BindingFlags.Public | BindingFlags.SetProperty
    if isassignable(T"System.Type", ty)
        invokemember(flags, obj, CLRObject(0), sym, val)
    else
        invokemember(flags, ty, obj, sym, val)
    end
end


text:
using Lime
using Test

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


text:
export VQEtrain

"""
    VQEtrain(set::Union{Array{Int64,1}, Array{Float64,1}};
            depth=4, niter=100, autoTrain=true, niterMax = 5000, Optimizer=:NADAM,
            StopTH = 2e-9, PerturbTH = 1e-7, Threshold = 0.5, PurtAmp = 0.05, ConvergeTH=1e-4)
    ->
    [Hvalues, isPerturbed, dEs, iConverge, EConverge, ECmin, ECmm]
The training function for VQE4NPP.
    \n== Arguments ==
    \n`set::Union{Array{Int64,1}, Array{Float64,1}}`: The target number set.
    \n`depth::Int64`: The depth (number of layers) of the quantum differentiable circuit.
    \n`niter::Int64`: The number of ierations.
    \n`autoTrain::Bool`: The option to enable auto training iterations.
    \n`niterMax::Float64`: The maximal number of ierations when applying `autoTain`.
    \n`Optimizer::Function`: The optimizer function for Gradient Descent. The default optimizer is NADAM from package Flux.jl and can be replaced by user-defined function in forms of `(args...)->f(args...)`.
    \n`StopTH::Float64`: The threshold of the derivative of the training curve to stop training.
    \n`PerturbTH::Float64`: The threshold of the derivative of the training curve to apply perturbations.
    \n`Threshold::Float64`: The threshold for determining whether the <H> is small enough.
    \n`PurtAmp::Float64`: The amplitude for perurbation.
    \n`ConvergeTH::Float64`: The threshold of the derivative of the training curve to determine whther the training curve has converged.
    \n== Output variables ==
    \n`Hvalues::Array{Float64,1}`: The <H> of each training step.
    \n`isPerturbed::Array{Float64,1}`: The mark of whether the <H> of ith step is gained under the condition where the differentiable circuit is perturbed. 1 means perturbed.
    \n`dEs::Array{Float64,1}`: The derivative of the training curve (<H>) of each training step.
    \n`iConverge::Int64`: The step when the curve converges.
    \n`EConverge::Float64`: The convergent energy corresponding to `iConverge`.
    \n`ECmin::Float64`: After the curve starts converging, the minimal value of <H> before the training stops.
    \n`ECmm::Float64`: A mean value of ECmin and its adjacent training values.
"""
function VQEtrain(set::Union{Array{Int64,1}, Array{Float64,1}};
                  depth::Int64=4,
                  niter::Int64=100,
                  autoTrain::Bool=true,
                  niterMax::Int64=5000,
                  Optimizer::Function=()->NADAM(),
                  StopTH::Float64=1e-8,
                  PerturbTH::Float64=1e-7,
                  Threshold::Float64=0.5,
                  PurtAmp::Float64=0.02,
                  ConvergeTH::Float64=1e-3,
                  showSteps::Bool=true)

    H = HofNPP(set)
    n = nqubits(H)
    i = 0
    k = 0
    h = 5
    di = 2
    sign = 1
    iConverge = 0
    EConverge = 0
    Hvalues = Float64[]
    len2Grad = 2h+di
    dEs = repeat([NaN],len2Grad)
    isPerturbed = zeros(len2Grad*2)

    circuit = variational_circuit(n,depth)
    dispatch!(circuit, :random)
    if autoTrain == true
        iDisplay = (i) -> i%200 == 0
    else
        iDisplay = (i) -> 70*i%niter == 0
    end
    while i < len2Grad*2
        _, grad = expect'(H, zero_state(n)=>circuit)
        pars = parameters(circuit)
        dispatch!(circuit, Flux.Optimise.update!(Optimizer(), pars, grad))
        append!(Hvalues, real.(expect(H, zero_state(n)=>circuit)))
        i = i + 1
    end
    while i < niter
        _, grad = expect'(H, zero_state(n)=>circuit)
        pars = parameters(circuit)
        dispatch!(circuit, Flux.Optimise.update!(Optimizer(), pars, grad))
        append!(Hvalues, real.(expect(H, zero_state(n)=>circuit)))
        append!(isPerturbed, floor(k/100))
        i=i+1
        cg, __ = curveGrad(i-len2Grad,Hvalues, h=5, di=2)
        append!(dEs, cg)
        iDisplay(i) && showSteps == true && println("Step $i, E = $(round(Hvalues[i], sigdigits=6)), dE = $(round(cg, sigdigits=6))")
        if iConverge == 0
            if abs(cg) <= ConvergeTH
                iConverge = i-len2Grad
                EConverge = Hvalues[iConverge]
                showSteps == true && println("Started to Coverge at Step $i,E = $(round(Hvalues[i], sigdigits=6))")
            end
        else
            if abs(cg) <= PerturbTH && k == 0 && Hvalues[i] > sum(set)%2 + Threshold
                t = sum(Hvalues[i-19:i-10]) / sum(Hvalues[i-9:i])-1
                if t > 0
                    sign = -sign
                end
                showSteps == true && println("Perturbation Activated")
                RXs = collect_blocks(RotationGate{1,Float64,XGate},circuit)
                dispatch!.(-,RXs,sign*PurtAmp*ones(length(RXs)))
                k = k + 100
            else
                k>0 && (k = k-1)
            end
        end
        if i == niter && autoTrain == true
            if Hvalues[i] > sum(set)%2 + Threshold && abs(cg) > StopTH && i < niterMax
                niter = niter + 100
            end
        end
    end
    append!(dEs,repeat([NaN],len2Grad))
    println("Total steps: $(length(Hvalues)), Final E: $(round(Hvalues[end], sigdigits=6))")
    ECmin, imin = findmin(Hvalues)
    ECmm = mean(Hvalues[imin:end])
    [Hvalues, isPerturbed, dEs, iConverge, EConverge, ECmin, ECmm]
end


text:
@testset "1011.capacity-to-ship-packages-within-d-days.jl" begin
    @test ship_within_days([1,2,3,1,1], 4) == 3
    @test ship_within_days([3,2,2,4,1,4], 3) == 6
    @test ship_within_days([1,2,3,4,5,6,7,8,9,10], 5) == 15
end


text:
""" Collection of tracer density models.

[`SersicModel`](@ref) : Representation of a Sersic surface density profile
"""
module Tracers

import Slomo.Models: DensityModel, mass, density, density2d
import Slomo.Utils: gamma, gamma_inc

export SersicModel

"""
    b_cb(n)

Compute the "b" parameter in the Sersic function, from the Ciotti & Bertin (1999)
approximation.  `n` is the Sersic index.
"""
function b_cb(n)
    -1. / 3 + 2. * n + 4 / (405. * n) + 46 / (25515. * n^2)
end

"""
    p_ln(n)

Compute the "p" parameter in fitting deprojected Sersic function, from
Lima Neto et al. (1999).  `n` is the Sersic index
"""
function p_ln(n)
    1. - 0.6097 / n + 0.05463 / n^2
end

"""
    s_sersic(R, Re, n, Mtot)

Compute the surface density profile for Sersic model.

* `R`: projected (2d) radii in kpc
* `Re`: effective radius in kpc
* `n`: index
* `Mtot`: total mass
"""
function s_sersic(R, Re, n, Mtot)
    bn = b_cb(n)
    s0 = Mtot * bn ^ (2 * n) / (2π * n * Re ^ 2 * gamma(2 * n))
    return s0 * @. exp(-bn * (R / Re) ^ (1.0 / n))
end

"""
    rho_sersic(r, Re, n, Mtot)

Compute the volume density profile for Sersic model.

* `r`: deprojected (3d) radii in kpc
* `Re`: effective radius in kpc
* `n`: index
* `Mtot`: total mass
"""
function rho_sersic(r, Re, n, Mtot)
    bn = b_cb(n)
    pn = p_ln(n)
    x = r / Re
    rho0 = Mtot * bn ^ (3 * n) / (4π * n * Re ^ 3 * gamma((3 - pn) * n))
    return rho0 * (@. (bn ^ n * x) ^ (-pn) * exp(-bn * x ^ (1.0 / n)))
end

"""
    M_sersic(r, Re, n, Mtot)

Compute the enclosed mass within deprojected radius for the Sersic model.

* `r`: deprojected (3d) radii in kpc
* `Re`: effective radius in kpc
* `n`: index
* `Mtot`: total mass
"""
function M_sersic(r, Re, n, Mtot)
    pn = p_ln(n)
    bn = b_cb(n)
    alpha = (3 - pn) * n
    x = @. bn * (r / Re) ^ (1.0 / n)
    return @. Mtot * gamma_inc(alpha, x)
end

"""
    SersicModel(Re, n, Mtot)

Representation of a Sersic surface density profile.

* `Re`: the effective radius in kpc.
* `n`: the Sersic index
* `Mtot`: the total mass associated to the profile.

For the purposes of using the `SersicModel` to track the tracer population's
density distribution, `Mtot` is somewhat arbitrary.  For using it as a component
of the total mass profile, this should take units of solar masses.
"""
struct SersicModel <: DensityModel
    Re::Float64
    n::Float64
    Mtot::Float64
end

SersicModel(Re, n) = SersicModel(Re, n, 1.0)
SersicModel() = SersicModel(10.0, 4.0)

function density2d(model::SersicModel, R)
    return s_sersic(R, model.Re, model.n, model.Mtot)
end

function density(model::SersicModel, r)
    return rho_sersic(r, model.Re, model.n, model.Mtot)
end

function mass(model::SersicModel, r)
    return M_sersic(r, model.Re, model.n, model.Mtot)
end

end


text:
#=
learning.jl; Provides the basic functions to do learning.
=#

train(y, hdvs) = Dict(c=>aggregate(hdvs[y.==c]) for c in unique(y))

predict(v::AbstractHDV, centers) = maximum((similarity(v, xcᵢ), yᵢ) for (yᵢ, xcᵢ) in centers)[2]

predict(hdvs::Vector{<:AbstractHDV}, centers) = [maximum((similarity(v, xcᵢ), yᵢ) for (yᵢ, xcᵢ) in centers)[2] for v in hdvs]

function retrain!(centers, y, hdvs; niters=10, verbose=true)
    @assert length(y) == length(hdvs)
    n_obs = length(y)
    wrong = zeros(Bool, n_obs)
    for iter in 1:niters
        verbose && print("Iteration $iter of $niters ...")
        # check all wrongly classified classes
        for (i, (yᵢ, hdv)) in enumerate(zip(y, hdvs))
            (predict(hdv, centers) != yᵢ) && (wrong[i] = true)
        end
        n_errors = sum(wrong)
        verbose && println(" found $n_errors classification errors")
        # aggregate all the mistaken vectors again in the centers
        for (yₖ, cₖ) in centers
            aggregatewith!(cₖ, hdvs[(y.==yₖ) .& wrong])
        end
        fill!(wrong, false)
    end
    return centers
end


text:
using Pkg
Pkg.activate(".")
using ConceptnetNumberbatch
using Serialization

local cptnet

# Load a serialized version of ConceptNetEnglish
fid = open("./_conceptnet_/numberbatch-en-17.06.txt.bin")
cptnet = deserialize(fid)
close(fid)

## Get an embedding matrix
phrase = "this is a phrase that containz some iwords"
ConceptnetNumberbatch.word_embeddings(cptnet, phrase, keep_size=false, search_mismatches=false)
@time embs=ConceptnetNumberbatch.word_embeddings(cptnet, phrase, keep_size=false, search_mismatches=false)
@time embs=ConceptnetNumberbatch.word_embeddings(cptnet, phrase, keep_size=false, search_mismatches=true)
println("Loaded $(size(embs, 2)) embedding vectors (out of $(length(split(phrase)))), $(size(embs,1)) elements each.")


text:
module TestCusp
include("preamble_plots.jl")

using Bifurcations: resolved_points

solver1 = init(Bifurcations.Examples.Cusp.prob)
solve!(solver1)

codim1_points = resolved_points(solver1)
@test length(codim1_points) == 2
@test codim1_points[1].point_type === Codim1.PointTypes.saddle_node
@test codim1_points[2].point_type === Codim1.PointTypes.saddle_node

sn_prob = BifurcationProblem(
    codim1_points[1],
    solver1,
    (@lens _.β₂),
    (-2.0, 2.0),
)
sn_solver = init(
    sn_prob;
)
solve!(sn_solver)

β₁s = [u[end-1] for sweep in sn_solver.super.sol.sweeps for u in sweep.u]
β₂s = [u[end]   for sweep in sn_solver.super.sol.sweeps for u in sweep.u]
@test all(@. isapprox(β₁s^2, (2 / (3 * √3) *  β₂s^(3 / 2))^2; atol=1e-6))
@test maximum(β₁s) > 1e-6
@test minimum(β₁s) < -1e-6
@test maximum(β₂s) > 2
@test minimum(β₂s) > -1e-6

codim2_points = resolved_points(sn_solver)
@test length(codim2_points) == 1
@test codim2_points[1].point_type === Codim2.PointTypes.cusp
β_cusp = codim2_points[1].u[end-1:end]
@test all(@. abs(β_cusp) < 1e-6)

plt = plot(sn_solver)
@test_nothrow nullshow(plt)

end  # module


text:
mutable struct Scene
    tNow::Float64
    dt::Float64
    lastTimeCalled::Float64
    paused::Bool
    contents::Array
end
Scene(dt, contents::Array) = Scene(0,dt,time(),false,contents)
Scene(contents::Array) = Scene(0.05, contents)

function evolve_tick(scene::Scene)
    scene.tNow += scene.dt
    for c in scene.contents
        evolve_tick(c)
    end
end

function evolve_tick(Any); end; # no op


text:
"""
    Problems
        Collection of trajectory optimization problems
"""
module Problems

using Altro
using TrajectoryOptimization  # import the API
using LinearAlgebra
using StaticArrays
using RobotZoo
using RobotDynamics
using Rotations
const Dynamics = RobotDynamics
const TO = TrajectoryOptimization

include("../problems/doubleintegrator.jl")
include("../problems/pendulum.jl")
include("../problems/acrobot.jl")
include("../problems/parallel_park.jl")
include("../problems/cartpole.jl")
include("../problems/dubins_car.jl")
include("../problems/bicycle.jl")
include("../problems/quadrotor.jl")
include("../problems/airplane.jl")

export
    doubleintegrator,
    pendulum,
    parallel_park,
    cartpole,
    doublependulum,
    acrobot,
    car_escape,
    car_3obs,
    quadrotor,
    quadrotor_maze,
    kuka_obstacles

export
    plot_escape,
    plot_car_3obj,
    quadrotor_maze_objects,
    kuka_obstacles_objects

end


text:
Base.@kwdef mutable struct ESARSA{O, T<:AbstractTraceUpdate} <: LearningUpdate
    lambda::Float64
    opt::O
    # trace::T = ReplacingTraces()
    trace::T
    e::IdDict = IdDict()
    prev_discounts::IdDict = IdDict()
end

function get_demon_parameters(lu::ESARSA, learner, demons, obs, state, action, next_obs, next_state, next_action, env_reward)
    C, next_discounts, _ = get(demons; state_t = obs, action_t = action, state_tp1 = next_obs, action_tp1 = next_action, reward = env_reward)
    target_pis = get_demon_pis(demons, learner.num_actions, state, obs)
    next_target_pis = get_demon_pis(demons, learner.num_actions, next_state, next_obs)
    C, next_discounts, target_pis, next_target_pis
end

function update!(lu::ESARSA,
                 learner::QLearner{M, LU},
                 demons,
                 obs,
                 next_obs,
                 state,
                 action,
                 next_state,
                 next_action,
                 is_terminal,
                 behaviour_pi_func,
                 reward) where {M<:AbstractMatrix, LU<:ESARSA}

    weights = learner.model
    λ = lu.lambda
    e = get!(()->zero(weights), lu.e, weights)::typeof(weights)

    # next_target_pis = behaviour_pi_func(next_state, next_obs)
    C, next_discounts, target_pis, next_target_pis = get_demon_parameters(lu, learner, demons, obs, state, action, next_obs, next_state, next_action, reward)
    b_πs = behaviour_pi_func(state, obs)

    discounts = get!(()->zero(next_discounts), lu.prev_discounts, learner)::typeof(next_discounts)

    #NOTE: With exploring starts an action could be taken that goes against b_π with prob 0.
    ρ = if b_πs[action] == 0
        zeros(length(target_pis[:,action]))
    else
        target_pis[:,action] ./ b_πs[action]
    end

    inds = get_action_inds(action, learner.num_actions, learner.num_demons)
    state_action_row_ind = inds

    update_trace!(lu.trace,
                  e,
                  state,
                  λ,
                  repeat(discounts, inner = learner.num_actions),
                  repeat(ρ, inner = learner.num_actions),
                  inds)

    next_preds = learner(next_state)
    pred = learner(state, action)

    Qs = reshape(next_preds, (learner.num_actions, learner.num_demons))'
    td_err = C .+ next_discounts .* sum(next_target_pis .* Qs, dims = 2) - pred
    td_err_across_demons = repeat(vec(td_err), inner=learner.num_actions)

    if lu.opt isa Auto
        next_state_action_row_ind = get_action_inds(next_action, learner.num_actions, learner.num_demons)
        state_discount = zero(e)
        state_discount[state_action_row_ind,:] .+= state'
        state_discount[next_state_action_row_ind,:] .-= next_discounts * next_state'
        abs_phi = abs.(e)
        update!(lu.opt, weights, e, td_err_across_demons, abs_phi .* max.(state_discount, abs_phi), learner.num_demons, learner.num_actions)
    else
        Flux.Optimise.update!(lu.opt, weights,  -(e .* td_err_across_demons))
    end
    discounts .= next_discounts
end

function update!(lu::ESARSA,
                 learner::Union{SRLearner,GPI},
                 demons,
                 obs,
                 next_obs,
                 state::SparseVector,
                 action,
                 next_state,
                 next_action,
                 is_terminal,
                 behaviour_pi_func,
                 env_reward)

    ψ = learner.ψ
    w = learner.r_w

    C, next_discounts, target_pis, next_target_pis =
        get_demon_parameters(lu,
                             learner,
                             demons,
                             obs,
                             state,
                             action,
                             next_obs,
                             next_state,
                             next_action,
                             env_reward)

    discounts = get!(()->zero(next_discounts), lu.prev_discounts, learner)::typeof(next_discounts)
    e_nz = get!(()->Int[], lu.e, learner)::Array{Int, 1}
    e_w_nz = get!(()->Int[], lu.e, learner)::Array{Int, 1}
    e_ψ = get!(()->zero(ψ), lu.e, ψ)::typeof(ψ)
    e_w = get!(()->zero(w), lu.e, w)::typeof(w)
    λ = lu.lambda

    next_active_state_action = get_active_action_state_vector(next_state, next_action,length(next_state), learner.num_actions)
    active_state_action = get_active_action_state_vector(state, action, length(state), learner.num_actions)

    projected_state = learner.feature_projector(obs, action, next_obs)

    (reward_C, SF_C) = C[1:learner.num_tasks] , C[learner.num_tasks + 1:end]
    (reward_discounts, SF_discounts) = discounts[1:learner.num_tasks], discounts[learner.num_tasks+1:end]
    (reward_next_discounts, SF_next_discounts) = next_discounts[1:learner.num_tasks], next_discounts[learner.num_tasks+1:end]
    (reward_target_pis, SF_target_pis) = target_pis[1:learner.num_tasks,:], target_pis[learner.num_tasks+1:end,:]
    (reward_next_target_pis, SF_next_target_pis) = next_target_pis[1:learner.num_tasks,:], next_target_pis[learner.num_tasks+1:end, :]

    b_πs = behaviour_pi_func(state, obs)
    #NOTE: With exploring starts an action could be taken that goes against b_π with prob 0.
    ρ_SF = if b_πs[action] == 0
        zeros(length(SF_target_pis[:,action]))
    else
        SF_target_pis[:,action] ./ b_πs[action]
    end
    ρ_reward = if b_πs[action] == 0
        zeros(length(reward_target_pis[:,action]))
    else
        reward_target_pis[:,action] ./ b_πs[action]
    end

    # Update Traces: See update_utils.jl
    update_trace!(lu.trace, e_ψ, active_state_action, λ, SF_discounts, ρ_SF)
    update_trace!(lu.trace, e_w, projected_state, λ, reward_discounts, ρ_reward)
    e_nz = e_nz ∪ active_state_action.nzind
    e_w_nz = e_w_nz ∪ projected_state.nzind

    pred = ψ * next_active_state_action
    reward_feature_backup = zeros(length(SF_C))
    for a in 1:learner.num_actions
        next_possible_active_state_action = get_active_action_state_vector(next_state, a, length(next_state), learner.num_actions)
        reward_feature_backup += SF_next_target_pis[:,a] .* (ψ * next_possible_active_state_action)
    end

    target = SF_C + SF_next_discounts .* reward_feature_backup
    td_err = target - ψ * active_state_action

    pred_err = reward_C - w * projected_state

    # This should always be true as this is immediate next step prediction which is equivalent to having discounts of 0 for all states
    @assert sum(reward_discounts) == 0
    #td_err is (336x1)
    # TD err is applied across rows

    if lu.opt isa Auto
        # next_state_action_row_ind = get_action_inds(next_action, learner.num_actions, learner.num_demons)
        state_discount = -SF_next_discounts * next_active_state_action'
        state_discount .+= active_state_action'
        abs_ϕ_ψ = if λ == 0.0
            abs.(repeat(active_state_action, outer=(1, length(td_err)))')
        else
            abs.(e_ψ)
        end
        z = abs_ϕ_ψ .* max.(abs_ϕ_ψ, state_discount)
        update!(lu.opt, ψ, e_ψ, td_err, z,  learner.num_demons - learner.num_tasks, 1)

        # Reward next_discounts is always 0, so don't need projected_next_state_action
        # state_discount_r = -reward_next_discounts * projected_next_state_action'
        state_discount_r = projected_state'
        abs_ϕ_w = if λ == 0.0
            abs.(repeat(projected_state, outer=(1, length(pred_err)))')
        else
            abs.(e_w)
        end
        z_r = abs_ϕ_w .* max.(abs_ϕ_w, state_discount_r)
        update!(lu.opt, w, e_w, pred_err, z_r, learner.num_tasks, 1)
    elseif lu.opt isa Flux.Descent
        α = lu.opt.eta
        if λ == 0.0
            ψ[:, active_state_action.nzind] .+= (α  * td_err) * active_state_action.nzval'
            w .= w .+ α * pred_err * projected_state'
        else
            ψ[:, e_nz] .+= (α  * td_err) .* e_ψ[:, e_nz]
            w[:, e_w_nz] .+= (α  * pred_err) .* e_w[:, e_w_nz]
        end
    else
        update!(lu.opt, ψ, - e_ψ .* td_err)
        update!(lu.opt, w, - e_w .* pred_err)
    end
    discounts .= next_discounts
end

function zero_eligibility_traces!(lu::ESARSA)
    for (k, v) ∈ lu.e
        if eltype(v) <: Integer
            lu.e[k] = Int[]
        else
            v .= 0
        end
    end
end

max_stars_repo_name: MatthewMcLeod/curiosity

max_issues_repo_name: MatthewMcLeod/curiosity

max_forks_repo_name: MatthewMcLeod/curiosity


text:
using LinearAlgebra
using PolynomialRoots
using Statistics

include("utils.jl")

struct Planar3RPR
    ρₘᵢₙ::Vector{Real}
    e::Vector{Real}
    f::Vector{Real}
    T::Matrix{Real}
    Θₑ::Vector{Real}
    Θ::Vector{Real}
    Fx
    Fy
    function Planar3RPR(ρₘᵢₙ, e, f, rot=0, P=[0, 0])
        @assert length(ρₘᵢₙ)==3 "Provide length of 3 Prismatic Links"
        @assert length(e)==3 "Provide length of 3 sides of End-effector triangle"
        @assert length(f)==3 "Provide length of 3 sides of Fixed/Base triangle"
        @assert ρₘᵢₙ.+√3e/2 > √3f/2 "ρₘᵢₙ + √3 E/2 should be > √3 F/2"
        T = [cosd(rot) -sind(rot) P[1];
             sind(rot)  cosd(rot) P[2];
                     0          0    1]
        Ae = acosd((e[1]^2 + e[3]^2 - e[2]^2) / (2 * e[1] * e[3]))
        Be = acosd((e[1]^2 + e[2]^2 - e[3]^2) / (2 * e[1] * e[2]))
        Ce = acosd((e[2]^2 + e[3]^2 - e[1]^2) / (2 * e[2] * e[3]))
        Θₑ = [Ae, Be, Ce]
        A = acosd((f[1]^2 + f[3]^2 - f[2]^2) / (2 * f[1] * f[3]))
        B = acosd((f[1]^2 + f[2]^2 - f[3]^2) / (2 * f[1] * f[2]))
        C = acosd((f[2]^2 + f[3]^2 - f[1]^2) / (2 * f[2] * f[3]))
        Θ = [A, B, C]
        F = T * [0.0 f[1] f[3] * cosd(Ae);
                 0.0  0.0 f[3] * sind(Ae);
                 1.0  1.0             1.0]
        Fx = F[1, :]
        Fy = F[2, :]
        return new(a, p, e, f, T, Θₑ, Θ, Fx, Fy)
    end
end

function Planar3RPR(ρₘᵢₙ::N, e::N, f::N, rot=0, P=[0, 0]) where N<:Real
    return Planar3RPR([ρₘᵢₙ, ρₘᵢₙ, ρₘᵢₙ], [e, e, e], [f, f, f], rot, P)
end


text:
module BoxLeastSquares

using Statistics

export BLS, autoperiod

# core BLS algorithm and struct
include("bls.jl")
# transit model
include("model.jl")
# pretty printing
include("printing.jl")
# plotting recipes
include("plotting.jl")

end


text:
function inverse(side::HorizonSide)
    return HorizonSide(mod(Int(side) + 2, 4))
end

function f(x::HorizonSide, y::HorizonSide)
    steps = 0
    while !isborder(robot, x) && !isborder(robot, y)
        steps += 1
        move!(robot, x)
        move!(robot, y)
        putmarker!(robot)
    end
    x = inverse(x)
    y = inverse(y)
    while steps > 0
        move!(robot, x)
        move!(robot, y)
        steps -= 1
    end
end

putmarker!(robot)
f(Nord, West)
f(Sud, Ost)
f(Ost, Nord)
f(West, Sud)

max_stars_repo_name: TheXlebushek/mirea-progs

max_issues_repo_name: TheXlebushek/mirea-progs

max_forks_repo_name: TheXlebushek/mirea-progs


text:
import CombinedParsers: ParserTypes, tryparsenext
using Nullables
import Nullables: isnull, Nullable

struct Greedy{Ps,A} <: CombinedParsers.AbstractToken{Any}
    pairs::Ps
    alt::A
end

export greedy
function greedy(tokens...; alt = [])
    Greedy([tokens...], alt)
end

function CombinedParsers.tryparsenext(tokf::Greedy, str, i, till, opts=TextParse.default_opts)
    T = result_type(tokf)
    sections=tokf.pairs
    RT(key, value) = if value[2] isa ParserTypes
        if Missing <: result_type(typeof(key))
            result_type(typeof(value[2]))
        else
            promote_type(result_type(typeof(key)), result_type(typeof(value[2])))
        end
    else
        result_type(typeof(key))
    end
    R = Dict([ value[1] => Vector{RT(key,value)}() for (key,value) in sections]...,
             [ key => Vector{result_type(typeof(value))}() for (key,value) in tokf.alt]...
             )
    hist = Any[]
    last_section = nothing
    last_content = nothing
    aggregator = :head
    function first_match(j)
        local repval, i__
        for (key, content) in sections
            repval, i__ = tryparsenext(key, str, j, till)
            !isnull(repval) && return key, content, repval, i__
        end
        return (first(sections)..., repval, i__)
    end
    head = nothing
    i_ = i ##isnull(1)
    while true
        key, content, r, i__ = first_match(i_)
        save = if isnull(r)
            cr, ci = if last_content === nothing || last_content === missing
                Nullable{T}(), i
            else
                tryparsenext(last_content, str, i_, till)
            end
            ai = 0
            while ai < lastindex(tokf.alt) && isnull(cr)
                ai = ai+1
                cr, ci = tryparsenext(tokf.alt[ai].second, str, i_, till)
            end
            if isnull(cr)
                return Nullable{T}(convert(T, R)), i_
            elseif ai == 0
                push!(hist, get(cr))
                i__ = ci
                false
            else
                aggregator != :head && append!(R[aggregator],hist)
                hist = [get(cr)]
                (aggregator, last_content) = tokf.alt[ai]
                last_section = ai
            end
        else
            if last_section !== nothing
                append!(R[aggregator],hist)
            end
            hist = get(r) !== missing ? [get(r)] : Vector{RT(key,content)}()
            aggregator, last_content = content
            last_section = key
        end
        i_ = i__
    end
    error("unreachable")
end


text:
abstract type AbstractPlanetEnvironment end

struct PlanetEnvironment <: AbstractPlanetEnvironment
    GravityModel::Function
    PlanetModel::Function
    AtmosphereFlag::Bool #Boolean to determine whether atmoshpere model should be looked at
    AtmosphereModel::Function
    WindFlag::Bool #Boolean to determine whether wind model should be looked at
    WindModel::Function
end


text:
import SpecialFunctions: gamma, loggamma, erf, erfc, erfcinv, erfi, erfcx,
                         dawson, digamma, trigamma, invdigamma, polygamma,
                         airyai, airyaiprime, airybi, airybiprime, besselj0,
                         besselj1, bessely0, bessely1, besselj, bessely, besseli,
                         besselk, hankelh1, hankelh2, polygamma, beta, logbeta

const monadic = [deg2rad, rad2deg, transpose, conj, asind, log1p, acsch,
                 acos, asec, acosh, acsc, cscd, log, tand, log10, csch, asinh,
                 abs2, cosh, sin, cos, atan, cospi, cbrt, acosd, acoth, acotd,
                 asecd, exp, acot, sqrt, sind, sinpi, asech, log2, tan, exp10,
                 sech, coth, asin, cotd, cosd, sinh, abs, csc, tanh, secd,
                 atand, sec, acscd, cot, exp2, expm1, atanh, real, gamma,
                 loggamma, erf, erfc, erfcinv, erfi, erfcx, dawson, digamma,
                 trigamma, invdigamma, polygamma, airyai, airyaiprime, airybi,
                 airybiprime, besselj0, besselj1, bessely0, bessely1]

const diadic = [max, min, hypot, atan, mod, rem, copysign,
                besselj, bessely, besseli, besselk, hankelh1, hankelh2,
                polygamma, beta, logbeta]
const previously_declared_for = Set([])

# TODO: it's not possible to dispatch on the symtype! (only problem is Parameter{})

assert_like(f, T) = nothing
function assert_like(f, T, a, b...)
    islike(a, T) || throw(ArgumentError("The function $f cannot be applied to $a which is not a $T-like object." *
                                        "Define `islike(::$(typeof(a)), ::Type{$T}) = true` to enable this."))
    assert_like(f, T, b...)
end

islike(a, T) = symtype(a) <: T

# TODO: keep domains tighter than this
function number_methods(T, rhs1, rhs2, options=nothing)
    exprs = []

    skip_basics = options !== nothing ? options == :skipbasics : false
    skips = Meta.isexpr(options, [:vcat, :hcat, :vect]) ? Set(options.args) : []
    basic_monadic = [-, +]
    basic_diadic = [+, -, *, /, \, ^]

    rhs2 = :($assert_like(f, Number, a, b); $rhs2)
    rhs1 = :($assert_like(f, Number, a); $rhs1)

    for f in (skip_basics ? diadic : vcat(basic_diadic, diadic))
        nameof(f) in skips && continue
        for S in previously_declared_for
            push!(exprs, quote
                      (f::$(typeof(f)))(a::$T, b::$S) = $rhs2
                      (f::$(typeof(f)))(a::$S, b::$T) = $rhs2
                  end)
        end

        # TODO: modularize and make another macro?
        expr = quote
            (f::$(typeof(f)))(a::$T, b::$T) = $rhs2
            (f::$(typeof(f)))(a::$T, b::Real)   = $rhs2
            (f::$(typeof(f)))(a::Real, b::$T)   = $rhs2
            (f::$(typeof(f)))(a::$T, b::Number) = $rhs2
            (f::$(typeof(f)))(a::Number, b::$T) = $rhs2
        end

        push!(exprs, expr)
    end

    for f in (skip_basics ? monadic : vcat(basic_monadic, monadic))
        nameof(f) in skips && continue
        push!(exprs, :((f::$(typeof(f)))(a::$T)   = $rhs1))
    end
    push!(exprs, :(push!($previously_declared_for, $T)))
    Expr(:block, exprs...)
end

macro number_methods(T, rhs1, rhs2, options=nothing)
    number_methods(T, rhs1, rhs2, options) |> esc
end

@number_methods(Sym, term(f, a), term(f, a, b), skipbasics)
@number_methods(Term, term(f, a), term(f, a, b), skipbasics)
@number_methods(Add, term(f, a), term(f, a, b), skipbasics)
@number_methods(Mul, term(f, a), term(f, a, b), skipbasics)
@number_methods(Pow, term(f, a), term(f, a, b), skipbasics)

for f in diadic
    @eval promote_symtype(::$(typeof(f)),
                   T::Type{<:Number},
                   S::Type{<:Number}) = promote_type(T, S)
end

for f in [+, -, *, \, /, ^]
    @eval promote_symtype(::$(typeof(f)),
                   T::Type{<:Number},
                   S::Type{<:Number}) = promote_type(T, S)
end

promote_symtype(::typeof(rem2pi), T::Type{<:Number}, mode) = T
Base.rem2pi(x::Symbolic, mode::Base.RoundingMode) = term(rem2pi, x, mode)

for f in monadic
    if f in [real]
        continue
    end
    @eval promote_symtype(::$(typeof(f)), T::Type{<:Number}) = promote_type(T, Real)
    @eval (::$(typeof(f)))(a::Symbolic)   = term($f, a)
end

Base.:*(a::AbstractArray, b::Symbolic{<:Number}) = map(x->x*b, a)
Base.:*(a::Symbolic{<:Number}, b::AbstractArray) = map(x->a*x, b)

for f in [identity, one, zero, *, +, -]
    @eval promote_symtype(::$(typeof(f)), T::Type{<:Number}) = T
end

promote_symtype(::typeof(Base.real), T::Type{<:Number}) = Real
Base.real(s::Symbolic) = islike(s, Real) ? s : term(real, s)

## Booleans

# binary ops that return Bool
for (f, Domain) in [(==) => Number, (!=) => Number,
                    (<=) => Real,   (>=) => Real,
                    (isless) => Real,
                    (<) => Real,   (> ) => Real,
                    (& ) => Bool,   (| ) => Bool,
                    xor => Bool]
    @eval begin
        promote_symtype(::$(typeof(f)), ::Type{<:$Domain}, ::Type{<:$Domain}) = Bool
        (::$(typeof(f)))(a::Symbolic{<:$Domain}, b::$Domain) = term($f, a, b, type=Bool)
        (::$(typeof(f)))(a::Symbolic{<:$Domain}, b::Symbolic{<:$Domain}) = term($f, a, b, type=Bool)
        (::$(typeof(f)))(a::$Domain, b::Symbolic{<:$Domain}) = term($f, a, b, type=Bool)
    end
end

for f in [!, ~]
    @eval begin
        promote_symtype(::$(typeof(f)), ::Type{<:Bool}) = Bool
        (::$(typeof(f)))(s::Symbolic{Bool}) = Term{Bool}(!, [s])
    end
end

# An ifelse node, ifelse is a built-in unfortunately
# So this uses IfElse.jl's ifelse that we imported
function ifelse(_if::Symbolic{Bool}, _then, _else)
    Term{Union{symtype(_then), symtype(_else)}}(ifelse, Any[_if, _then, _else])
end
promote_symtype(::typeof(ifelse), _, ::Type{T}, ::Type{S}) where {T,S} = Union{T, S}
Base.@deprecate cond(_if, _then, _else) ifelse(_if, _then, _else)

# Specially handle inv and literal pow
Base.inv(x::Symbolic) = Base.:^(x, -1)
Base.literal_pow(::typeof(^), x::Symbolic, ::Val{p}) where {p} = Base.:^(x, p)


text:
using Documenter
using MartaCT

makedocs(
    sitename = "MartaCT",
    format = Documenter.HTML(
        prettyurls = "local" ∉ ARGS,
    ),
    modules = [MartaCT],
    pages = [
        "Home" => "index.md",
        "Getting Started" => [
            "Installation.md",
        ],
        "Manual" => [
            "Tutorial.md",
            "Quantum Optics" => "QuantumOptics.md",
        ],
        "API Reference" => [
            "api/Algorithms.md",
            "api/CTImages.md",
            "api/Geometry.md",
            "api/CTIO.md",
            "api/CTPlots.md",
            "api/Simulations.md",
            "api/Calibration.md",
            "api/TestImages.md",
            "api/CTScan.md",
            "api/Info.md",
            "api/Interpolation.md",
            "api/Applicative.md",
            "api/Monads.md",
            "api/Utils.md",
        ],
        "Index" => "GeneralIndex.md",
    ]
)


text:
@symbol_func function K_P(cur_reactor::AbstractReactor)
  cur_K = 5.0 + cur_reactor.Q

  cur_K /= 5.0 * cur_reactor.Q

  cur_K *= K_F(cur_reactor)

  cur_K
end


text:
"""
A [logarithmic number system](https://en.wikipedia.org/wiki/Logarithmic_number_system).

Provides the signed [`Logarithmic`](@ref) and unsigned [`ULogarithmic`](@ref) types, which represent real numbers and positive real numbers respectively.
"""
module LogarithmicNumbers

using Requires, Random

import Base:
	exp, log, *, /, ^, inv, +, -, prod, sum, angle, show,
	write, read, float, big, unsigned, signed, widen,
	typemin, typemax, zero, one, iszero, isone, isinf,
	isfinite, isnan, sign, signbit, abs, ==, isequal, <, ≤,
	cmp, isless, nextfloat, prevfloat, rand, promote_rule,
	conj, real, imag

export ULogarithmic, Logarithmic, CLogarithmic

include("types.jl")
include("constructors.jl")
include("ulogarithmic.jl")
include("logarithmic.jl")
include("clogarithmic.jl")
include("promotion.jl")
include("arithmetic.jl")
include("random.jl")
include("io.jl")
include("init.jl")

# type aliases
for T in (Float64, Float32, Float16, BigFloat)
  ULogT = Symbol(:ULog, T)
  LogT = Symbol(:Log, T)
  CLogT = Symbol(:CLog, T)
  @eval begin
      const $(ULogT) = $(ULogarithmic{T})
      const $(LogT) = $(Logarithmic{T})
      const $(CLogT) = $(CLogarithmic{T})
      export $(ULogT), $(LogT), $(CLogT)
  end
end

end # module


text:
using Pkg
Pkg.activate(".")

Pkg.add(name = "PlanktonIndividuals", version = "0.1.9")
Pkg.add(name = "Oceananigans", version = "0.55.0")
pkg"add Random, Statistics, Printf, YAML, Serialization, SeawaterPolynomials"

using Random, Statistics, Printf, YAML, Serialization
using Oceananigans
using Oceananigans.Units: day
using SeawaterPolynomials
using PlanktonIndividuals

#params_path = "params_autotroph.yml" # parameters for autotrophs
params_path = "params_mixotroph.yml" # parameters for mixotrophs
paths_path = "paths_nut_init.yml"
par_path = "../Input/par_doc.bin"
nday = 30*12
println(params_path)

Ogrid = RegularRectilinearGrid(size = (25, 1, 50), extent = (100, 4, 200),topology = (Periodic, Periodic, Bounded),)

buoyancy = SeawaterBuoyancy(
    equation_of_state = SeawaterPolynomials.TEOS10EquationOfState(),
    constant_salinity = 35.0 # psu
)

N² = 1e-5 # s⁻²
α = SeawaterPolynomials.thermal_expansion(20, 35, 0, buoyancy.equation_of_state)
g = buoyancy.gravitational_acceleration
ρᵣ = buoyancy.equation_of_state.reference_density
∂T∂z = ρᵣ * N² / (α * g)
bottom_temperature_boundary_condition = GradientBoundaryCondition(∂T∂z)

peak_outgoing_radiation = 200 # Watts / m²
heat_capacity = 3991 # J / kg / ᵒC
reference_density = buoyancy.equation_of_state.reference_density # kg m⁻³
peak_outgoing_flux = peak_outgoing_radiation / (reference_density * heat_capacity)
Qᵇ = α * g * peak_outgoing_flux

@inline diurnal_cycle(t, day) = max(0, - cos(2π * t / day))
@inline nocturnal_cycle(t, day) = max(0, cos(2π * t / day))
@inline outgoing_flux(x, y, t, p) = p.peak * nocturnal_cycle(t, p.day)

surface_temperature_boundary_condition = FluxBoundaryCondition(outgoing_flux, parameters=(day=day, peak=peak_outgoing_flux))

T_bcs = TracerBoundaryConditions(Ogrid, bottom = bottom_temperature_boundary_condition,
                                 top = surface_temperature_boundary_condition)

light_attenuation_scale = 20 # m
surface_solar_insolation = 400 # Watts / m²
surface_solar_temperature_flux = surface_solar_insolation / (reference_density * heat_capacity)

@inline daylight(z, t, λ, day) = exp(z / λ) * diurnal_cycle(t, day)
@inline solar_flux_divergence(z, t, Qᴵ, λ, day) = Qᴵ / λ * daylight(z, t, λ, day)
@inline diurnal_solar_flux_divergence(x, y, z, t, p) =
    max(0, solar_flux_divergence(z, t, p.surface_flux, p.attenuation, p.day))

interior_heating = Forcing(diurnal_solar_flux_divergence,
                                 parameters = (surface_flux = surface_solar_temperature_flux,
                                               attenuation = light_attenuation_scale,
                                               day = day))

Omodel = IncompressibleModel(architecture = CPU(),
                                    grid = Ogrid,
                                 closure = AnisotropicMinimumDissipation(),
                                coriolis = FPlane(f=-1e-4),
                                 tracers = :T, #tracer_names,
                                buoyancy = buoyancy,
                     boundary_conditions = (T=T_bcs,),
                                 forcing = (T=interior_heating,),
)

initial_temperature(x, y, z) = (20 + ∂T∂z * z + ∂T∂z * Ogrid.Lz * 1e-4 * randn() * exp(z / (8 * Ogrid.Δz)))

set!(Omodel, T = initial_temperature,)

wizard = TimeStepWizard(cfl=0.2, Δt=1.0, max_change=1.1, max_Δt=20.0)

simulation = Simulation(Omodel, iteration_interval = 1, stop_time=86400, Δt = wizard)

run!(simulation)

resultpath = PrepRunDir()
phy_grid = read_Ogrids(Omodel.grid, resultpath);
RunParam.nTime = 1440*nday
RunParam.ΔT = 60

# update parameters
phy_params = YAML.load(open(params_path))
update_params!(RunParam.params,phy_params)

paths = YAML.load(open(paths_path))
phy_model = PI_Model(phy_grid, RunParam,
                     PAR = read_IR_input(RunParam.ΔT, phy_grid, par_path),
                     nutrients = load_nut_initials(paths, phy_grid));

Nsimulation = Simulation(Omodel, Δt=20.0, stop_time=86400, iteration_interval=1)

for i in 1:RunParam.nTime
    vel_field =[]
    for j in 1:3
        Nsimulation.stop_time += 20
        run!(Nsimulation)
        u = Omodel.velocities.u.data.parent
        v = Omodel.velocities.v.data.parent
        w = Omodel.velocities.w.data.parent
        vel = PlanktonIndividuals.velocity(u, v, w)
        push!(vel_field,vel)
    end
    vel_itps = (generate_vel_itp(phy_model.grid, vel_field[1]),
                generate_vel_itp(phy_model.grid, vel_field[2]),
                generate_vel_itp(phy_model.grid, vel_field[3]))
    PI_advectRK4!(phy_model, RunParam.ΔT, vel_itps)
    PI_TimeStep!(phy_model, RunParam.ΔT, vel_field[end], resultpath)

    phyts_sp = sort_species(phy_model.individuals.phytos, phy_model.params["P_Nsp"])
    write_species_dynamics(phy_model.t, phyts_sp, resultpath)

end

open(resultpath*"diags_spcs.bin", "w") do io
    serialize(io, phy_model.diags.spcs)
end
open(resultpath*"diags_pop.bin", "w") do io
    serialize(io, phy_model.diags.pop)
end
open(resultpath*"diags_tr.bin", "w") do io
    serialize(io, phy_model.diags.tr)
end


text:
using FermiCG
using Printf
using Test

@testset "Clusters" begin
    atoms = []
    push!(atoms,Atom(1,"H",[0,0,0]))
    push!(atoms,Atom(2,"H",[0,0,1]))
    push!(atoms,Atom(3,"H",[0,0,2]))
    push!(atoms,Atom(4,"H",[0,0,3]))
    push!(atoms,Atom(5,"H",[0,0,4]))
    push!(atoms,Atom(6,"H",[0,0,5]))
    basis = "6-31g"
    basis = "sto-3g"
    mol     = Molecule(0,1,atoms,basis)

    mf = FermiCG.pyscf_do_scf(mol)
    nbas = size(mf.mo_coeff)[1]
    ints = FermiCG.pyscf_build_ints(mol,mf.mo_coeff, zeros(nbas,nbas));
    e_fci, d1_fci, d2_fci = FermiCG.pyscf_fci(ints,1,1)
    @printf(" FCI Energy: %12.8f\n", e_fci)

    C = mf.mo_coeff
    Cl = FermiCG.localize(mf.mo_coeff,"lowdin",mf)
    FermiCG.pyscf_write_molden(mol,Cl,filename="lowdin.molden")
    S = FermiCG.get_ovlp(mf)
    U =  C' * S * Cl
    println(" Build Integrals")
    flush(stdout)
    ints = FermiCG.orbital_rotation(ints,U)
    println(" done.")
    flush(stdout)

    clusters    = [(1:2),(3:4),(5:6)]
    #clusters    = [(1:4),(5:8),(9:12)]
    init_fspace = [(1,1),(1,1),(1,1)]

    max_roots = 20

    clusters = [Cluster(i,collect(clusters[i])) for i = 1:length(clusters)]
    display(clusters)

    cluster_bases = FermiCG.compute_cluster_eigenbasis(ints, clusters, verbose=1)

    #
    #   Just test that the sum of all the basis vector matrices is reproduced
    println("")
    tst1 = 0
    for ci in clusters
        display(cluster_bases[ci.idx])
        for (sector,vecs) in cluster_bases[ci.idx].basis
            tst1 += sum(abs.(vecs))
        end
    end
    println(tst1)
    @test isapprox(tst1, 66.05063700792823, atol=1e-10)

    # now try with restrictions on fock space, and dimensions
    cluster_bases = Vector{ClusterBasis}()
    max_roots=2
    cluster_bases = FermiCG.compute_cluster_eigenbasis(ints, clusters, verbose=1,
                                                       max_roots=2, init_fspace=init_fspace, delta_elec=1)

    println("")
    tst1 = 0
    for ci in clusters
        display(cluster_bases[ci.idx])
        for (sector,vecs) in cluster_bases[ci.idx].basis
            tst1 += sum(abs.(vecs))
        end
    end
    println(tst1)
    @test isapprox(tst1, 43.466561233596934, atol=1e-10)

end

max_stars_repo_name: nmayhall-vt/FermiCG

max_issues_repo_name: nmayhall-vt/FermiCG

max_forks_repo_name: nmayhall-vt/FermiCG


text:
#Julia Control Flow

function larger(x, y)
    if (x>y)
        return x
    end
    return y
end
println(larger(7,8))


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct VpnClientParameters <: SwaggerModel
    processorArchitecture::Any # spec type: Union{ Nothing, String } # spec name: processorArchitecture
    authenticationMethod::Any # spec type: Union{ Nothing, AuthenticationMethod } # spec name: authenticationMethod
    radiusServerAuthCertificate::Any # spec type: Union{ Nothing, String } # spec name: radiusServerAuthCertificate
    clientRootCertificates::Any # spec type: Union{ Nothing, Vector{String} } # spec name: clientRootCertificates

    function VpnClientParameters(;processorArchitecture=nothing, authenticationMethod=nothing, radiusServerAuthCertificate=nothing, clientRootCertificates=nothing)
        o = new()
        validate_property(VpnClientParameters, Symbol("processorArchitecture"), processorArchitecture)
        setfield!(o, Symbol("processorArchitecture"), processorArchitecture)
        validate_property(VpnClientParameters, Symbol("authenticationMethod"), authenticationMethod)
        setfield!(o, Symbol("authenticationMethod"), authenticationMethod)
        validate_property(VpnClientParameters, Symbol("radiusServerAuthCertificate"), radiusServerAuthCertificate)
        setfield!(o, Symbol("radiusServerAuthCertificate"), radiusServerAuthCertificate)
        validate_property(VpnClientParameters, Symbol("clientRootCertificates"), clientRootCertificates)
        setfield!(o, Symbol("clientRootCertificates"), clientRootCertificates)
        o
    end
end # type VpnClientParameters

const _property_map_VpnClientParameters = Dict{Symbol,Symbol}(Symbol("processorArchitecture")=>Symbol("processorArchitecture"), Symbol("authenticationMethod")=>Symbol("authenticationMethod"), Symbol("radiusServerAuthCertificate")=>Symbol("radiusServerAuthCertificate"), Symbol("clientRootCertificates")=>Symbol("clientRootCertificates"))
const _property_types_VpnClientParameters = Dict{Symbol,String}(Symbol("processorArchitecture")=>"String", Symbol("authenticationMethod")=>"AuthenticationMethod", Symbol("radiusServerAuthCertificate")=>"String", Symbol("clientRootCertificates")=>"Vector{String}")
Base.propertynames(::Type{ VpnClientParameters }) = collect(keys(_property_map_VpnClientParameters))
Swagger.property_type(::Type{ VpnClientParameters }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_VpnClientParameters[name]))}
Swagger.field_name(::Type{ VpnClientParameters }, property_name::Symbol) =  _property_map_VpnClientParameters[property_name]

const _allowed_VpnClientParameters_processorArchitecture = ["Amd64", "X86"]

function check_required(o::VpnClientParameters)
    true
end

function validate_property(::Type{ VpnClientParameters }, name::Symbol, val)
    if name === Symbol("processorArchitecture")
        Swagger.validate_param(name, "VpnClientParameters", :enum, val, _allowed_VpnClientParameters_processorArchitecture)
    end
end


text:
using Spring
using CUDA

Spring.do_synchronize(::Type{CuArray{Float64,2}}) = CUDA.synchronize()

function go()
    sp=SpringParam(ls=0.1,ms=1,ks=2.5,ns=1000)
    ip=InitParam(λ=2.0,shift=0.1,pos=sp.ls*sp.ns/2)
    ap=AnimParam(δt=0.1,nδt=5000,nδtperframe=100)

    CUDA.allowscalar(false)
    # V=Array{Float64,2}
    V=CuArray{Float64,2}

    @time animate_spring2D(sp,ip,ap,V)
end
go()


text:
import Images
using Color, Base.Test

# Comparison of each element in arrays with a scalar
approx_equal(ar, v) = all(abs(ar.-v) .< sqrt(eps(v)))
approx_equal(ar::Images.AbstractImage, v) = approx_equal(Images.data(ar), v)

# arithmetic
img = convert(Images.Image, zeros(3,3))
@assert Images.limits(img) == (0,1)
img2 = (img .+ 3)/2
@assert all(img2 .== 1.5)
@assert Images.limits(img2) == (1.5,2.0)
img3 = 2img2
@assert all(img3 .== 3)
img3 = copy(img2)
img3[img2 .< 4] = -1
@assert all(img3 .== -1)
img = convert(Images.Image, rand(3,4))
A = rand(3,4)
img2 = img .* A
@assert all(Images.data(img2) == Images.data(img).*A)
@assert Images.limits(img2) == (0,1)
img2 = convert(Images.Image, A)
img2 = img2 .- 0.5
img3 = 2img .* img2
@assert Images.limits(img3) == (-1, 1)
img2 = img ./ A
@assert Images.limits(img2) == (0, Inf)
img2 = (2img).^2
@assert Images.limits(img2) == (0, 4)
imgu = convert(Images.Image, Uint8[1 240; 10 128])  # from #101
@assert Images.limits(2imgu) == (0x00,0xff)

# scaling, ssd
img = convert(Images.Image, fill(typemax(Uint16), 3, 3))
scalei = Images.scaleinfo(Uint8, img)
img8 = scale(scalei, img)
@assert all(img8 .== typemax(Uint8))
mxA = -1.0
while mxA < 0
    A = randn(3,3)
    mxA = maximum(A)
end
offset = 30.0
img = convert(Images.Image, A .+ offset)
scalei = Images.ScaleMinMax{Uint8, Float64}(offset, offset+mxA, 100/mxA)
imgs = scale(scalei, img)
@assert minimum(imgs) == 0
@assert maximum(imgs) == 100
@assert eltype(imgs) == Uint8
imgs = Images.imadjustintensity(img, [])
mnA = minimum(A)
@assert Images.ssd(imgs, (A.-mnA)/(mxA-mnA)) < eps()
A = reshape(1:9, 3, 3)
B = scale(Images.ClipMin(Float32, 3), A)
@assert eltype(B) == Float32 && B == [3 4 7; 3 5 8; 3 6 9]
B = scale(Images.ClipMax(Uint8, 7), A)
@assert eltype(B) == Uint8 && B == [1 4 7; 2 5 7; 3 6 7]

# Reductions
let
    A = rand(5,5,3)
    img = Images.colorim(A, "RGB")
    img["limits"] = (0.0, 1.0)
    s12 = sum(img, (1,2))
    @test colorspace(s12) == "RGB"
    @test limits(s12) == (0.0,25.0)
    s3 = sum(img, (3,))
    @test colorspace(s3) == "Unknown"
    @test limits(s3) == (0.0,3.0)
end

# Array padding
let A = [1 2; 3 4]
    @test Images.padarray(A, (0,0), (0,0), "replicate") == A
    @test Images.padarray(A, (1,2), (2,0), "replicate") == [1 1 1 2; 1 1 1 2; 3 3 3 4; 3 3 3 4; 3 3 3 4]
    @test Images.padarray(A, [2,1], [0,2], "circular") == [2 1 2 1 2; 4 3 4 3 4; 2 1 2 1 2; 4 3 4 3 4]
    @test Images.padarray(A, (1,2), (2,0), "symmetric") == [2 1 1 2; 2 1 1 2; 4 3 3 4; 4 3 3 4; 2 1 1 2]
    @test Images.padarray(A, (1,2), (2,0), "value", -1) == [-1 -1 -1 -1; -1 -1 1 2; -1 -1 3 4; -1 -1 -1 -1; -1 -1 -1 -1]
    A = [1 2 3; 4 5 6]
    @test Images.padarray(A, (1,2), (2,0), "reflect") == [6 5 4 5 6; 3 2 1 2 3; 6 5 4 5 6; 3 2 1 2 3; 6 5 4 5 6]
    A = [1 2; 3 4]
    @test Images.padarray(A, (1,1)) == [1 1 2 2; 1 1 2 2; 3 3 4 4; 3 3 4 4]
    @test Images.padarray(A, (1,1), "replicate", "both") == [1 1 2 2; 1 1 2 2; 3 3 4 4; 3 3 4 4]
    @test Images.padarray(A, (1,1), "circular", "pre") == [4 3 4; 2 1 2; 4 3 4]
    @test Images.padarray(A, (1,1), "symmetric", "post") == [1 2 2; 3 4 4; 3 4 4]
    A = ["a" "b"; "c" "d"]
    @test Images.padarray(A, (1,1)) == ["a" "a" "b" "b"; "a" "a" "b" "b"; "c" "c" "d" "d"; "c" "c" "d" "d"]
end

# filtering
EPS = 1e-14
for T in (Float64, Int)
    A = zeros(T,3,3); A[2,2] = 1
    kern = rand(3,3)
    @test maximum(abs(Images.imfilter(A, kern) - rot180(kern))) < EPS
    kern = rand(2,3)
    @test maximum(abs(Images.imfilter(A, kern)[1:2,:] - rot180(kern))) < EPS
    kern = rand(3,2)
    @test maximum(abs(Images.imfilter(A, kern)[:,1:2] - rot180(kern))) < EPS
end
for T in (Float64, Int)
    # Separable kernels
    A = zeros(T,3,3); A[2,2] = 1
    kern = rand(3).*rand(3)'
    @test maximum(abs(Images.imfilter(A, kern) - rot180(kern))) < EPS
    kern = rand(2).*rand(3)'
    @test maximum(abs(Images.imfilter(A, kern)[1:2,:] - rot180(kern))) < EPS
    kern = rand(3).*rand(2)'
    @test maximum(abs(Images.imfilter(A, kern)[:,1:2] - rot180(kern))) < EPS
end
A = zeros(3,3); A[2,2] = 1
kern = rand(3,3)
@test maximum(abs(Images.imfilter_fft(A, kern) - rot180(kern))) < EPS
kern = rand(2,3)
@test maximum(abs(Images.imfilter_fft(A, kern)[1:2,:] - rot180(kern))) < EPS
kern = rand(3,2)
@test maximum(abs(Images.imfilter_fft(A, kern)[:,1:2] - rot180(kern))) < EPS

@assert approx_equal(Images.imfilter(ones(4,4), ones(3,3)), 9.0)
@assert approx_equal(Images.imfilter(ones(3,3), ones(3,3)), 9.0)
@assert approx_equal(Images.imfilter(ones(3,3), [1 1 1;1 0.0 1;1 1 1]), 8.0)
img = convert(Images.Image, ones(4,4))
@assert approx_equal(Images.imfilter(img, ones(3,3)), 9.0)
A = zeros(5,5,3); A[3,3,[1,3]] = 1
@assert Images.colordim(A) == 3
kern = rand(3,3)
kernpad = zeros(5,5); kernpad[2:4,2:4] = kern
Af = Images.imfilter(A, kern)

@test_approx_eq Af cat(3, rot180(kernpad), zeros(5,5), rot180(kernpad))
Aimg = permutedims(convert(Images.Image, A), [3,1,2])
@test_approx_eq Images.imfilter(Aimg, kern) permutedims(Af, [3,1,2])
@assert approx_equal(Images.imfilter(ones(4,4),ones(1,3),"replicate"), 3.0)

@assert approx_equal(Images.imfilter_gaussian(ones(4,4), [5,5]), 1.0)

A = zeros(Int, 9, 9); A[5, 5] = 1
@test maximum(abs(Images.imfilter_LoG(A, [1,1]) - imlog(1.0))) < EPS

# restriction
A = reshape(uint16(1:60), 4, 5, 3)
B = Images.restrict(A, (1,2))
@test_approx_eq B cat(3, [ 0.96875  4.625   5.96875;
                           2.875   10.5    12.875;
                           1.90625  5.875   6.90625],
                         [ 8.46875  14.625 13.46875;
                          17.875    30.5   27.875;
                           9.40625  15.875 14.40625],
                         [15.96875  24.625 20.96875;
                          32.875    50.5   42.875;
                          16.90625  25.875 21.90625])
A = reshape(1:60, 5, 4, 3)
B = Images.restrict(A, (1,2,3))
@test_approx_eq B cat(3, [ 2.6015625  8.71875 6.1171875;
                           4.09375   12.875   8.78125;
                           3.5390625 10.59375 7.0546875],
                         [10.1015625 23.71875 13.6171875;
                          14.09375   32.875   18.78125;
                          11.0390625 25.59375 14.5546875])

# color conversion
gray = linspace(0.0,1.0,5) # a 1-dimensional image
gray8 = iround(Uint8, 255*gray)
gray32 = [uint32(g)<<16 | uint32(g)<<8 | uint32(g) for g in gray8]
imgray = Images.Image(gray, ["colordim"=>0, "colorspace"=>"Gray"])
buf = Images.uint32color(imgray)
@assert buf == gray32
rgb = [RGB(g, g, g) for g in gray]
buf = Images.uint32color(rgb)
@assert buf == gray32
img = Images.Image(gray32, ["colordim"=>0, "colorspace"=>"RGB24"])
buf = Images.uint32color(img)
@assert buf == gray32
rgb = repeat(gray, outer=[1,3])
img = Images.Image(rgb, ["colordim"=>2, "colorspace"=>"RGB"])
buf = Images.uint32color(img)
@assert buf == gray32
rgb = repeat(gray', outer=[3,1])
img = Images.Image(rgb, ["colordim"=>1, "colorspace"=>"RGB"])
buf = Images.uint32color(img)
@assert buf == gray32
ovr = Images.Overlay((gray, 0*gray), (RGB(1,0,1), RGB(0,1,0)), ((0,1),(0,1)))
buf = Images.uint32color(ovr)
nogreen = [uint32(g)<<16 | uint32(g) for g in gray8]
@assert buf == nogreen
ovr = Images.Overlay((gray, gray), (RGB(1,0,1), RGB(0,1,0)), ((0,1),(0,1)))
ovr.visible[2] = false
buf = Images.uint32color(ovr)
@assert buf == nogreen

# erode/dilate
A = zeros(4,4,3)
A[2,2,1] = 0.8
A[4,4,2] = 0.6
Ae = Images.erode(A)
@assert Ae == zeros(size(A))
Ad = Images.dilate(A)
Ar = [0.8 0.8 0.8 0;
      0.8 0.8 0.8 0;
      0.8 0.8 0.8 0;
      0 0 0 0]
Ag = [0 0 0 0;
      0 0 0 0;
      0 0 0.6 0.6;
      0 0 0.6 0.6]
@assert Ad == cat(3, Ar, Ag, zeros(4,4))
Ae = Images.erode(Ad)
Ar = [0.8 0.8 0 0;
      0.8 0.8 0 0;
      0 0 0 0;
      0 0 0 0]
Ag = [0 0 0 0;
      0 0 0 0;
      0 0 0 0;
      0 0 0 0.6]
@assert Ae == cat(3, Ar, Ag, zeros(4,4))

# opening/closing
A = zeros(4,4,3)
A[2,2,1] = 0.8
A[4,4,2] = 0.6
Ao = Images.opening(A)
@assert Ao == zeros(size(A))
A = zeros(10,10)
A[4:7,4:7] = 1
B = copy(A)
A[5,5] = 0
Ac = Images.closing(A)
@assert Ac == B

# label_components
A = [true  true  false true;
     true  false true  true]
lbltarget = [1 1 0 2;
             1 0 2 2]
lbltarget1 = [1 2 0 4;
              1 0 3 4]
@assert Images.label_components(A) == lbltarget
@assert Images.label_components(A, [1]) == lbltarget1
connectivity = [false true  false;
                true  false true;
                false true  false]
@assert Images.label_components(A, connectivity) == lbltarget
connectivity = trues(3,3)
lbltarget2 = [1 1 0 1;
              1 0 1 1]
@assert Images.label_components(A, connectivity) == lbltarget2

# phantoms

P = [ 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0;
      0.0  0.0  1.0  0.2  0.2  1.0  0.0  0.0;
      0.0  0.0  0.2  0.3  0.3  0.2  0.0  0.0;
      0.0  0.0  0.2  0.0  0.2  0.2  0.0  0.0;
      0.0  0.0  0.2  0.0  0.0  0.2  0.0  0.0;
      0.0  0.0  0.2  0.2  0.2  0.2  0.0  0.0;
      0.0  0.0  1.0  0.2  0.2  1.0  0.0  0.0;
      0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0 ]

Q = Images.shepp_logan(8)
@assert norm((P-Q)[:]) < 1e-10

P = [ 0.0  0.0  0.0   0.0   0.0   0.0   0.0  0.0;
      0.0  0.0  2.0   1.02  1.02  2.0   0.0  0.0;
      0.0  0.0  1.02  1.03  1.03  1.02  0.0  0.0;
      0.0  0.0  1.02  1.0   1.02  1.02  0.0  0.0;
      0.0  0.0  1.02  1.0   1.0   1.02  0.0  0.0;
      0.0  0.0  1.02  1.02  1.02  1.02  0.0  0.0;
      0.0  0.0  2.0   1.02  1.02  2.0   0.0  0.0;
      0.0  0.0  0.0   0.0   0.0   0.0   0.0  0.0 ]

Q = Images.shepp_logan(8,highContrast=false)
@assert norm((P-Q)[:]) < 1e-10