text:
using PyCall

println("Running build.jl for the PicoQuant package.")

# Change that to whatever packages you need.
const PACKAGES = ["qiskit==0.19.2", "pylatexenc==2.1"]

# Use eventual proxy info
proxy_arg=String[]
if haskey(ENV, "http_proxy")
    push!(proxy_arg, "--proxy")
    push!(proxy_arg, ENV["http_proxy"])
end

# Import pip
try
    pyimport("pip")
catch
    # If it is not found, install it
    println("Pip not found on your sytstem. Downloading it.")
    get_pip = joinpath(dirname(@__FILE__), "get-pip.py")
    download("https://bootstrap.pypa.io/get-pip.py", get_pip)
    run(`$(PyCall.python) $(proxy_arg) $get_pip --user`)
end

println("Installing required python packages using pip")
run(`$(PyCall.python) $(proxy_arg) -m pip install --upgrade pip setuptools`)
run(`$(PyCall.python) $(proxy_arg) -m pip install $(PACKAGES)`)


text:
export create_paramops, create_curls, create_πcmps
export create_iF′ₗgen, create_iHₗgen, create_iEₗgen
export create_βF′ₜgen, create_βHₜgen, create_βEₜgen
export create_A
export calc_mode

function create_stretched_∆ls(mdl::Model)
    ωpml = mdl.ωpml
    grid = mdl.grid
    Npml = mdl.Npml
    s∆l = create_stretched_∆l(ωpml, grid, Npml)
    s∆l⁻¹ = invert_∆l(s∆l)

    boundft = mdl.boundft
    gtₑ = ft2gt.(EE, boundft)
    gtₘ = ft2gt.(HH, boundft)

    s∆lₑ = t_ind(s∆l, gtₑ)  # stretched ∆l's centered at E-field plane locations
    s∆lₘ = t_ind(s∆l, gtₘ)  # stretched ∆l's centered at H-field plane locations
    s∆lₑ⁻¹ = t_ind(s∆l⁻¹, gtₑ)
    s∆lₘ⁻¹ = t_ind(s∆l⁻¹, gtₘ)

    return (s∆lₑ, s∆lₘ), (s∆lₑ⁻¹, s∆lₘ⁻¹)
end

function create_paramops(mdl::Model{K,Kₑ,Kₘ,Kₑₜ,Kₘₜ}) where {K,Kₑ,Kₘ,Kₑₜ,Kₘₜ}
    (s∆lₑ, s∆lₘ), (s∆lₑ⁻¹, s∆lₘ⁻¹) = create_stretched_∆ls(mdl)
    calc_matparams!(mdl)  # assignment and smoothing; implemented for each specialized alias of Model

    boundft = mdl.boundft
    isbloch = mdl.grid.isbloch
    e⁻ⁱᵏᴸ = create_e⁻ⁱᵏᴸ(mdl)

    isfwd_inₑ = boundft.!=EE
    isfwd_inₘ = boundft.!=HH

    Pεₜ = Kₑₜ==1 ? create_paramop(mdl.εₜarr; mdl.order_cmpfirst) :
                   create_paramop(mdl.εₜarr, isfwd_inₑ, s∆lₘ, s∆lₑ⁻¹, isbloch, e⁻ⁱᵏᴸ; mdl.order_cmpfirst)
    Pμₜ = Kₘₜ==1 ? create_paramop(mdl.μₜarr; mdl.order_cmpfirst) :
                   create_paramop(mdl.μₜarr, isfwd_inₘ, s∆lₑ, s∆lₘ⁻¹, isbloch, e⁻ⁱᵏᴸ; mdl.order_cmpfirst)

    # The following lines work even if the model does not deal with the corresponding
    # longitudinal components.
    Pεₗ = create_paramop(mdl.εₗarr; mdl.order_cmpfirst)  # 0×0 matrix for Kₑₗ = 0 (no longitudinal E-field)
    Pμₗ = create_paramop(mdl.μₗarr; mdl.order_cmpfirst)  # 0×0 matrix for Kₘₗ = 0 (no longitudinal H-field)

    return (Pεₜ,Pμₜ), (Pεₗ,Pμₗ)
end

function create_curls(mdl::Model)
    _, (s∆lₑ⁻¹, s∆lₘ⁻¹) = create_stretched_∆ls(mdl)

    boundft = mdl.boundft

    cmpₛ = mdl.cmpₛ
    cmpₑₜ, cmpₘₜ = mdl.cmpₑₜ, mdl.cmpₘₜ
    cmpₑₗ, cmpₘₗ = mdl.cmpₑₗ, mdl.cmpₘₗ

    isbloch = mdl.grid.isbloch
    e⁻ⁱᵏᴸ = create_e⁻ⁱᵏᴸ(mdl)

    isfwdₑ = boundft.==EE
    isfwdₘ = boundft.==HH

    # The following lines work even if the model does not deal with the corresponding
    # longitudinal components.
    ∇̽ₑₜ = create_curl(isfwdₑ, s∆lₘ⁻¹, isbloch, e⁻ⁱᵏᴸ; cmp_shp=cmpₛ, cmp_out=cmpₘₗ, cmp_in=cmpₑₜ, mdl.order_cmpfirst)
    ∇̽ₘₜ = create_curl(isfwdₘ, s∆lₑ⁻¹, isbloch, e⁻ⁱᵏᴸ; cmp_shp=cmpₛ, cmp_out=cmpₑₗ, cmp_in=cmpₘₜ, mdl.order_cmpfirst)

    ∇̽ₑₗ = create_curl(isfwdₑ, s∆lₘ⁻¹, isbloch, e⁻ⁱᵏᴸ; cmp_shp=cmpₛ, cmp_out=cmpₘₜ, cmp_in=cmpₑₗ, mdl.order_cmpfirst)
    ∇̽ₘₗ = create_curl(isfwdₘ, s∆lₑ⁻¹, isbloch, e⁻ⁱᵏᴸ; cmp_shp=cmpₛ, cmp_out=cmpₑₜ, cmp_in=cmpₘₗ, mdl.order_cmpfirst)

    return (∇̽ₑₜ,∇̽ₘₜ), (∇̽ₑₗ,∇̽ₘₗ)
end

function create_πcmps(mdl::Model{K}) where {K}
    grid = mdl.grid
    cmpₑₜ, cmpₘₜ = mdl.cmpₑₜ, mdl.cmpₘₜ

    permute = reverse(SVec(ntuple(identity, Val(K))))
    scale₀ = SFloat(1,-1)

    πcmpₑ = create_πcmp(grid.N, permute, scale=scale₀[cmpₑₜ]; mdl.order_cmpfirst)
    πcmpₘ = create_πcmp(grid.N, permute, scale=scale₀[cmpₘₜ]; mdl.order_cmpfirst)

    return πcmpₑ, πcmpₘ
end

create_iHₗgen(ω, Ps::Tuple22{AbsMatNumber}, ∇̽s::Tuple22{AbsMatNumber}) = create_iF′ₗgen(EE, ω, Ps, ∇̽s)
create_iEₗgen(ω, Ps::Tuple22{AbsMatNumber}, ∇̽s::Tuple22{AbsMatNumber}) = create_iF′ₗgen(HH, ω, Ps, ∇̽s)

create_iF′ₗgen(ft, ω, Ps::Tuple22{AbsMatNumber}, ∇̽s::Tuple22{AbsMatNumber}) =
    create_iF′ₗgen(ft, ω, Ps[2], ∇̽s[1])

# Create the operator that generates im * F′ₗ from Fₜ, where F′ is the field complementary
# to F (e.g., F′ = H-field for F = E-field).
function create_iF′ₗgen(ft::FieldType,  # type of input field Fₜ, not output field F′ₗ
                        ω::Number,
                        Pₗs::Tuple2{AbsMatNumber},
                        ∇̽ₜs::Tuple2{AbsMatNumber})
    nft = Int(ft)
    nft′ = alter(nft)

    α = ft==EE ? -1/ω : 1/ω

    P′ₗ = Pₗs[nft′]
    ∇̽ₜ = ∇̽ₜs[nft]

    L = α .* (P′ₗ \ ∇̽ₜ)

    return L
end

create_βHₜgen(ω, Ps, ∇̽s, πcmps) = create_βF′ₜgen(EE, ω, Ps, ∇̽s, πcmps)
create_βEₜgen(ω, Ps, ∇̽s, πcmps) = create_βF′ₜgen(HH, ω, Ps, ∇̽s, πcmps)

# Create the operator that generates β * F′ₜ, where F′ is the field complementary to F (e.g.,
# F′ = H-field for F = E-field).
create_βF′ₜgen(ft, ω, Ps, ∇̽s, πcmps) = sum(create_βF′ₜgen_comps(ft, ω, Ps, ∇̽s, πcmps))

function create_βF′ₜgen_comps(ft::FieldType,  # type of input field Fₜ, not output field F′ₗ
                              ω::Number,
                              Ps::Tuple22{AbsMatNumber},
                              ∇̽s::Tuple22{AbsMatNumber},
                              πcmps::Tuple2{AbsMatNumber})
    nft = Int(ft)
    nft′ = alter(nft)

    α = ft==EE ? ω : -ω

    Pₜs, Pₗs = Ps
    ∇̽ₜs, ∇̽ₗs = ∇̽s

    Pₜ = Pₜs[nft]
    ∇̽′ₗ = ∇̽ₗs[nft′]
    πcmp = πcmps[nft]

    βF′ₜgen₁ = (α .* πcmp) * Pₜ

    if length(Pₗs[nft′]) > 0
        iF′ₗgen = create_iF′ₗgen(ft, ω, Pₗs, ∇̽ₜs)
        βF′ₜgen₂ = (πcmp * ∇̽′ₗ) * iF′ₗgen
    else
        βF′ₜgen₂ = spzeros(eltype(βF′ₜgen₁), size(βF′ₜgen₁)...)
    end

    return βF′ₜgen₁, βF′ₜgen₂
end

function create_A(ft::FieldType,  # type of input field Fₜ
                  ω::Number,
                  Ps::Tuple22{AbsMatNumber},
                  ∇̽s::Tuple22{AbsMatNumber},
                  πcmps::Tuple2{AbsMatNumber})
    ft′ = alter(ft)

    # Old implementation.  Though more straightforward, this implementation does not
    # gurantee the vector calculus identity βFₜgen₂ * βF′ₜgen₂ = 0 due to rounding errors.
    #
    # βF′ₜgen = create_βF′ₜgen(ft, ω, Ps, ∇̽s, πcmps)
    # βFₜgen = create_βF′ₜgen(ft′, ω, Ps, ∇̽s, πcmps)
    # A = βFₜgen * βF′ₜgen

    βF′ₜgen₁, βF′ₜgen₂ = create_βF′ₜgen_comps(ft, ω, Ps, ∇̽s, πcmps)
    βFₜgen₁, βFₜgen₂ = create_βF′ₜgen_comps(ft′, ω, Ps, ∇̽s, πcmps)

    A = βFₜgen₁ * βF′ₜgen₁ + βFₜgen₁ * βF′ₜgen₂ + βFₜgen₂ * βF′ₜgen₁  # βFₜgen₂ * βF′ₜgen₂ = 0 always

    return A
end

function calc_mode(mdl::Model, ω::Real, βguess::Number;
                   nmode::Int=1,  # mode number
                   fₜref::AbsVecNumber=ComplexF[])  # reference field w.r.t. which optimal phase angle is calculated
    ## Create the eigenequation and solve it.
    Ps = create_paramops(mdl)
    ∇̽s = create_curls(mdl)
    πcmps = create_πcmps(mdl)

    ft_eq = mdl.ft_eq
    A = create_A(ft_eq, ω, Ps, ∇̽s, πcmps)

    β², f = eigs(A, nev=nmode, sigma=βguess^2)

    fₜ = f[:,nmode]
    β = .√β²[nmode]

    # See RN - NLSE > December 12, 2021
    if length(fₜref)==length(fₜ)
        α = fₜref⋅fₜ
        θ = -angle(α)
        fₜ .*= exp(im*θ)
    end

    E, H = complete_fields(β, fₜ, ft_eq, ω, Ps, ∇̽s, πcmps, mdl)
    normalize!(E, H, mdl)

    return (β=β, E=E, H=H, f=fₜ)
end

ext: jl


text:
@memoize function version()
    buf = Vector{Cchar}(undef, NVML_SYSTEM_NVML_VERSION_BUFFER_SIZE)
    nvmlSystemGetNVMLVersion(pointer(buf), length(buf))

    # the version string is too long for Julia to handle, e.g. 11.450.36.06,
    # so split off the driver part into the build suffix
    ver = unsafe_string(pointer(buf))
    parts = parse.(Int, split(ver, '.'))
    return VersionNumber(parts[1], 0, 0, (), Tuple(parts[2:end]))
end

@memoize function driver_version()
    buf = Vector{Cchar}(undef, NVML_SYSTEM_DRIVER_VERSION_BUFFER_SIZE)
    nvmlSystemGetDriverVersion(pointer(buf), length(buf))
    return VersionNumber(unsafe_string(pointer(buf)))
end

@memoize function cuda_driver_version()
    ref = Ref{Cint}()
    nvmlSystemGetCudaDriverVersion_v2(ref)
    major, ver = divrem(ref[], 1000)
    minor, patch = divrem(ver, 10)
    return VersionNumber(major, minor, patch)
end


text:
function equals(obj::Object, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function get_absolute_accuracy(obj::BaseAbstractUnivariateIntegrator)
    return jcall(obj, "getAbsoluteAccuracy", jdouble, ())
end

function get_class(obj::Object)
    return jcall(obj, "getClass", Class, ())
end

function get_evaluations(obj::BaseAbstractUnivariateIntegrator)
    return jcall(obj, "getEvaluations", jint, ())
end

function get_iterations(obj::BaseAbstractUnivariateIntegrator)
    return jcall(obj, "getIterations", jint, ())
end

function get_maximal_iteration_count(obj::BaseAbstractUnivariateIntegrator)
    return jcall(obj, "getMaximalIterationCount", jint, ())
end

function get_minimal_iteration_count(obj::BaseAbstractUnivariateIntegrator)
    return jcall(obj, "getMinimalIterationCount", jint, ())
end

function get_relative_accuracy(obj::BaseAbstractUnivariateIntegrator)
    return jcall(obj, "getRelativeAccuracy", jdouble, ())
end

function hash_code(obj::Object)
    return jcall(obj, "hashCode", jint, ())
end

function integrate(obj::BaseAbstractUnivariateIntegrator, arg0::jint, arg1::UnivariateFunction, arg2::jdouble, arg3::jdouble)
    return jcall(obj, "integrate", jdouble, (jint, UnivariateFunction, jdouble, jdouble), arg0, arg1, arg2, arg3)
end

function notify(obj::Object)
    return jcall(obj, "notify", void, ())
end

function notify_all(obj::Object)
    return jcall(obj, "notifyAll", void, ())
end

function to_string(obj::Object)
    return jcall(obj, "toString", JString, ())
end

function wait(obj::Object)
    return jcall(obj, "wait", void, ())
end

function wait(obj::Object, arg0::jlong)
    return jcall(obj, "wait", void, (jlong,), arg0)
end

function wait(obj::Object, arg0::jlong, arg1::jint)
    return jcall(obj, "wait", void, (jlong, jint), arg0, arg1)
end


text:
corpus_dir = joinpath(@__DIR__, "corpus", "tests")
corpus_files = filter(x -> splitext(x)[2] == ".json", readdir(corpus_dir))

@testset "corpus" begin

@testset "roundtrip $(splitext(file)[1])" for file in corpus_files
    doc = JSON3.read(read(joinpath(corpus_dir, file)))
    t = parse(UInt8, doc["bson_type"])
    if haskey(doc, "valid")
        @testset """roundtrip $file $(test_case["description"]) $validator""" for test_case in doc["valid"], validator in [UncheckedBSONValidator(), LightBSONValidator(), StrictBSONValidator()]
            src = hex2bytes(test_case["canonical_bson"])
            x = BSONReader(src, StrictBSONValidator())[Any]
            dst = UInt8[]
            writer = BSONWriter(dst)
            writer[] = x
            close(writer)
            @test dst == src
        end
    end
    if haskey(doc, "decodeErrors")
        @testset """decode error $file $(test_case["description"])""" for test_case in doc["decodeErrors"]
            src = hex2bytes(test_case["bson"])
            @test_throws Exception BSONReader(src, StrictBSONValidator())[Any]
        end
    end
end

end


text:
"""
    AbstractSimplex{D, T, I<:Integer} <: AbstractCell{I}

An abstract type for representing simplices. A simplex is an `AbstractCell` with an integer
index. A `D`-simplex has `D + 1` vertices.

If a type implements the interface below, default implementations of [`boundary`](@ref),
[`coboundary`](@ref), and [`vertices`](@ref) are provided.

# Interface

* `AbstractSimplex{0}(::I, ::T)`
* `AbstractSimplex{1}(::I, ::T)`
* [`birth(::AbstractSimplex)`](@ref)
* [`index(::AbstractSimplex)`](@ref)
* [`sign(::AbstractSimplex)`](@ref)
* [`Base.:-(::AbstractSimplex)`](@ref)
* [`unsafe_simplex`](@ref)
* [`unsafe_cofacet`](@ref)

"""
abstract type AbstractSimplex{D,T,I} <: AbstractCell{D,T,I} end

function Base.show(io::IO, sx::AbstractSimplex{D}) where {D}
    sgn = sign(sx) == 1 ? :+ : :-
    name = nameof(typeof(sx))
    return print(io, "$sgn$name{$D}($(vertices(sx)), $(birth(sx)))")
end
function Base.show(io::IO, ::MIME"text/plain", sx::AbstractSimplex{D}) where {D}
    sgn = sign(sx) == 1 ? :+ : :-
    name = nameof(typeof(sx))
    return print(
        io,
        "$D-dimensional $name(index=$(index(sx)), birth=$(birth(sx))):\n  $sgn",
        vertices(sx),
    )
end

###
### Vertices and indices
###
"""
    _binomial(n, ::Val{K})

Binomial coefficients for small, statically known values of `K`, where `n` and `K` are
always positive. Using `Val(K)` allows the compiler to optimize away the loop.
"""
_binomial(::I, ::Val{0}) where {I} = one(I)
_binomial(n, ::Val{1}) = n
function _binomial(n::I, ::Val{K}) where {K,I}
    x = nn = I(n - K + 1)
    nn += one(I)
    for rr in I(2):I(K)
        x = div(x * nn, rr)
        nn += one(I)
    end
    return I(x)
end

"""
    _first_vertex(index, ::Val{K})

Use binary search to find index of first vertex in `(K - 1)`-dimensional simplex with index
`index`.
"""
function _first_vertex(index::I, ::Val{K}) where {I,K}
    lo = I(K - 1)
    hi = I(K + 100)
    while _binomial(hi, Val(K)) ≤ index
        lo = hi
        hi <<= 0x01
        hi + one(I) < lo && throw(OverflowError("simplex overflowed! This is a bug"))
    end
    return _first_vertex(index, Val(K), hi + one(I), lo)
end
function _first_vertex(index::I, ::Val{K}, hi::I, lo::I=I(K - 1)) where {I,K}
    while lo < hi - one(I)
        m = lo + ((hi - lo) >>> 0x01)
        if _binomial(m, Val(K)) ≤ index
            lo = m
        else
            hi = m
        end
    end
    return lo
end
function _first_vertex(index::I, ::Val{1}, ::I=I(0), ::I=I(0)) where {I}
    return index
end
function _first_vertex(index::I, ::Val{2}, ::I=I(0), ::I=I(0)) where {I}
    # This is https://oeis.org/A002024
    return floor(I, (√(8 * index + 1) + 1) / 2)
end

"""
    _vertices(index::I, ::Val{N})::NTuple{N,I}

Get the vertices of simplex represented by index.
"""
@inline function _vertices(index::I, ::Val{K}) where {I,K}
    index -= one(I)
    vk = _first_vertex(index, Val(K))
    index -= _binomial(vk, Val(K))
    return tuple(vk + one(I), _vertices(index, Val(K - 1), vk)...)
end
@inline function _vertices(index::I, ::Val{K}, prev) where {I,K}
    vk = _first_vertex(index, Val(K), prev)
    index -= _binomial(vk, Val(K))
    return tuple(vk + one(I), _vertices(index, Val(K - 1), vk)...)
end
@inline _vertices(index::I, ::Val{1}, _) where {I} = (index + one(I),)
@inline _vertices(index, ::Val{1}) = (index,)

"""
    index(vertices)

Calculate the index from tuple or static vector of vertices. The index is equal to

```math
(i_d, i_{d-1}, ..., i_1) \\mapsto \\sum_{k=1}^{d+1} \\binom{i_k - 1}{k},
```

where ``i_k`` are the simplex vertex indices.

```jldoctest
julia> index((6,2,1))
11
```
"""
index(vertices::SVector) = index(Tuple(vertices))
index(vertices::NTuple) = _index(vertices, one(eltype(vertices)))

@inline _index(::NTuple{0}, acc) = acc
@inline function _index(vertices::NTuple{N,I}, acc::I) where {N,I}
    acc += _binomial(first(vertices) - one(I), Val(N))
    return _index(TupleTools.tail(vertices), acc)
end

"""
    index_overflow_check(vertices[, message])

Check if index overflows for a particular collection of vertices. Throw an error if it does.
"""
function index_overflow_check(
    vertices,
    message="simplex $vertices overflows in $(I)! " *
            "Try using a bigger index type in your filtration.",
)
    idx = index(vertices)
    idx_big = index(BigInt.(vertices))
    if idx ≠ idx_big
        throw(OverflowError(message))
    end
end

vertices(sx::AbstractSimplex) = _vertices(index(sx), Val(length(sx)))
Base.length(::Type{<:AbstractSimplex{D}}) where {D} = D + 1

###
### Coboundaries
###
function coboundary(
    filtration::AbstractFiltration, simplex::AbstractSimplex, ::Val{A}=Val(true)
) where {A}
    return Coboundary{A}(filtration, simplex)
end

struct Coboundary{A,D,I,F,S}
    filtration::F
    simplex::S
    vertices::NTuple{D,I}

    function Coboundary{A}(
        filtration::F, simplex::AbstractSimplex{D,<:Any,I}
    ) where {A,D,I,F}
        S = typeof(simplex)
        return new{A,D + 1,I,F,S}(filtration, simplex, Tuple(vertices(simplex)))
    end
end

function Base.iterate(
    ci::Coboundary{A,D,I}, (v, k)=(I(nv(ci.filtration) + 1), D)
) where {A,D,I}
    @inbounds while true
        v -= one(I)
        while k > 0 && v == ci.vertices[end + 1 - k]
            A || return nothing
            v -= one(I)
            k -= 1
        end
        v > 0 || return nothing
        new_vertices = TupleTools.insertafter(ci.vertices, D - k, (v,))
        sx = unsafe_cofacet(ci.filtration, ci.simplex, new_vertices, v)
        if !isnothing(sx)
            _sx::simplex_type(ci.filtration, D) = sx
            if isodd(k)
                _sx = -_sx
            end
            return _sx, (v, k)
        end
    end
end

struct SparseCoboundary{A,D,I,F,S}
    filtration::F
    simplex::S
    vertices::NTuple{D,I}
    ptrs_begin::SVector{D,I}
    ptrs_end::SVector{D,I}

    function SparseCoboundary{A}(
        filtration::F, simplex::AbstractSimplex{D,<:Any,I}
    ) where {A,D,I,F}
        verts = vertices(simplex)
        adj = adjacency_matrix(filtration)
        colptr = adj.colptr
        ptrs_begin = colptr[SVector(verts) .+ 1]
        ptrs_end = colptr[SVector(verts)]
        return new{A,D + 1,I,F,typeof(simplex)}(
            filtration, simplex, verts, ptrs_begin, ptrs_end
        )
    end
end

@propagate_inbounds @inline function _next_common(
    ptrs::SVector{D,I}, ptrs_end::SVector{D,I}, rowval
) where {D,I}
    ptrs = ptrs .- one(I)
    for i in 1:D
        ptrs[i] < ptrs_end[i] && return zero(ptrs), 0
    end
    m = I(rowval[ptrs[2]])
    i = 1
    while true
        ptrs_i = ptrs[i]
        row = I(rowval[ptrs_i])
        while row > m
            ptrs -= SVector(ntuple(isequal(i), Val(D)))
            ptrs_i -= one(I)
            ptrs_i < ptrs_end[i] && return zero(ptrs), zero(eltype(rowval))
            row = I(rowval[ptrs_i])
        end
        i = ifelse(row == m, i + 1, 1)
        i > D && return ptrs, m
        m = row
    end
end

function Base.iterate(
    it::SparseCoboundary{A,D,I}, (ptrs, k)=(it.ptrs_begin, D)
) where {A,D,I}
    adj = adjacency_matrix(it.filtration)
    rowval = adj.rowval
    nzval = adj.nzval
    @inbounds while true
        ptrs, v = _next_common(ptrs, it.ptrs_end, rowval)
        if iszero(first(ptrs))
            return nothing
        elseif k > 0 && v == it.vertices[end + 1 - k]
            k -= 1
        else
            while k > 0 && v < it.vertices[end + 1 - k]
                k -= 1
            end
            !A && k ≠ D && return nothing

            new_vertices = TupleTools.insertafter(it.vertices, D - k, (v,))
            new_edges = nzval[ptrs]
            sx = unsafe_cofacet(it.filtration, it.simplex, new_vertices, v, new_edges)
            if !isnothing(sx)
                _sx::simplex_type(it.filtration, D) = sx
                if isodd(k)
                    _sx = -sx
                end
                return _sx, (ptrs, k)
            end
        end
    end
end

###
### Boundary
###
function boundary(filtration::AbstractFiltration, simplex::AbstractSimplex)
    return Boundary(filtration, simplex)
end

struct Boundary{D,I,F,S}
    filtration::F
    simplex::S
    vertices::NTuple{D,I}

    function Boundary(filtration::F, simplex::AbstractSimplex{D,<:Any,I}) where {D,I,F}
        S = typeof(simplex)
        return new{D + 1,I,F,S}(filtration, simplex, Tuple(vertices(simplex)))
    end
end

function Base.iterate(bi::Boundary{D,I}, k=1) where {D,I}
    while k ≤ D
        facet_vertices = TupleTools.deleteat(bi.vertices, k)
        k += 1
        sx = unsafe_simplex(bi.filtration, Val(D - 2), facet_vertices)
        if !isnothing(sx)
            _sx::simplex_type(bi.filtration, D - 2) = sx
            if isodd(k)
                _sx = -_sx
            end
            return _sx, k
        end
    end
    return nothing
end


text:
using MuZero
using Base.Threads: @spawn

include("game.jl")
include("inference_tests.jl")
cd(@__DIR__)
config = getconfig(Tac)
config.batch_size = 32
config.root_dirichlet_alpha=0.3
config.num_actors = 3
config.weight_decay = 1e-4
config.lr_init=0.002f0
config.save_interval = 10000
config.num_simulations = 100
config.num_unroll_steps = 1
config.training_steps = 200001

p_tests = getpolicytests(Tac)
p_tests.period = 250

learner = Learner(
    Tac;
    dir="output",
    config=config,
    tests=[p_tests]
)

@spawn selfplay!(learner)
@spawn selfplay!(learner)
@spawn selfplay!(learner)
length(learner.replay_buffer)

@show @time fit!(learner; steps=3)
@spawn fit!(learner)

learner.stats.losses |> length
learner |> gettraindf |> display
learner |> gettraindf |> mean1000 |> display
learner |> gettraindfall |> mean1000 |> display

learner |> show_pv_tests
learner |> show_reward_tests

#@learnertest 4854 |> svg
#@learnersample 33

selfplay!(learner; limit=10)

learner |> show_game

# still a fence post problem here
s = learner |> randtrainingsample
learner |> randtrainingsample |> svg

learner |> randtrainingbatch

using Base.Threads: @spawn
for i = 1:num_actors(learner)
    @spawn selfplay!(learner)
end

length(learner.replay_buffer)

@show get_steps(learner), learner.stats.losses, learner.stats.accuracies, learner.times

@time fit!(learner, steps=1)

learner |> gettraindfall |> display

@show learner.stats.losses, learner.stats.accuracies, learner.times

get_steps(learner), get_accuracy(learner)


text:
using ParserCombinator
import Base.==

reader_table = Dict{Symbol, Function}()

expr         = Delayed()
floaty_dot   = p"[-+]?[0-9]*\.[0-9]+([eE][-+]?[0-9]+)?[Ff]" > (x -> parse(Float32, x[1:end-1]))
floaty_nodot = p"[-+]?[0-9]*[0-9]+([eE][-+]?[0-9]+)?[Ff]" > (x -> parse(Float32, x[1:end-1]))
floaty       = floaty_dot | floaty_nodot
white_space  = p"([\s\n\r]*(?<!\\);[^\n\r$]+[\n\r\s$]*|[\s\n\r]+)"
opt_ws       = white_space | e""

doubley      = p"[-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)?[dD]" > (x -> parse(Float64, x[1:end-1]))

inty         = p"[-+]?\d+" > (x -> parse(Int, x))

uchary       = p"\\(u[\da-fA-F]{4})" > (x -> first(unescape_string(x)))
achary       = p"\\[0-7]{3}" > (x -> unescape_string(x)[1])
chary        = p"\\." > (x -> x[2])

stringy      = p"(?<!\\)\".*?(?<!\\)\"" > (x -> x[2:end-1]) #_0[2:end-1] } #r"(?<!\\)\".*?(?<!\\)"
booly        = p"(true|false)" > (x -> x == "true" ? true : false)
symboly      = p"[^\d(){}#'`,@~;~\[\]^\s][^\s()#'`,@~;^{}~\[\]]*" > Symbol
macrosymy    = p"@[^\d(){}#'`,@~;~\[\]^\s][^\s()#'`,@~;^{}~\[\]]*" > Symbol

sexpr        = E"(" + ~opt_ws + Repeat(expr + ~opt_ws) + E")" |> (x -> s_expr(x))
hashy        = E"#{" + ~opt_ws + Repeat(expr + ~opt_ws) + E"}" |> (x -> Set(x))
curly        = E"{" + ~opt_ws + Repeat(expr + ~opt_ws) + E"}" |> (x -> Dict(x[i] => x[i+1] for i = 1:2:length(x)))
dispatchy    = E"#" + symboly + ~opt_ws + expr |> (x -> reader_table[x[1]](x[2]))
bracket      = E"[" + ~opt_ws + Repeat(expr + ~opt_ws) + E"]" |> (x -> s_expr(x)) # TODO: not quite right
quot         = E"'" + expr > (x -> sx(:quote, x))
quasi        = E"`" + expr > (x -> sx(:quasi, x))
tildeseq     = E"~@" + expr > (x -> sx(:splice_seq, x))
tilde        = E"~" + expr > (x -> sx(:splice, x))

expr.matcher = doubley | floaty | inty | uchary | achary | chary | stringy | booly | symboly |
               macrosymy | dispatchy | sexpr | hashy | curly | bracket |
               quot | quasi | tildeseq | tilde

top_level    = Repeat(~opt_ws + expr) + ~opt_ws + Eos()

function read(str)
  x = parse_one(str, expr)
  x[1]
end


text:
module CompositeGrids
using StaticArrays

include("old/grid.jl")
export Grid

include("grid/simple.jl")
const SimpleGrid = SimpleG # alias for older convention
const AbstractGrid, OpenGrid, ClosedGrid = SimpleGrid.AbstractGrid, SimpleGrid.OpenGrid, SimpleGrid.ClosedGrid
export SimpleG, SimpleGrid, AbstractGrid, OpenGrid, ClosedGrid, denseindex

include("grid/composite.jl")
const CompositeGrid = CompositeG
export CompositeG, CompositeGrid

include("grid/interpolate.jl")
export Interp

end # module


text:
@testset "Moving max-min functions" begin

    len = 1000
    Random.seed!(0)
    T = Float64
    vec = rand(T, len)
    window = 15

    mx_ref = similar(vec)
    mn_ref = similar(vec)

    for i = 1:len
        mx_ref[i] = maximum(vec[max(1, i-window+1) : i])
        mn_ref[i] = minimum(vec[max(1, i-window+1) : i])
    end
    range_ref = mx_ref .- mn_ref

    @testset "movmaxmin" begin
        mx, mn = movmaxmin(vec, window)
        @test mx == mx_ref
        @test mn == mn_ref
        fill!(mx, 0)
        fill!(mn, 0)
        movmaxmin!(mx, mn, vec, window)
        @test mx == mx_ref
        @test mn == mn_ref
    end

    @testset "movrange" begin
        range = movrange(vec, window)
        @test range == range_ref
        fill!(range, 0)
        movrange!(range, vec, window)
        @test range == range_ref
        fill!(range, 0)
    end

    @testset "movmax" begin
        mx = movmax(vec, window)
        @test mx == mx_ref
        fill!(mx, 0)
        movmax!(mx, vec, window)
        @test mx == mx_ref
    end

    @testset "movmin" begin
        mn = movmin(vec, window)
        @test mn == mn_ref
        fill!(mn, 0)
        movmin!(mn, vec, window)
        @test mn == mn_ref
    end

end
# @test_throws BoundsError first(buf)


text:
using Colors, FixedPointNumbers
using Test
using ColorTypes: eltype_default, parametric3

struct C3{T} <: Color{T,3}
    c1::T; c2::T; c3::T;
end
struct C4{T} <: Color{T,4}
    c1::T; c2::T; c3::T; c4::T;
end
function ColorTypes._convert(::Type{Cdest}, ::Type{Odest}, ::Type{Osrc}, c) where {Cdest<:AbstractRGB,Odest,Osrc<:C4}
    Cdest(c.c1 + c.c4, c.c2 + c.c4, c.c3 + c.c4)
end
function ColorTypes._convert(::Type{Cdest}, ::Type{Odest}, ::Type{Osrc}, c) where {Cdest<:C4,Odest,Osrc<:AbstractRGB}
    r, g, b = red(c), green(c), blue(c)
    c4 = min(r, g, b)
    Cdest(r - c4, g - c4, b - c4, c4)
end

@testset "Conversion" begin
    r8(x) = reinterpret(N0f8, x)

    # srgb_compand / invert_srgb_compand
    @test Colors.srgb_compand(0.5) ≈ 0.7353569830524494 atol=eps()
    @test Colors.invert_srgb_compand(0.7353569830524494) ≈ 0.5 atol=eps()
    @test Colors.invert_srgb_compand(0.735357f0) ≈ 0.5f0 atol=eps(Float32)

    fractional_types = (RGB, BGR, XRGB, RGBX)  # types that support Fractional

    redF64 = RGB{Float64}(1,0,0)
    redF32 = RGB{Float32}(1,0,0)
    red24 = reinterpret(RGB24, 0x00ff0000)
    red32 = reinterpret(ARGB32, 0xffff0000)
    @testset "type check: RGB, RGB{$T}" for T in (Float64, Float32, N0f8)
        c = RGB(one(T), zero(T), zero(T))
        @test eltype(c) == T
        @test convert(RGB{Float64}, c) === redF64
        @test convert(RGB{T}, redF64) == c
    end
    @test RGB(1,0,0) == redF64
    @test RGB(UInt8(1), 0, 0) == redF64
    @test RGB(UInt8(1), UInt8(0), UInt8(0)) == redF64 # != colorant"#010000"
    @test convert(RGB, red24) == redF64

    @testset "type check: RGB24-->$Cto" for Cto in parametric3
        @test typeof(convert(Cto, red24)) == Cto{eltype_default(Cto)}
        @test typeof(convert(Cto{Float64}, red24)) == Cto{Float64}
    end

    @testset "type check: C{Float}-->C{Float}" for Cfrom in parametric3, Tfrom in (Float32, Float64)
        c = convert(Cfrom{Tfrom}, redF64)
        @test typeof(c) == Cfrom{Tfrom}
        @testset "$Cfrom{$Tfrom}-->$Cto" for Cto in parametric3
            c1 = convert(Cto, c)
            @test eltype(c1) == Tfrom
        end
        @testset "$Cfrom{$Tfrom}-->$Cto{$Tto}" for Cto in parametric3, Tto in (Float32, Float64)
            c2 = convert(Cto{Tto}, c)
            @test typeof(c2) == Cto{Tto}
        end
    end

    normed_types = (N0f8, N6f10, N4f12, N2f14, N0f16)
    # Test conversion from Normed types
    @testset "type check: C{Normed}-->C{Float}" for Cfrom in fractional_types, Tfrom in normed_types
        c = convert(Cfrom{Tfrom}, redF64)
        @test typeof(c) == Cfrom{Tfrom}
        @testset "$Cfrom{$Tfrom}-->$Cto" for Cto in parametric3
            eltype_default(Cto) <: FixedPoint && continue
            c1 = convert(Cto, c)
            @test eltype(c1) == eltype_default(Cto)
        end
        @testset "$Cfrom{$Tfrom}-->$Cto{$Tto}" for Cto in parametric3, Tto in (Float32, Float64)
            c2 = convert(Cto{Tto}, c)
            @test typeof(c2) == Cto{Tto}
        end
    end

    # Test conversion to Normed types
    @testset "type check: C{Float}-->C{Normed}" for Cfrom in parametric3, Tfrom in (Float32, Float64)
        c = convert(Cfrom{Tfrom}, redF64)
        @test typeof(c) == Cfrom{Tfrom}
        @testset "$Cfrom{$Tfrom}-->$Cto{$Tto}" for Cto in fractional_types, Tto in normed_types
            c2 = convert(Cto{Tto}, c)
            @test typeof(c2) == Cto{Tto}
        end
    end

    ac = RGBA(redF64)

    @test convert(RGB, ac) == RGB(1,0,0)
    @test convert(RGB{N0f8}, ac) == RGB{N0f8}(1,0,0)
    @test convert(RGBA{N0f8}, ac) == RGBA{N0f8}(1,0,0,1)
    @test convert(HSVA, ac) == HSVA{Float64}(convert(HSV, redF64), 1.0)
    @test convert(HSVA{Float32}, ac) == HSVA{Float32}(convert(HSV{Float32}, redF64), 1.0f0)
    @test convert(RGBA, redF64) == ac

    @test convert(ARGB32, ac) == reinterpret(ARGB32, 0xffff0000)
    @test convert(ARGB32, RGB(1,0,0)) == reinterpret(ARGB32, 0xffff0000)
    ac32 = convert(ARGB32, RGB(1,0,0), 0.5)
    @test ac32 == reinterpret(ARGB32, 0x80ff0000)
    @test convert(ARGB32, ac32) === ac32
    @test convert(ARGB32, BGRA(1,0,0,0.5)) == reinterpret(ARGB32, 0x80ff0000)
    @test reinterpret(UInt32, convert(ARGB32, ac)) == 0xffff0000
    @test convert(RGB24, RGB(r8(0xff),r8(0x00),r8(0x00))) == reinterpret(RGB24, 0x00ff0000)
    @test reinterpret(UInt32, convert(RGB24, RGB(r8(0xff),r8(0x00),r8(0x00)))) == 0x00ff0000
    redhsv = convert(HSV, redF64)
    @test convert(RGB24, red24) === red24
    @test convert(RGB24, redhsv) == reinterpret(RGB24, 0x00ff0000)
    @test_throws ArgumentError convert(RGB24,  RGB(0, 1.1, 0))
    @test_throws ArgumentError convert(ARGB32, RGBA(0, 1.1, 0, 0.8))
    @test_throws ArgumentError convert(ARGB32, RGBA(0, 0.8, 0, 1.1))

    @test convert(RGB{N0f8}, red24) == RGB{N0f8}(1,0,0)
    @test convert(RGBA{N0f8}, red32) == RGBA{N0f8}(1,0,0,1)
    @test convert(HSVA{Float64}, red32) == HSVA{Float64}(0, 1, 1, 1)

    # whitepoint conversions
    @test isa(convert(XYZ, convert(Lab, redF64), Colors.WP_DEFAULT), XYZ{Float64})
    @test isa(convert(XYZ{Float32}, convert(Lab, redF64), Colors.WP_DEFAULT), XYZ{Float32})
    @test isa(convert(XYZ, convert(Luv, redF64), Colors.WP_DEFAULT), XYZ{Float64})
    @test isa(convert(XYZ{Float32}, convert(Luv, redF64), Colors.WP_DEFAULT), XYZ{Float32})
    @test isa(convert(Lab, convert(XYZ, redF64), Colors.WP_DEFAULT), Lab{Float64})
    @test isa(convert(Lab{Float32}, convert(XYZ, redF64), Colors.WP_DEFAULT), Lab{Float32})
    @test isa(convert(Luv, convert(XYZ, redF64), Colors.WP_DEFAULT), Luv{Float64})
    @test isa(convert(Luv{Float32}, convert(XYZ, redF64), Colors.WP_DEFAULT), Luv{Float32})

    #59
    @test Colors.xyz_to_uv(XYZ{Float64}(0.0, 0.0, 0.0)) === (0.0, 0.0)
    @test Colors.xyz_to_uv(XYZ{Float64}(0.0, 1.0, 0.0)) === (0.0, 0.6)
    @test Colors.xyz_to_uv(XYZ{Float64}(0.0, 1.0, 1.0)) === (0.0, 0.5)
    @test Colors.xyz_to_uv(XYZ{Float64}(1.0, 0.0, 1.0)) === (1.0, 0.0)
    @test Colors.xyz_to_uv(XYZ{Float64}(1.0, 0.0, 0.0)) === (4.0, 0.0)

    # ColorTypes.jl issue #40
    @test convert(HSL, RGB{N0f8}(0.678, 0.847, 0.902)) ≈ HSL{Float32}(194.73685f0,0.5327105f0,0.7901961f0) atol=100eps(Float32)

    # YIQ
    @test convert(YIQ, RGB(1,0,0)) == YIQ{Float32}(0.299, 0.595716, 0.211456)
    @test convert(YIQ, RGB(0,1,0)) == YIQ{Float32}(0.587, -0.274453, -0.522591)
    @test convert(YIQ, RGB(0,0,1)) == YIQ{Float32}(0.114, -0.321263, 0.311135)
    @test convert(RGB, YIQ(1.0,0.0,0.0)) == RGB(1,1,1)
    v = 0.5957
    @test convert(RGB, YIQ(0.0,1.0,0.0)) == RGB(0.9563*v,0,0)
    v = -0.5226
    @test convert(RGB, YIQ(0.0,0.0,-1.0)) == RGB(0,-0.6474*v,0)

    # Gray
    @test convert(Gray{N0f8}, Gray{N0f8}(0.1)) == Gray{N0f8}(0.1)
    @test convert(Gray{N0f8}, Gray(0.1))     == Gray{N0f8}(0.1)
    @test convert(Gray{N0f8}, Gray24(0.1))   == Gray{N0f8}(0.1)
    @test convert(Gray24, Gray{N0f8}(0.1))   == Gray24(0.1)

    @test convert(RGB{N0f8}, Gray{N0f8}(0.1)) == RGB{N0f8}(0.1,0.1,0.1)
    @test convert(RGB{N0f8}, Gray24(0.1))   == RGB{N0f8}(0.1,0.1,0.1)

    grayN0f16 = Gray{N0f16}(0.8)
    @test convert(RGB, grayN0f16) == RGB{N0f16}(0.8,0.8,0.8)
    @test convert(RGB{Float32}, grayN0f16) == RGB{Float32}(0.8,0.8,0.8)

    @testset "$C-->Gray" for C in parametric3
        c = convert(C, RGB(1,1,1))
        g1 = convert(Gray, c)
        @test isa(g1, Gray)
        @test gray(g1) ≈ 1 atol=0.01
        g2 = convert(Gray{Float64}, c)
        @test typeof(g2) == Gray{Float64}
        @test gray(g2) ≈ 1 atol=0.01
    end
    # Issue #377
    @test convert(Gray, RGB24(1,0,0)) === convert(Gray, RGB(1,0,0)) === Gray{N0f8}(0.298)
    @test convert(Gray24, RGB(1,0,0)) === Gray24(0.298)
    # Check for roundoff error
    for N in (N0f8, N0f16, N0f32)
        @test convert(Gray{N}, RGB{N}(1,1,1)) === Gray{N}(1)
    end
    @test gray(convert(Gray{N0f64}, RGB{N0f64}(1,1,1))) ≈ 1.0

    # Images issue #382
    @test convert(Gray, RGBA(1,1,1,1)) == Gray(N0f8(1))

    # https://github.com/timholy/Images.jl/pull/445#issuecomment-189866806
    @test convert(Gray, RGB{N0f8}(0.145,0.145,0.145)) == Gray{N0f8}(0.145)

    # More AbstractRGB
    r4 = RGBX(1,0,0)
    @test convert(RGB, r4) == RGB(1,0,0)
    @test convert(RGB{N0f8}, r4) == RGB{N0f8}(1,0,0)
    @test convert(RGBX{N0f8}, r4) == RGBX{N0f8}(1,0,0)
    @test convert(RGBX{Float32}, r4) == RGBX{Float32}(1,0,0)
    @test convert(BGR{Float32}, r4) == BGR{Float32}(1,0,0)

    # Issue #257
    c = RGB{Float16}(0.9473,0.962,0.9766)
    convert(HSI, c)
    hsi = convert(HSI, c)
    @test hsi.i > 0.96 && hsi.h ≈ 210

    # {HSV, HSL, HSI} --> RGB (issue #379)
    @testset "HSx --> RGB" begin
        @test convert(RGB, HSV{Float32}( 780, 1, 1)) === RGB{Float32}(1,1,0)
        @test convert(RGB, HSV{Float32}(   0, 1, 1)) === RGB{Float32}(1,0,0)
        @test convert(RGB, HSV{Float32}(-780, 1, 1)) === RGB{Float32}(1,0,1)
        @test convert(RGB, HSV{Float32}(30,  2, .5)) === RGB{Float32}(.5,.25,0)
        @test convert(RGB, HSV{Float32}(30, .5, -1)) === RGB{Float32}(0,0,0)
        @test convert(RGB{Float64}, HSV{BigFloat}(-360120, 2, 1)) === RGB{Float64}(0,0,1)

        @test convert(RGB, HSL{Float32}( 780, 1, .5)) === RGB{Float32}(1,1,0)
        @test convert(RGB, HSL{Float32}(   0, 1, .5)) === RGB{Float32}(1,0,0)
        @test convert(RGB, HSL{Float32}(-780, 1, .5)) === RGB{Float32}(1,0,1)
        @test convert(RGB, HSL{Float32}(30,  2, .25)) === RGB{Float32}(.5,.25,0)
        @test convert(RGB, HSL{Float32}(30, .5,  -1)) === RGB{Float32}(0,0,0)
        @test convert(RGB{Float64}, HSL{BigFloat}(-360120, 2, .5)) === RGB{Float64}(0,0,1)

        @test convert(RGB, HSI{Float32}( 780, .5, .5)) ≈ RGB{Float32}(.625,.625,.25)
        @test convert(RGB, HSI{Float32}(   0, .5, .5)) ≈ RGB{Float32}(1,.25,.25)
        @test convert(RGB, HSI{Float32}(-780, .5, .5)) ≈ RGB{Float32}(.625,.25,.625)
        @test convert(RGB, HSI{Float32}(30,  2, .25)) ≈ RGB{Float32}(.5,.25,0)
        @test convert(RGB, HSI{Float32}(30, .5,  -1)) ≈ RGB{Float32}(0,0,0)
        @test convert(RGB{Float64}, HSI{BigFloat}(-360120, .5, .5)) ≈ RGB{Float64}(.25,.25,1)
    end

    @testset "custom types" begin
        # issue #465
        @test_throws ErrorException convert(RGB,   C3{Float32}(1, 2, 3))
        @test_throws ErrorException convert(RGB24, C3{Float64}(1, 2, 3))
        @test_throws ErrorException convert(Lab,   C3{Float16}(1, 2, 3))

        @test @inferred(convert(RGB,   C4{Float32}(0.1, 0.2, 0.3, 0.4))) ≈ RGB{Float32}(0.5, 0.6, 0.7)
        @test @inferred(convert(RGB24, C4{Float64}(0.1, 0.2, 0.3, 0.4))) ≈ RGB24(0.5, 0.6, 0.7)
        @test @inferred(convert(ALab,  C4{Float16}(0.1, 0.2, 0.3, 0.4))) ≈ ALab{Float16}(62.12, -2.86, -16.25, 1.0)

        @test @inferred(convert(C4, RGB{Float32}(0.5, 0.6, 0.7))) ≈ C4{Float32}(0.0, 0.1, 0.2, 0.5)
        @test @inferred(convert(C4{N0f8}, ARGB32(0.5, 0.6, 0.7))) ≈ C4{N0f16}(0.0, 0.1, 0.2, 0.5)
        @test @inferred(convert(C4, Lab{Float16}(62.12, -2.86, -16.25))) ≈ C4{Float16}(0.0, 0.1, 0.2, 0.5)
    end

    # Test accuracy of conversion
    include("test_conversions.jl")

    # Since `colordiff`(e.g. `DE_2000`) involves a color space conversion, it is
    # not suitable for evaluating the conversion itself. On the other hand,
    # since the tolerance varies from component to component, a homogeneous
    # error evaluation function (e.g. a simple sum of differences) is also not
    # appropriate. Therefore, a series of `diffnorm`, which returns the
    # normalized Euclidean distance, is defined as follows. They are just for
    # testing purposes as the cyclicity of hue is ignored.
    sqd(a, b, s=1.0) = ((float(a) - float(b))/s)^2
    function diffnorm(a::T, b::T) where {T<:Color3} # RGB,XYZ,xyY,LMS
        sqrt(sqd(comp1(a), comp1(b)) + sqd(comp2(a), comp2(b)) + sqd(comp3(a),comp3(b)))/sqrt(3)
    end
    function diffnorm(a::T, b::T) where {T<:Union{HSV,HSL,HSI}}
        sqrt(sqd(a.h, b.h, 360) + sqd(a.s, b.s) + sqd(comp3(a), comp3(b)))/sqrt(3)
    end
    function diffnorm(a::T, b::T) where {T<:Union{Lab,Luv}}
        sqrt(sqd(a.l, b.l, 100) + sqd(comp2(a), comp2(b), 200) + sqd(comp3(a), comp3(b), 200))/sqrt(3)
    end
    function diffnorm(a::T, b::T) where {T<:Union{LCHab,LCHuv}}
        sqrt(sqd(a.l, b.l, 100) + sqd(a.c, b.c, 100) + sqd(a.h, b.h, 360))/sqrt(3)
    end
    function diffnorm(a::T, b::T) where {T<:Union{DIN99,DIN99d,DIN99o}} # csconv has no DIN99 case
        sqrt(sqd(a.l, b.l, 100) + sqd(a.a, b.a, 100) + sqd(a.b, b.b, 100))/sqrt(3)
    end
    function diffnorm(a::T, b::T) where {T<:YIQ}
        sqrt(sqd(a.y, b.y) + sqd(a.i, b.i, 1.2) + sqd(a.q, b.q, 1.2))/sqrt(3)
    end
    function diffnorm(a::T, b::T) where {T<:YCbCr}
        sqrt(sqd(a.y, b.y, 219) + sqd(a.cb, b.cb, 224) + sqd(a.cr, b.cr, 224))/sqrt(3)
    end

    for C in ColorTypes.parametric3
        y = convert(C, RGB(1.0,1.0,0.0))
        b = convert(C, RGB(0.1,0.1,0.2))
        diffnorm(y, b) < 0.5 && @warn("`diffnorm` for $C may be broken")
    end

    function convcompare(from, to, tol; showfailure::Bool=false)
        errmax = 0.0
        for i = 1:length(from)
            t = to[i]
            f = convert(typeof(t), from[i])
            diff = diffnorm(t, f)
            if showfailure && diff>tol
                original = from[i]
                @show original f t diff
            end
            errmax = max(errmax, diff)
        end
        errmax > tol && @warn("Error on conversions from $(eltype(from)) to $(eltype(to)), relative error = $errmax")
        errmax <= tol
    end

    @testset "accuracy test: $Cfrom-->$Cto" for Cto in keys(csconv), Cfrom in keys(csconv)
        f = csconv[Cfrom]
        t = csconv[Cto]
        @test convcompare(f, t, 2e-3, showfailure=false)
    end
end


text:
# ------------------------------------------------------------------------------------------
# # Introduction to DataFrames
# **[Bogumił Kamiński](http://bogumilkaminski.pl/about/), May 23, 2018**
#
# Let's get started by loading the `DataFrames` package.
# ------------------------------------------------------------------------------------------

using DataFrames

# ------------------------------------------------------------------------------------------
# ## Constructors and conversion
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ### Constructors
#
# In this section, you'll see many ways to create a `DataFrame` using the `DataFrame()`
# constructor.
#
# First, we could create an empty DataFrame,
# ------------------------------------------------------------------------------------------

DataFrame() # empty DataFrame

# ------------------------------------------------------------------------------------------
# Or we could call the constructor using keyword arguments to add columns to the
# `DataFrame`.
# ------------------------------------------------------------------------------------------

DataFrame(A=1:3, B=rand(3), C=randstring.([3,3,3]))

# ------------------------------------------------------------------------------------------
# We can create a `DataFrame` from a dictionary, in which case keys from the dictionary will
# be sorted to create the `DataFrame` columns.
# ------------------------------------------------------------------------------------------

x = Dict("A" => [1,2], "B" => [true, false], "C" => ['a', 'b'])
DataFrame(x)

# ------------------------------------------------------------------------------------------
# Rather than explicitly creating a dictionary first, as above, we could pass `DataFrame`
# arguments with the syntax of dictionary key-value pairs.
#
# Note that in this case, we use symbols to denote the column names and arguments are not
# sorted. For example, `:A`, the symbol, produces `A`, the name of the first column here:
# ------------------------------------------------------------------------------------------

DataFrame(:A => [1,2], :B => [true, false], :C => ['a', 'b'])

# ------------------------------------------------------------------------------------------
# Here we create a `DataFrame` from a vector of vectors, and each vector becomes a column.
# ------------------------------------------------------------------------------------------

DataFrame([rand(3) for i in 1:3])

# ------------------------------------------------------------------------------------------
#  For now we can construct a single `DataFrame` from a `Vector` of atoms, creating a
# `DataFrame` with a single row. In future releases of DataFrames.jl, this will throw an
# error.
# ------------------------------------------------------------------------------------------

DataFrame(rand(3))

# ------------------------------------------------------------------------------------------
# Instead use a transposed vector if you have a vector of atoms (in this way you effectively
# pass a two dimensional array to the constructor which is supported).
# ------------------------------------------------------------------------------------------

DataFrame(transpose([1, 2, 3]))

# ------------------------------------------------------------------------------------------
# Pass a second argument to give the columns names.
# ------------------------------------------------------------------------------------------

DataFrame([1:3, 4:6, 7:9], [:A, :B, :C])

# ------------------------------------------------------------------------------------------
# Here we create a `DataFrame` from a matrix,
# ------------------------------------------------------------------------------------------

DataFrame(rand(3,4))

# ------------------------------------------------------------------------------------------
# and here we do the same but also pass column names.
# ------------------------------------------------------------------------------------------

DataFrame(rand(3,4), Symbol.('a':'d'))

# ------------------------------------------------------------------------------------------
# We can also construct an uninitialized DataFrame.
#
# Here we pass column types, names and number of rows; we get `missing` in column :C because
# `Any >: Missing`.
# ------------------------------------------------------------------------------------------

DataFrame([Int, Float64, Any], [:A, :B, :C], 1)

# ------------------------------------------------------------------------------------------
# Here we create a `DataFrame`, but column `:C` is #undef and Jupyter has problem with
# displaying it. (This works OK at the REPL.)
#
# This will be fixed in next release of DataFrames!
# ------------------------------------------------------------------------------------------

DataFrame([Int, Float64, String], [:A, :B, :C], 1)

# ------------------------------------------------------------------------------------------
# To initialize a `DataFrame` with column names, but no rows use
# ------------------------------------------------------------------------------------------

DataFrame([Int, Float64, String], [:A, :B, :C], 0)

# ------------------------------------------------------------------------------------------
# This syntax gives us a quick way to create homogenous `DataFrame`.
# ------------------------------------------------------------------------------------------

DataFrame(Int, 3, 5)

# ------------------------------------------------------------------------------------------
# This example is similar, but has nonhomogenous columns.
# ------------------------------------------------------------------------------------------

DataFrame([Int, Float64], 4)

# ------------------------------------------------------------------------------------------
# Finally, we can create a `DataFrame` by copying an existing `DataFrame`.
#
# Note that `copy` creates a shallow copy.
# ------------------------------------------------------------------------------------------

y = DataFrame(x)
z = copy(x)
(x === y), (x === z), isequal(x, z)

# ------------------------------------------------------------------------------------------
# ### Conversion to a matrix
#
# Let's start by creating a `DataFrame` with two rows and two columns.
# ------------------------------------------------------------------------------------------

x = DataFrame(x=1:2, y=["A", "B"])

# ------------------------------------------------------------------------------------------
# We can create a matrix by passing this `DataFrame` to `Matrix`.
# ------------------------------------------------------------------------------------------

Matrix(x)

# ------------------------------------------------------------------------------------------
# This would work even if the `DataFrame` had some `missing`s:
# ------------------------------------------------------------------------------------------

x = DataFrame(x=1:2, y=[missing,"B"])

Matrix(x)

# ------------------------------------------------------------------------------------------
# In the two previous matrix examples, Julia created matrices with elements of type `Any`.
# We can see more clearly that the type of matrix is inferred when we pass, for example, a
# `DataFrame` of integers to `Matrix`, creating a 2D `Array` of `Int64`s:
# ------------------------------------------------------------------------------------------

x = DataFrame(x=1:2, y=3:4)

Matrix(x)

# ------------------------------------------------------------------------------------------
# In this next example, Julia correctly identifies that `Union` is needed to express the
# type of the resulting `Matrix` (which contains `missing`s).
# ------------------------------------------------------------------------------------------

x = DataFrame(x=1:2, y=[missing,4])

Matrix(x)

# ------------------------------------------------------------------------------------------
# Note that we can't force a conversion of `missing` values to `Int`s!
# ------------------------------------------------------------------------------------------

Matrix{Int}(x)

# ------------------------------------------------------------------------------------------
# ### Handling of duplicate column names
#
# We can pass the `makeunique` keyword argument to allow passing duplicate names (they get
# deduplicated)
# ------------------------------------------------------------------------------------------

df = DataFrame(:a=>1, :a=>2, :a_1=>3; makeunique=true)

# ------------------------------------------------------------------------------------------
# Otherwise, duplicates will not be allowed in the future.
# ------------------------------------------------------------------------------------------

df = DataFrame(:a=>1, :a=>2, :a_1=>3)

# ------------------------------------------------------------------------------------------
# A constructor that is passed column names as keyword arguments is a corner case.
# You cannot pass `makeunique` to allow duplicates here.
# ------------------------------------------------------------------------------------------

df = DataFrame(a=1, a=2, makeunique=true)

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
using Kalman, Plots

import Base:+
+(p::Mark, i::CartesianIndex{2}) = p + Mark(i.I...)
+(i::CartesianIndex{2}, p::Mark) = p + i

const Image = Matrix{Gray{N0f16}}

const speed = 0.1719606156256064 # mean downwards speed along the rows of the image in mm per hour
# const ρ = (0.8, 0.5) # how much do you trust each speed. This has been calculated from a bunch of datasets
# const ρ = 0.8 # how much do you trust each speed. This has been calculated from a bunch of datasets
# const Qρ = (0.38629498696108994*(1 - ρ[1]^2), 0.2584574519006643*(1 - ρ[2]^2))
# const Qρ = 0.32*(1 - ρ^2)

const weight_disk = disk(weight_radius)
const intensity_disk = disk(intensity_radius)

function image_feedback(img::Image, p::Mark)
    ind = mark2cart(p)
    μ = mean(img[i + ind] for i in weight_disk)
    S = sum(max(0, img[i + ind] - μ) for i in weight_disk)
    sum(max(0, img[i + ind] - μ)*(i + p) for i in weight_disk)/S
end

function initiatemodel(vrow::Float64)
    p0 = SVector(NaN, NaN, vrow, 0.0)
    P0 = 10*@SMatrix eye(4)
    # A = SMatrix{4, 4, Float64}([1 0 1 0
    #                             0 1 0 1
    #                             0 0 ρ 0
    #                             0 0 0 ρ])
    A = SMatrix{4, 4, Float64}([1 0 1 0
                                0 1 0 1
                                0 0 1 0
                                0 0 0 1])
    # b = SVector(0, 0, (1-ρ)*vrow, 0)
    b = SVector(0., 0, 0, 0)
    # Q = 1*SMatrix{4,4, Float64}(diagm([1, 1, Qρ, Qρ]))
    Q = 1*SMatrix{4,4, Float64}(diagm([1, 1, 5, 5]))
    y = SVector(NaN, NaN)
    C = @SMatrix eye(2, 4)
    R = 2*@SMatrix eye(2)
    LinearHomogSystem(p0, P0, A, b, Q, y, C, R)
end

function getintensity(img::Image, p::Mark)
    ind = CartesianIndex(round.(Int, p)...)
    sum(img[i + ind] for i in intensity_disk)
end

struct Track
    lengths::Vector{Float64}
    coordinates::Vector{Mark}
    times::Vector{Float64}
    intensities::Vector{Vector{Float64}}
    id::Int
    color::RGB{Float64}
    function Track(p::Mark, st::CalibStage, id::Int, color::RGB{Float64})
        lengths = Float64[]
        coordinates = Mark[]
        times = Float64[]
        intensities = Vector{Float64}[]
        ntl = length(st.timelapse)
        sizehint!(lengths, ntl)
        sizehint!(coordinates, ntl)
        sizehint!(times, ntl)
        sizehint!(intensities, ntl)
        push!(lengths, 0.0)
        push!(coordinates, p)
        push!(times, st.timelapse[1].time)
        img = load(st.timelapse[1].path)
        I = getintensity(img, p)
        push!(intensities, [I])
        new(lengths, coordinates, times, intensities, id, color)
    end
end

function updatetrack(t::Track, img::Image, p::Mark, time::Float64)
    push!(t.lengths, t.lengths[end] + norm(t.coordinates[end] - p))
    push!(t.coordinates, p)
    push!(t.times, time)
    I = [getintensity(img, i) for i in t.coordinates]
    push!(t.intensities, I)
end

mutable struct Root
    model::LinearHomogSystem
    x::SVector{4, Float64}
    P::SMatrix{4, 4, Float64}
    grow::Bool
    function Root(p::Mark, vrow::Float64)
        model = initiatemodel(vrow)
        x = SVector(p..., model.x0[3:4]...)
        P = model.P0
        new(model, x, P, true)
    end
end

function correctpredict!(r::Root, img::Image, t1::Float64, t2::Float64)
    x, Ppred, A = Kalman.predict!(t1, r.x, r.P, t2, r.model)
    if outside(x[1:2])
        r.grow = false
        return nothing
    end
    y = image_feedback(img, Mark(x[1:2]...))
    _, obs, C, R = Kalman.observe!(t1, x, r.P, t2, y, r.model)
    x, P, yres, S, K = Kalman.correct!(x, Ppred, obs, C, R, r.model)
    if outside(r.x[1:2])
        r.grow = false
        return nothing
    end
    r.x = x
    r.P = P
    return nothing
end

function trackroot(st::CalibStage, startpoints::Vector{Mark})
    isempty(startpoints) && return Track[]
    vrow = speed/st.speed
    roots = [Root(p, vrow) for p in startpoints]
    n = length(startpoints)
    inds = 1:n
    m = 5
    C(g::ColorGradient) = RGB[g[z] for z=linspace(0,1,m)]
    colors = distinguishable_colors(n+m, C(cgrad(:inferno)))[m+1:end]
    tracks = [Track(p, st, i, c) for (i, p, c) in zip(inds, startpoints, colors)]
    for (tl1, tl2) in zip(st.timelapse[2:end], [st.timelapse[3:end]; st.timelapse[end]])
        img = load(tl1.path)
        for (r, t) in zip(roots, tracks)
            if r.grow
                correctpredict!(r, img, tl1.time, tl2.time)
                updatetrack(t, img, Mark(r.x[1:2]), tl2.time)
            end
        end
    end
    return tracks
end

# end # module


text:
export @snoopc

using Serialization

"""
```
@snoopc "compiledata.csv" begin
    # Commands to execute, in a new process
end
```
causes the julia compiler to log all functions compiled in the course
of executing the commands to the file "compiledata.csv". This file
can be used for the input to `SnoopCompile.read`.
"""
macro snoopc(flags, filename, commands)
    return :(snoopc($(esc(flags)), $(esc(filename)), $(QuoteNode(commands))))
end
macro snoopc(filename, commands)
    return :(snoopc(String[], $(esc(filename)), $(QuoteNode(commands))))
end

function snoopc(flags, filename, commands)
    println("Launching new julia process to run commands...")
    # addprocs will run the unmodified version of julia, so we
    # launch it as a command.
    code_object = """
            using Serialization
            while !eof(stdin)
                Core.eval(Main, deserialize(stdin))
            end
            """
    process = open(`$(Base.julia_cmd()) $flags --eval $code_object`, stdout, write=true)
    serialize(process, quote
        let io = open($filename, "w")
            ccall(:jl_dump_compiles, Nothing, (Ptr{Nothing},), io.handle)
            try
                $commands
            finally
                ccall(:jl_dump_compiles, Nothing, (Ptr{Nothing},), C_NULL)
                close(io)
            end
        end
        exit()
    end)
    wait(process)
    println("done.")
    nothing
end


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

mutable struct ApplicationGatewaySslPredefinedPolicyPropertiesFormat <: SwaggerModel
    cipherSuites::Any # spec type: Union{ Nothing, Vector{CipherSuitesEnum} } # spec name: cipherSuites
    minProtocolVersion::Any # spec type: Union{ Nothing, ProtocolsEnum } # spec name: minProtocolVersion

    function ApplicationGatewaySslPredefinedPolicyPropertiesFormat(;cipherSuites=nothing, minProtocolVersion=nothing)
        o = new()
        validate_property(ApplicationGatewaySslPredefinedPolicyPropertiesFormat, Symbol("cipherSuites"), cipherSuites)
        setfield!(o, Symbol("cipherSuites"), cipherSuites)
        validate_property(ApplicationGatewaySslPredefinedPolicyPropertiesFormat, Symbol("minProtocolVersion"), minProtocolVersion)
        setfield!(o, Symbol("minProtocolVersion"), minProtocolVersion)
        o
    end
end # type ApplicationGatewaySslPredefinedPolicyPropertiesFormat

const _property_map_ApplicationGatewaySslPredefinedPolicyPropertiesFormat = Dict{Symbol,Symbol}(Symbol("cipherSuites")=>Symbol("cipherSuites"), Symbol("minProtocolVersion")=>Symbol("minProtocolVersion"))
const _property_types_ApplicationGatewaySslPredefinedPolicyPropertiesFormat = Dict{Symbol,String}(Symbol("cipherSuites")=>"Vector{CipherSuitesEnum}", Symbol("minProtocolVersion")=>"ProtocolsEnum")
Base.propertynames(::Type{ ApplicationGatewaySslPredefinedPolicyPropertiesFormat }) = collect(keys(_property_map_ApplicationGatewaySslPredefinedPolicyPropertiesFormat))
Swagger.property_type(::Type{ ApplicationGatewaySslPredefinedPolicyPropertiesFormat }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ApplicationGatewaySslPredefinedPolicyPropertiesFormat[name]))}
Swagger.field_name(::Type{ ApplicationGatewaySslPredefinedPolicyPropertiesFormat }, property_name::Symbol) =  _property_map_ApplicationGatewaySslPredefinedPolicyPropertiesFormat[property_name]

function check_required(o::ApplicationGatewaySslPredefinedPolicyPropertiesFormat)
    true
end

function validate_property(::Type{ ApplicationGatewaySslPredefinedPolicyPropertiesFormat }, name::Symbol, val)
end


text:
using MeshCatMechanisms
using LinearAlgebra
using RigidBodyDynamics
using Plots
using GeometryTypes
using CoordinateTransformations
import TrajectoryOptimization: hold_trajectory, Trajectory, total_time
import RigidBodyDynamics: transform
include("N_plots.jl")
include("../kuka_visualizer.jl")
model, obj = Dynamics.kuka
n,m = model.n, model.m
nn = m   # Number of positions

function plot_sphere(vis::MechanismVisualizer,frame::CartesianFrame3D,center,radius,mat,name="")
    geom = HyperSphere(Point3f0(center), convert(Float32,radius))
    setelement!(vis,frame,geom,mat,name)
end

function plot_cylinder(vis::MechanismVisualizer,frame::CartesianFrame3D,c1,c2,radius,mat,name="")
    geom = Cylinder(Point3f0(c1),Point3f0(c2),convert(Float32,radius))
    setelement!(vis,frame,geom,mat,name)
end

# Create Mechanism
kuka = parse_urdf(Dynamics.urdf_kuka,remove_fixed_tree_joints=false)
kuka_visuals = URDFVisuals(Dynamics.urdf_kuka)
state = MechanismState(kuka)
world = root_frame(kuka)

# Create Visualizer
vis = Visualizer()
mvis = MechanismVisualizer(kuka, kuka_visuals, vis[:base])
open(vis)

# Default solver options
opts = SolverOptions()
opts.cost_tolerance = 1e-6
opts.constraint_tolerance = 1e-5
dircol_options = Dict("tol"=>opts.cost_tolerance,"constr_viol_tol"=>opts.constraint_tolerance)

##########################################################
#            Full body Obstacle Avoidance                #
##########################################################

# Collision Bubbles
function kuka_points(kuka,plot=false)
    bodies = [3,4,5,6]
    radii = [0.1,0.12,0.09,0.09]
    ee_body,ee_point = Dynamics.get_kuka_ee(kuka)
    ee_radii = 0.05

    points = Point3D[]
    frames = CartesianFrame3D[]
    for (idx,radius) in zip(bodies,radii)
        body = findbody(kuka,"iiwa_link_$idx")
        frame = default_frame(body)
        point = Point3D(frame,0.,0.,0.)
        plot ? plot_sphere(mvis,frame,0,radius,body_collision,"body_$idx") : nothing
        # plot ? setelement!(vis,point,radius,"body_$idx") : nothing
        push!(points,point)
        push!(frames,frame)
    end
    frame = default_frame(ee_body)
    plot ? plot_sphere(mvis,frame,0,ee_radii,body_collision,"end_effector") : nothing
    # plot ? setelement!(vis,ee_point,ee_radii,"end_effector") : nothing
    push!(points,ee_point)
    push!(radii,ee_radii)
    return points, radii, frames
end

function calc_kuka_points(x::Vector{T},points) where T
    state = state_cache[T]
    set_configuration!(state,x[1:nn])
    [transform(state,point,world).v for point in points]
end

function generate_collision_constraint(kuka::Mechanism, circles, cylinders=[])
    # Specify points along the arm
    points, radii = kuka_points(kuka)
    num_points = length(points)
    nCircle = length(circles)
    nCylinder = length(cylinders)
    num_obstacles = nCircle + nCylinder

    function cI_obstacles(c,x,u)
        nn = length(u)

        # Get current world location of points along arm
        arm_points = calc_kuka_points(x[1:nn],points)

        C = reshape(view(c,1:num_points*num_obstacles),num_points,num_obstacles)
        for i = 1:nCircle
            c_obstacle = view(C,1:num_points,i)
            for (p,kuka_point) in enumerate(arm_points)
                c_obstacle[p] = sphere_constraint(circles[i],kuka_point,radii[p]+circles[i][4])
            end
        end
        for j = 1:nCylinder
            i = j + nCircle
            c_obstacle = view(C,1:num_points,i)
            for (p,kuka_point) in enumerate(arm_points)
                c_obstacle[p] = circle_constraint(cylinders[j],kuka_point,radii[p]+cylinders[j][3])
            end
        end
    end
end

function addcircles!(vis,circles)
    world = root_frame(kuka)
    for (i,circle) in enumerate(circles)
        p = Point3D(world,collect(circle[1:3]))
        setelement!(vis,p,circle[4],"obs$i")
    end
end

function addcylinders!(vis,cylinders,height=1.5)
    for (i,cyl) in enumerate(cylinders)
        plot_cylinder(vis,world,[cyl[1],cyl[2],0],[cyl[1],cyl[2],height],cyl[3],blue_,"cyl_$i")
    end
end
ee_fun = Dynamics.get_kuka_ee_postition_fun(kuka)

# Add obstacles
d = 0.25
circles2 = Vector{Float64}[]
push!(circles2,[d,0.0,1.2,0.2])
push!(circles2,[0,-d,0.4,0.15])
push!(circles2,[0,-d,1.2,0.15])
addcircles!(mvis,circles2)

cylinders = [[d,-d,0.08],[d,d,0.08],[-d,-d,0.08]]
addcylinders!(mvis,cylinders)

points,radii,frames = kuka_points(kuka,false)

# Plot Goal
plot_sphere(mvis,world,xf_ee,0.04,green_,"goal")
plot_sphere(mvis,world,x0_ee,0.04,red_,"start")

# Generate constraint function
state_cache = StateCache(kuka)
num_obstacles = length(circles2)+length(cylinders)
c = zeros(length(points)*num_obstacles)
cI_arm_obstacles = generate_collision_constraint(kuka,circles2,cylinders)
cI_arm_obstacles(c,x0,zeros(m))

# Build the Objective
x0 = zeros(n)
x0[2] = pi/2
x0[3] = pi/2
x0[4] = pi/2
xf = zeros(n)
xf[1] = pi/2
xf[4] = pi/2
Q = Diagonal([ones(7); ones(7)*100])
Qf = 10.0*Diagonal(I,n)
R = 1e-2*Diagonal(I,m)
tf = 5.0
obj_uncon = LQRObjective(Q, R, Qf, tf, x0, xf)
xf_ee = ee_fun(xf)
x0_ee = ee_fun(x0)

set_configuration!(mvis, x0[1:7])
set_configuration!(mvis, xf[1:7])
costfun = LQRCost(Q,R,Qf,xf)
obj_obs_arm = ConstrainedObjective(costfun,tf,x0,xf,cI=cI_arm_obstacles,u_min=-80,u_max=80)

# Solve the Problem
solver = Solver(model, obj_obs_arm, N=41)
solver.opts.verbose = true
solver.opts.penalty_scaling = 50
solver.opts.penalty_initial = 0.01
solver.opts.cost_tolerance_intermediate = 1e-2
solver.opts.cost_tolerance = 1e-4
solver.opts.constraint_tolerance = 1e-4
solver.opts.bp_reg_initial = 0
solver.opts.square_root = true
solver.opts.use_nesterov = false
solver.opts.iterations_outerloop = 50
solver.opts.iterations = 300
solver.opts.penalty_max = 1e8
solver.opts.outer_loop_update_type = :default
U0_hold = hold_trajectory(solver, kuka, x0[1:7])
res_obs, stats_obs = solve(solver,U0_hold)
stats_obs["iterations"]
stats_obs["max_mu_iteration"]
stats_obs["runtime"]
evals(solver,:f) / stats_obs["iterations"]
res = res_obs

# Visualize
set_configuration!(mvis, x0[1:7])
X_interp, U_interp = TrajectoryOptimization.interp_traj(201,5.,res_obs.X,res_obs.U)
animate_trajectory(mvis, X_interp, 0.01)
t = range(0,length=201,step=0.01)
q = [X_interp[1:7,i] for i = 1:201]
setanimation!(mvis,t,q)

stats0 = copy(stats_obs)
res0 = deepcopy(res_obs)
plot(stats0["c_max"],yscale=:log10)
plot!(stats_obs["c_max"])

plot(stats0["cost"],yscale=:log10)
plot!(stats_obs["cost"])

plot(res_obs.C[1:N-1],ylim=(0,0.1))

function plot_ghost_trajectory(urdf,K,N,α=0.5)
    mech = parse_urdf(urdf)
    visuals = URDFVisuals(urdf)
    state = MechanismState(mech)
    vis_el = visual_elements(mech, visuals)
    set_alpha!(vis_el,α)

    ks = round.(Int,range(1,N,length=K))
    for i = 1:K
        mvis2 = MechanismVisualizer(state, vis[Symbol("shadow$i")])
        MeshCatMechanisms._set_mechanism!(mvis2, vis_el)
        MeshCatMechanisms._render_state!(mvis2)
        set_configuration!(mvis2,res_obs.X[ks[i]][1:7])
    end
end

function set_alpha!(visuals::Vector{VisualElement}, α)
    for el in visuals
        c = el.color
        c_new = RGBA(red(c),green(c),blue(c),α)
        el.color = c_new
    end
end

plot_ghost_trajectory(Dynamics.urdf_kuka,10,solver.N,0.0)

# Solve Dircol
X0 = rollout(solver,U0_hold)
dircol_options = Dict("tol"=>1e-3,"constr_viol_tol"=>1e-2,"max_iter"=>20000,"mu_init"=>0.1)
X0 = line_trajectory(solver)
res_d, stats_d = solve_dircol(solver,to_array(res_obs.X),to_array(res_obs.U),options=dircol_options)
solver.opts.verbose = false
res_d, stats_d = solve_dircol(solver,X0,U0_hold,options=dircol_options)
evals(solver,:f) / stats_d["iterations"]

X_d, U_d = TrajectoryOptimization.interp_traj(201,5.,res_d.X,res_d.U)
animate_trajectory(mvis, X_d, 0.05)

settransform!(vis["/Cameras/default"], compose(Translation(-1.5, 0., 0.),LinearMap(RotX(0)*RotZ(0))))
settransform!(vis["/Cameras/default"], compose(Translation(0, 1.5, 0.),LinearMap(RotX(0)*RotZ(-pi/2))))

# new ghost plot
function plot_ghost_trajectory_2(urdf,traj_idx,α=0.5)
    mech = parse_urdf(urdf)
    visuals = URDFVisuals(urdf)
    state = MechanismState(mech)
    vis_el = visual_elements(mech, visuals)
    set_alpha!(vis_el,α)

    for i = 1:length(traj_idx)
        # if i == length(traj_idx)
        #     set_alpha!(vis_el,1.0)
        # end
        mvis2 = MechanismVisualizer(state, vis[Symbol("shadow$i")])
        MeshCatMechanisms._set_mechanism!(mvis2, vis_el)
        MeshCatMechanisms._render_state!(mvis2)
        set_configuration!(mvis2,res_obs.X[traj_idx[i]][1:7])
    end
end

traj_idx = [1;7;11;18;41]
set_configuration!(mvis, res_obs.X[traj_idx[end]][1:7])
plot_ghost_trajectory_2(Dynamics.urdf_kuka,traj_idx,0.65)
# note


text:
function smolyak_weights(y::AbstractArray{T,1},nodes::Union{Array{T,1},Array{T,2}},multi_index::Array{S,2}) where {T<:AbstractFloat,S<:Integer}

  interpolation_matrix = zeros(size(nodes,1),size(nodes,1))

  unique_multi_index = sort(unique(multi_index))
  unique_orders      = m_i(unique_multi_index).-1

  # Below we do the following things:

  #   Generate the polynomial terms for each order
  #   Generate the unique polynomial terms introduced at each higher order
  #   Combine the polynomial terms to construct a row of the interpolation matrix
  #   Iterate over the nodes, doing the above for steps at each iteration, to compute all rows of the interpolation matrix

  base_polynomials        = Array{Array{T,2},1}(undef,length(unique_orders))
  unique_base_polynomials = Array{Array{T,2},1}(undef,length(unique_orders))

  @inbounds for k = 1:size(nodes,1)

    # Construct the base polynomials

    for i = 1:length(unique_orders)
      base_polynomials[i] = chebyshev_polynomial(unique_orders[i],nodes[k,:])
    end

    # Compute the unique polynomial terms from the base polynomials

    for i = length(unique_orders):-1:2
      unique_base_polynomials[i] = base_polynomials[i][:,size(base_polynomials[i-1],2)+1:end]
    end
    unique_base_polynomials[1] = base_polynomials[1]

    # Construct a row of the interplation matrix

    l = 1
    @inbounds for j = 1:size(multi_index,1)
      new_polynomials = unique_base_polynomials[multi_index[j,1]][1,:]
      for i = 2:size(multi_index,2)
        new_polynomials = kron(new_polynomials,unique_base_polynomials[multi_index[j,i]][i,:])
      end
      m = length(new_polynomials)
      interpolation_matrix[k,l:l+m-1] = new_polynomials
      l += m
    end

  end

  weights = interpolation_matrix\y

  return weights

end

function smolyak_weights(y::AbstractArray{T,1},nodes::Union{Array{T,1},Array{T,2}},multi_index::Array{S,2},domain::Union{Array{T,1},Array{T,2}}) where {T<:AbstractFloat,S<:Integer}

  # Normalize nodes to the [-1.0 1.0] interval

  d = size(multi_index,2)
  nodes = copy(nodes)
  for i = 1:d
    nodes[:,i] = normalize_node(nodes[:,i],domain[:,i])
  end

  weights = smolyak_weights(y,nodes,multi_index)

  return weights

end

function smolyak_weights_threaded(y::AbstractArray{T,1},nodes::Union{Array{T,1},Array{T,2}},multi_index::Array{S,2}) where {T<:AbstractFloat,S<:Integer}

  interpolation_matrix = zeros(size(nodes,1),size(nodes,1))

  unique_multi_index = sort(unique(multi_index))
  unique_orders      = m_i(unique_multi_index).-1

  # Below we do the following things:

  #   Generate the polynomial terms for each order
  #   Generate the unique polynomial terms introduced at each higher order
  #   Combine the polynomial terms to construct a row of the interpolation matrix
  #   Iterate over the nodes, doing the above for steps at each iteration, to compute all rows of the interpolation matrix

  @inbounds @sync @qthreads for k = 1:size(nodes,1)

    base_polynomials        = Array{Array{T,2},1}(undef,length(unique_orders))
    unique_base_polynomials = Array{Array{T,2},1}(undef,length(unique_orders))

    # Construct the base polynomials

    for i = 1:length(unique_orders)
      base_polynomials[i] = chebyshev_polynomial(unique_orders[i],nodes[k,:])
    end

    # Compute the unique polynomial terms from the base polynomials

    for i = length(unique_orders):-1:2
      unique_base_polynomials[i] = base_polynomials[i][:,size(base_polynomials[i-1],2)+1:end]
    end
    unique_base_polynomials[1] = base_polynomials[1]

    # Construct a row of the interplation matrix

    l = 1
    @inbounds for j = 1:size(multi_index,1)
      new_polynomials = unique_base_polynomials[multi_index[j,1]][1,:]
      for i = 2:size(multi_index,2)
        new_polynomials = kron(new_polynomials,unique_base_polynomials[multi_index[j,i]][i,:])
      end
      m = length(new_polynomials)
      interpolation_matrix[k,l:l+m-1] = new_polynomials
      l += m
    end

  end

  weights = interpolation_matrix\y

  return weights

end

function smolyak_weights_threaded(y::AbstractArray{T,1},nodes::Union{Array{T,1},Array{T,2}},multi_index::Array{S,2},domain::Union{Array{T,1},Array{T,2}}) where {T<:AbstractFloat,S<:Integer}

  # Normalize nodes to the [-1.0 1.0] interval

  d = size(multi_index,2)
  nodes = copy(nodes)
  for i = 1:d
    nodes[:,i] = normalize_node(nodes[:,i],domain[:,i])
  end

  weights = smolyak_weights_threaded(y,nodes,multi_index)

  return weights

end

function smolyak_inverse_interpolation_matrix(nodes::Union{Array{T,1},Array{T,2}},multi_index::Array{S,2}) where {T<:AbstractFloat,S<:Integer}

  interpolation_matrix = zeros(size(nodes,1),size(nodes,1))

  unique_multi_index = sort(unique(multi_index))
  unique_orders      = m_i(unique_multi_index).-1

  # Below we do the following things:

  #   Generate the polynomial terms for each order
  #   Generate the unique polynomial terms introduced at each higher order
  #   Combine the polynomial terms to construct a row of the interpolation matrix
  #   Iterate over the nodes, doing the above for steps at each iteration, to compute all rows of the interpolation matrix

  base_polynomials        = Array{Array{T,2},1}(undef,length(unique_orders))
  unique_base_polynomials = Array{Array{T,2},1}(undef,length(unique_orders))

  @inbounds for k = 1:size(nodes,1)

    # Construct the base polynomials

    for i = 1:length(unique_orders)
      base_polynomials[i] = chebyshev_polynomial(unique_orders[i],nodes[k,:])
    end

    # Compute the unique polynomial terms from the base polynomials

    for i = length(unique_orders):-1:2
      unique_base_polynomials[i] = base_polynomials[i][:,size(base_polynomials[i-1],2)+1:end]
    end
    unique_base_polynomials[1] = base_polynomials[1]

    # Construct the first row of the interplation matrix

    l = 1
    @inbounds for j = 1:size(multi_index,1)
      new_polynomials = unique_base_polynomials[multi_index[j,1]][1,:]
      for i = 2:size(multi_index,2)
        new_polynomials = kron(new_polynomials,unique_base_polynomials[multi_index[j,i]][i,:])
      end
      m = length(new_polynomials)
      interpolation_matrix[k,l:l+m-1] = new_polynomials
      l += m
    end

  end

  inverse_interpolation_matrix = inv(interpolation_matrix)

  return inverse_interpolation_matrix

end

function smolyak_inverse_interpolation_matrix(nodes::Union{Array{T,1},Array{T,2}},multi_index::Array{S,2},domain::Union{Array{T,1},Array{T,2}}) where {T<:AbstractFloat,S<:Integer}

  # Normalize nodes to the [-1.0 1.0] interval

  d = size(multi_index,2)
  nodes = copy(nodes)
  for i = 1:d
    nodes[:,i] = normalize_node(nodes[:,i],domain[:,i])
  end

  inverse_interpolation_matrix = smolyak_inverse_interpolation_matrix(nodes,multi_index)

  return inverse_interpolation_matrix

end

function smolyak_inverse_interpolation_matrix_threaded(nodes::Union{Array{T,1},Array{T,2}},multi_index::Array{S,2}) where {T<:AbstractFloat,S<:Integer}

  interpolation_matrix = zeros(size(nodes,1),size(nodes,1))

  unique_multi_index = sort(unique(multi_index))
  unique_orders      = m_i(unique_multi_index).-1

  # Below we do the following things:

  #   Generate the polynomial terms for each order
  #   Generate the unique polynomial terms introduced at each higher order
  #   Combine the polynomial terms to construct a row of the interpolation matrix
  #   Iterate over the nodes, doing the above for steps at each iteration, to compute all rows of the interpolation matrix

  @inbounds @sync @qthreads for k = 1:size(nodes,1)

    base_polynomials        = Array{Array{T,2},1}(undef,length(unique_orders))
    unique_base_polynomials = Array{Array{T,2},1}(undef,length(unique_orders))

    # Construct the base polynomials

    for i = 1:length(unique_orders)
      base_polynomials[i] = chebyshev_polynomial(unique_orders[i],nodes[k,:])
    end

    # Compute the unique polynomial terms from the base polynomials

    for i = length(unique_orders):-1:2
      unique_base_polynomials[i] = base_polynomials[i][:,size(base_polynomials[i-1],2)+1:end]
    end
    unique_base_polynomials[1] = base_polynomials[1]

    # Construct the first row of the interplation matrix

    l = 1
    @inbounds for j = 1:size(multi_index,1)
      new_polynomials = unique_base_polynomials[multi_index[j,1]][1,:]
      for i = 2:size(multi_index,2)
        new_polynomials = kron(new_polynomials,unique_base_polynomials[multi_index[j,i]][i,:])
      end
      m = length(new_polynomials)
      interpolation_matrix[k,l:l+m-1] = new_polynomials
      l += m
    end

  end

  inverse_interpolation_matrix = inv(interpolation_matrix)

  return inverse_interpolation_matrix

end

function smolyak_inverse_interpolation_matrix_threaded(nodes::Union{Array{T,1},Array{T,2}},multi_index::Array{S,2},domain::Union{Array{T,1},Array{T,2}}) where {T<:AbstractFloat,S<:Integer}

  # Normalize nodes to the [-1.0 1.0] interval

  d = size(multi_index,2)
  nodes = copy(nodes)
  for i = 1:d
    nodes[:,i] = normalize_node(nodes[:,i],domain[:,i])
  end

  inverse_interpolation_matrix = smolyak_inverse_interpolation_matrix_threaded(nodes,multi_index)

  return inverse_interpolation_matrix

end

function smolyak_weights(y::AbstractArray{T,1},inverse_interpolation_matrix::Array{T,2}) where {T<:AbstractFloat}

  weights = inverse_interpolation_matrix*y

  return weights

end


text:
# Autogenerated wrapper script for primecount_jll for i686-linux-musl-cxx03
export libprimecount, libprimesieve, primecount

using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("primecount")
JLLWrappers.@declare_library_product(libprimecount, "libprimecount.so.6")
JLLWrappers.@declare_library_product(libprimesieve, "libprimesieve.so.9")
JLLWrappers.@declare_executable_product(primecount)
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libprimecount,
        "lib/libprimecount.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libprimesieve,
        "lib/libprimesieve.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        primecount,
        "bin/primecount",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
#!/usr/bin/env julia

# Given some model parameters, plot the channel maps and
# integrated spectrum

using ArgParse

s = ArgParseSettings()
@add_arg_table s begin
    "--config"
    help = "a YAML configuration file"
    default = "config.yaml"
end

parsed_args = parse_args(ARGS, s)

import YAML
config = YAML.load(open(parsed_args["config"]))

using DiskJockey.constants
using DiskJockey.image
using DiskJockey.model
using HDF5

import PyPlot.plt
using LaTeXStrings
import Images

species = config["species"]
transition = config["transition"]
lam0 = lam0s[species*transition]
model = config["model"]

# Choose a diverging colorscheme for in front of/ behind the plane
cmap = plt[:get_cmap]("coolwarm")

# Use the same configuration as plot_chmaps, but we'll use different scaling.
function plot_tausurf(img::image.TausurfImage; fname="tausurf.png")

    data = img.data ./ AU # Convert from cm to AU

    vvmax = maxabs(data) # [AU]

    norm = PyPlot.matplotlib[:colors][:Normalize](-vvmax, vvmax)

    (im_ny, im_nx) = size(img.data)[1:2] # y and x dimensions of the image

    xxx = ((Float64[i for i=0:im_nx-1] + 0.5) - im_nx/2.) * img.pixsize_x ./ AU
    yyy = ((Float64[i for i=0:im_ny-1] + 0.5) - im_ny/2.) * img.pixsize_y ./ AU

    ext = (xxx[1], xxx[end], yyy[1], yyy[end])

    # Figure out how many plots we'll have.
    ncols = 8
    nrows = ceil(Int, nlam/ncols)

    xx = 1.5 * 9
    dx = 1.5
    yy = (nrows + 1) * 1.5
    dy = 1.5

    fig, ax = plt[:subplots](nrows=nrows, ncols=ncols, figsize=(xx, yy))

    for row=1:nrows
        for col=1:ncols
            iframe = col + (row - 1) * ncols

            if col != 1 || row != nrows
                ax[row, col][:xaxis][:set_ticklabels]([])
                ax[row, col][:yaxis][:set_ticklabels]([])
            else
                ax[row, col][:set_xlabel](L"$xx$ (AU)")
                ax[row, col][:set_ylabel](L"$yy$ (AU)")
            end

            if iframe > nlam
                # Plot a blank square if we run out of channels
                ax[row, col][:imshow](zeros((im_ny, im_nx)), cmap=cmap, vmin=0, vmax=20, extent=ext, origin="lower")

            else
                frame = data[:,:,iframe]

                im = ax[row, col][:imshow](frame, extent=ext, interpolation="none", origin="lower", cmap=cmap, norm=norm)

                if iframe==1
                    # Plot the colorbar
                    cax = fig[:add_axes]([(xx - 0.35 * dx)/xx, (yy - 1.5 * dy)/yy, (0.1 * dx)/xx, dy/yy])
                    cbar = fig[:colorbar](mappable=im, cax=cax)

                    cbar[:ax][:tick_params](labelsize=6)
                    fig[:text](0.99, (yy - 1.7 * dy)/yy, "AU", size=8, ha="right")
                end

                ax[row, col][:annotate](@sprintf("%.1f", vels[iframe]), (0.1, 0.8), xycoords="axes fraction", size=8)
            end

        end
    end

    fig[:subplots_adjust](hspace=0.00, wspace=0.00, top=(yy - 0.5 * dy)/yy, bottom=(0.5 * dy)/yy, left=(0.5 * dx)/xx, right=(xx - 0.5 * dy)/xx)

    plt[:savefig](fname, dpi=600)
end

# Read the image
im = tauread()
global nlam = length(im.lams)
global vels = c_kms * (im.lams .- lam0)/lam0

plot_tausurf(im)

max_stars_repo_name: elnjensen/DiskJockey

max_issues_repo_name: elnjensen/DiskJockey

max_forks_repo_name: elnjensen/DiskJockey


text:
@testset "Markov1" begin
    m = Markov()
    @tr m begin
        :up => :down, 1.0
        :down => :up, 100.0
    end
    @init m begin
        :up, 1.0
    end
    @reward m :avail begin
        :up, 1.0
    end
    Q, initv, rwd, = generate(m)
    println(initv)
    println(Q)
    println(rwd)
end

@testset "Markov1_1" begin
    m = Markov()
    @tr m begin
        :up => :down, 1.0
        :down => :up, 100.0
    end
    @init m begin
        :up, 1.0
    end
    @reward m :avail begin
        :up, 1.0
    end
    Q, initv, rwd, = generate(m, modeltype=:DenseCTMC)
    println(initv)
    println(Q)
    println(rwd)
end

@testset "Markov2" begin
    @vars lam1 lam2
    m = Markov()
    @tr m begin
        :up => :down, lam1
        :down => :up, lam2
    end
    @init m begin
        :up, 1.0
    end
    @reward m :avail begin
        :up, 1
    end
    Q, initv, rwd, = generate(m)
    println(initv)
    println(Q)
    println(rwd)

    avail = dot(sprob(Q), rwd[:avail])

    @bind begin
        lam1 = 1.0
        lam2 = 100.0
    end

    a = seval(avail)
    println(a)
    da1 = seval(avail, :lam1)
    println(da1)
    da2 = seval(avail, :lam2)
    println(da2)
    da12 = seval(avail, (:lam1,:lam2))
    println(da12)
end

@testset "Markov3" begin
    @markov reliab(lam1, lam2) begin
        @tr begin
            :up => :down, lam1
            :down => :up, lam2
        end
        @init begin
            :up, 1.0
        end
        @reward :avail begin
            :up, 1
        end
    end

    @bind begin
        lam1 = 1.0
        lam2 = 100.0
    end
    m = reliab(lam1, lam2)
    println(ctmc)

    model = ctmc(m)
    avail = exrss(model, reward=:avail)

    a = seval(avail)
    println(a)
    da1 = seval(avail, :lam1)
    println(da1)
    da2 = seval(avail, :lam2)
    println(da2)
    da12 = seval(avail, (:lam1,:lam2))
    println(da12)
end

@testset "Markov4" begin
    @markov reliab(lam1, lam2) begin
        @tr :up => :down, lam1
        @tr :down => :up, lam2
        @init :up, 1.0
        @reward :avail :up, 1
    end
    m = reliab(1, 100)
    println(m)

    avail = exrss(ctmc(m), reward=:avail)
    println(avail)
end

@testset "Markov5" begin
    @markov reliab(lam1, lam2) begin
        @tr :up => :down, lam1
        @tr :down => :up, lam2
        @init :up, 1.0
        @reward :avail :up, 1
    end
    m = reliab(1, 100)
    println(m)

    avail = exrss(ctmc(m), reward=:avail)
    println(avail)

    tavail = exrt(LinRange(0, 1, 10), ctmc(m), reward=:avail)
    println(tavail)
end

@testset "Markov6" begin
    @markov reliab(lam1, lam2) begin
        @tr :up => :down, lam1
        @tr :down => :up, lam2
        @init :up, 1.0
        @reward :avail :up, 1
    end
    m = reliab(1, 100)
    println(m)

    avail = exrss(ctmc(m), reward=:avail)
    println(avail)

    tavail = cexrt(LinRange(0, 1, 10), ctmc(m), reward=:avail)
    println(tavail)
end

@testset "Markov7" begin
    @markov reliab(lam1, lam2) begin
        @tr :up => :down, lam1
        @tr :down => :up, lam2
        @init :up, 1.0
        @reward :avail :up, 1
    end
    m = reliab(1, 100)
    println(m)

    avail = exrss(ctmc(m), reward=:avail)
    println(avail)

    tavail = cexrt(0.5, ctmc(m), reward=:avail)
    println(tavail)
end

@testset "Markov8" begin
    @markov reliab(lam1, lam2) begin
        @tr :up => :down, lam1
        @tr :down => :up, lam2
        @init :up, 1.0
        @reward :avail :up, 1
    end

    @bind begin
        lam1 = 1.0
        lam2 = 100.0
    end

    m = reliab(lam1, lam2)
    println(m)

    avail = exrss(ctmc(m), reward=:avail)
    println(avail)

    tavail = exrt(1.0, ctmc(m), reward=:avail)
    println(tavail)

    a = seval(tavail)
    println(a)
    da1 = seval(tavail, :lam1)
    println(da1)
    da2 = seval(tavail, :lam2)
    println(da2)
    da12 = seval(tavail, (:lam1,:lam2))
    println(da12)
end

@testset "Markov9" begin
    e = @macroexpand @markov mm1k(lam,mu,k) begin
        for i = 0:k-1
            @tr Symbol(i) => Symbol(i+1), lam
        end
        for i = 1:k
            @tr Symbol(i) => Symbol(i-1), mu
        end
        @init @s(0), 1.0
        for i = 0:k
            @reward :len Symbol(i), i
        end
    end
    println(e)
end

@testset "hybrid2" begin
    @markov midplane(lam, mu) begin
        @tr begin
            :up => :down, lam
            :down => :up, mu
        end
        @reward :r begin
            :up, 1
        end
    end

    @bind begin
        lam = 0.4
        mu = 1.5
        delta = 1.5
    end

    m = ctmc(midplane(lam, mu))
    println(seval(sprob(m.Q), lam))
    println(seval([lam, mu], delta))
    println(seval(sprob(m.Q), delta))
    println(seval(sprob(m.Q), (lam, delta)))
end

@testset "m1" begin
    @markov markov1(lamA, lamB) begin
        @tr :AB  => :AX, lamB
        @tr :AB => :XB, lamA
        @tr :AX => :XX1, lamA
        @tr :XB => :XX2, lamB
        @init :AB, 1.0
        @reward :failureA begin ### define the reward vector labeled by :faulreA
            :XB,  1
            :XX1, 1
            :XX2, 1
        end
        @reward :failureAdash begin  ## define the reward vector labeled by :faulreAdash
            :XX2,  1
        end
    end
    @bind begin
        lamA = 0.1
        lamB = 0.5
    end
    m1 = ctmc(markov1(lamA, lamB))
    ts = 10.0
    fa = exrt(ts, m1, reward=:failureA)
    display(fa)
end


text:
mutable struct NfLattice{T}
  basis::Vector{T}
  discriminant::fmpz
  is_minkowski_exact::Bool
  minkowski_gram_exact::fmpz_mat
  minkowski_gram_scaled::Tuple{Int, fmpz_mat}
  function NfLattice{T}(v::Vector{T}, discriminant::fmpz) where {T <: NumFieldElem}
    n = new{T}()
    n.basis = v
    n.discriminant = discriminant
    is_minkowski_exact = false
    return n
  end
end

function lattice(v::Vector{T}, disc::fmpz; isexact::Bool = false) where T <: NumFieldElem
  @assert !isempty(v)
  L = NfLattice{T}(v, abs(disc))
  L.is_minkowski_exact = isexact
  return L
end

function dim(L::NfLattice)
  return length(L.basis)
end

function basis(L::NfLattice)
  return L.basis
end

function discriminant(L::NfLattice)
  return L.discriminant
end

function nf(L::NfLattice{T}) where T <: NumFieldElem
  return parent(basis(L)[1])
end

function minkowski_matrix(L::NfLattice, p::Int)
  return minkowski_matrix(basis(L), p)
end

function minkowski_gram_mat_scaled(L::NfLattice, p::Int)
  if L.is_minkowski_exact
    L.minkowski_gram_exact = _exact_minkowski_matrix(basis(L))
    return L.minkowski_gram_exact
  end
  K = nf(L)
  if isdefined(L, :minkowski_gram_scaled) && L.minkowski_gram_scaled[1] >= p
    A = deepcopy(O.minkowski_gram_mat_scaled[1])
    shift!(A, p - O.minkowski_gram_mat_scaled[2])
  else
    c = minkowski_matrix(L, p)
    B = basis(L)
    d = zero_matrix(FlintZZ, length(B), absolute_degree(K))
    A = zero_matrix(FlintZZ, length(B), length(B))
    round_scale!(d, c, p)
    ccall((:fmpz_mat_gram, libflint), Nothing, (Ref{fmpz_mat}, Ref{fmpz_mat}), A, d)
    shift!(A, -p)
    L.minkowski_gram_scaled = (p, deepcopy(A))
  end
  for i=1:absolute_degree(K)
    fmpz_mat_entry_add_ui!(A, i, i, UInt(nrows(A)))
  end
  return A
end

function weighted_minkowski_gram_scaled(L::NfLattice, v::fmpz_mat, prec::Int)
  c = deepcopy(minkowski_matrix(L, prec))
  mult_by_2pow_diag!(c, v)
  d = zero_matrix(FlintZZ, nrows(c), ncols(c))
  round_scale!(d, c, prec)
  g = zero_matrix(FlintZZ, nrows(c), nrows(c))
  ccall((:fmpz_mat_gram, libflint), Nothing, (Ref{fmpz_mat}, Ref{fmpz_mat}), g, d)
  shift!(g, -prec)
  for i=1:n
    fmpz_mat_entry_add_ui!(g, i, i, UInt(nrows(d)))
  end
  return g
end

function lll(L::NfLattice, weights::fmpz_mat = zero_matrix(FlintZZ, 1, 1); starting_prec::Int = 100)
  if L.is_minkowski_exact
    M = _exact_minkowski_matrix(basis(L))
    l, v = lll_gram_with_transform(M)
    return FakeFmpqMat(l), v
  end
  prec = starting_prec
  local l1::FakeFmpqMat
  local v::fmpz_mat
  while true
    if prec > 2^18
      error("Something wrong in short_elem")
    end
    try
      l1, v = _lll(L, weights, prec)
      break
    catch e
      if !(e isa LowPrecisionLLL || e isa InexactError)
        rethrow(e)
      end
      prec *= 2
    end
  end
  return l1, v
end

function _lll(L::NfLattice, weights::fmpz_mat, prec::Int)
  if iszero(weights)
    d = minkowski_gram_mat_scaled(L, prec)
    sv = fmpz(0)
  else
    d = weighted_minkowski_gram_scaled(L, weights, prec)
    sv = max(fmpz(0), sum(weights[1, i] for i=1:ncols(weights)))
  end
  n = dim(L)
  g = identity_matrix(FlintZZ, n)
  ctx = Nemo.lll_ctx(0.99, 0.51, :gram)
  ccall((:fmpz_lll, libflint), Nothing, (Ref{fmpz_mat}, Ref{fmpz_mat}, Ref{Nemo.lll_ctx}), d, g, ctx)
  l, t = d, g
  if nbits(maximum(abs, t)) >  div(prec, 2)
    throw(LowPrecisionLLL())
  end
  n = nrows(l)
  disc = discriminant(L) * fmpz(2)^(2*sv)
  di = root(disc, n)+1
  di *= fmpz(2)^(div(n+1,2)) * fmpz(2)^prec
  if compare_index(l, 1, 1, di) > 0
    throw(LowPrecisionLLL())
  end
  pr = prod_diagonal(l)
  if pr > fmpz(2)^(div(n*(n-1), 2)) * disc * fmpz(2)^(n*prec)
    throw(LowPrecisionLLL())
  end
  return FakeFmpqMat(l, fmpz(2)^prec), t
end

function lll_basis(L::NfLattice{T}) where T
  K = nf(L)
  l, t = lll(L)
  B = basis(L)
  new_basis = Vector{T}(undef, dim(L))
  tmp = K()
  for i = 1:length(new_basis)
    a = K()
    for j = 1:ncols(l)
      tmp = B[j]*t[i, j]
      a = add!(a, a, tmp)
    end
    new_basis[i] = a
  end
  return new_basis
end

_abs_disc(O::NfRelOrd) = absolute_discriminant(O)
_abs_disc(I::NfRelOrdIdl) = absolute_discriminant(order(I))*absolute_norm(I)

function lll_basis(OL::T) where T <: Union{NfRelOrdIdl, NfRelOrd}
  L = nf(OL)
  B = pseudo_basis(OL, copy = false)
  ideals = Dict{typeof(B[1][2]), Vector{elem_type(base_field(L))}}()
  for i = 1:length(B)
    if haskey(ideals, B[i][2])
      continue
    end
    ideals[B[i][2]] = lll_basis(B[i][2])
  end
  abs_bas = Vector{elem_type(L)}(undef, absolute_degree(L))
  ind = 1
  for i = 1:length(B)
    J = ideals[B[i][2]]
    for j = 1:length(J)
      abs_bas[ind] = J[j]*B[i][1]
      ind += 1
    end
  end
  isexact = false
  if istotally_real(L)
    isexact = true
  end
  V = lattice(abs_bas, _abs_disc(OL), isexact = isexact)
  return lll_basis(V)
end


text:
using ExprParsers
using Documenter

makedocs(;
    modules=[ExprParsers],
    authors="Stephan Sahm <stephan.sahm@gmx.de> and contributors",
    repo="https://github.com/schlichtanders/ExprParsers.jl/blob/{commit}{path}#L{line}",
    sitename="ExprParsers.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://schlichtanders.github.io/ExprParsers.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Manual" => "manual.md",
        "Library" => "library.md",
    ],
)

deploydocs(;
    repo="github.com/schlichtanders/ExprParsers.jl",
)


text:
"""
Simulator object. Holds all relevant variables and interfaces with AdaStress.
"""
Base.@kwdef mutable struct Simulator <: AdaStress.GrayBox
	n::Z                        = 2
    step::Z                     = 0
    num_steps::Z                = 50
	dt::R                       = 1.0
	seed::Z                     = 0
    randomize::Bool             = true
	enc::Encounter              = Encounter(n=n, randomize=randomize, seed=seed)
	acs::Vector{Aircraft}       = @array n Aircraft()
	cass::Vector{<:CAS}         = @array n RuleBasedCAS()
	prs::Vector{PilotResponse}  = @array n PilotResponse()
	trs::Vector{<:Distribution} = @array n MvNormal([10.0, 3.0, 3.0])
	nmac::Tuple{R,R}            = (50.0, 15.0)
    reduce_symmetries::Bool     = true
    metrics::Metrics            = Metrics()
    logging::Bool               = false
    log::Dict{String,Any}       = Dict{String, Any}()
end

function AdaStress.reset!(sim::Simulator)
	sim.step = 0
	initialize(sim.enc)
	initialize.(sim.acs, sim.enc.acs)
	initialize.(sim.cass)
	initialize.(sim.prs)
    initialize(sim.metrics, sim.acs, sim.nmac)
    initialize(sim.log, sim)
end

function AdaStress.environment(sim::Simulator)
    AdaStress.Environment(Symbol("cmd_$i") => sim.trs[i] for i in 1:sim.n)
end

function AdaStress.observe(sim::Simulator)
    if sim.n == 2 && sim.reduce_symmetries
        kernel(sim.acs...)
    else
        reduce(vcat, observation.(sim.acs))
    end
end

function AdaStress.step!(sim::Simulator, x::AdaStress.EnvironmentValue)
	sim.step += 1
	cmds = [Command(x[Symbol("cmd_$i")]) for i in 1:sim.n]
    update.(sim.cass, Ref(sim.acs), 1:sim.n)
    update.(sim.prs, cmds, sim.cass)
	update.(sim.acs, (pr -> pr.cmd).(sim.prs), sim.dt)
    update(sim.metrics, sim.acs, sim.nmac)
	update(sim.log, sim)
end

AdaStress.isterminal(sim::Simulator) = sim.step >= sim.num_steps

AdaStress.isevent(sim::Simulator) = sim.metrics.d["is_nmac"]

AdaStress.distance(sim::Simulator) = sim.metrics.d["d"]

AdaStress.flatten(dist::MvNormal, val::Vector{<:Real}) = sqrt(invcov(dist)) * val

AdaStress.unflatten(dist::MvNormal, val::Vector{<:Real}) = sqrt(cov(dist)) * val


text:
"""
Function prints numbers with related to them colors in `printstyled`
"""
function palette(bold::Bool=true)

    print("\n    Choose your color:\n\n")
    for i = 1:255

        if ((i-16)%12 == 0)&(i>15)

            print("\n")

        end

        printstyled(string(i) * repeat(" ",6 - length(string(i)) - 1*(i<16) );color = i, bold = bold )

        if (i == 15)

            print("\n")

        end

        if ((i-16)%6 == 5)&(i>15)

            print("   ")

        end

    end

end

"""
Use it when you really need to be noticed:

# Examples
```julia-repl
julia> printclowned("I want to be noticed!!!\\n\\nRight now!!!)")
```
"""
function printclowned(text::String)

    for t ∈ text

        printstyled(t, color = rand(1:255), bold = rand( (true,false)))

    end

end


text:
using Statistics
using CairoMakie, ColorSchemes

import Base: parse, +, /

struct Params
    N :: Int
    L :: Int
end

function parse(::Type{Params}, s::AbstractString)
    element = split(s,';')
    return Params(
        parse(Int,element[1]),
        parse(Int,element[2]),
    )
end

struct Duration
    hour   :: Int
    minute :: Int
    second :: Float64
end

function parse(::Type{Duration}, s::AbstractString)
    index = (
        o = 0,
        h = findfirst('h', s),
        m = findfirst('m', s),
        s = findfirst('s', s),
    )

    inc(::Nothing) = 1
    inc(x::Int64)  = x+1

    h = index.h === nothing ? 0 : parse(Int, s[inc(index.o):index.h-1])
    m = index.m === nothing ? 0 : parse(Int, s[inc(index.h):index.m-1])
    s = index.s === nothing ? 0 : parse(Float64, s[inc(index.m):index.s-1])

    return Duration(h,m,s)
end

seconds(t::Duration) = t.second + 60*t.minute + 3600*t.hour

function +(t1::Duration, t2::Duration)
    s = t1.second + t2.second
    m = t1.second + t2.second + (s÷60)
    h = t1.hour + t2.hour + (m÷60)

    return Duration(h, m%60, s%60)
end

function /(t::Duration, x::Int)
    h = t.hour ÷ x
    m = 60*(t.hour % x) + t.minute ÷ x
    s = 60*(t.minute % x) + t.second ÷ x

    m += s % 60
    s %= 60

    h += m % 60
    m %= 60

    return Duration(h, m, s)
end

function benchmarks(io)
    duration = Dict{Params,Array{Duration,1}}()
    for line in eachline(io)
        key, val = strip.(split(line, "=>"))
        param = parse(Params, key)
        dtime = parse(Duration, val)

        if param in keys(duration)
            push!(duration[param], dtime)
        else
            duration[param] = Duration[dtime]
        end
    end

    return duration
end

field(data, key) = Set(getfield.(keys(data), key)) |> collect |> sort
function plots(benchmark)
    L = field(benchmark, :L)
    N = field(benchmark, :N)

    colors = cgrad(:matter, length(L), categorical=true)

    fig  = Figure(font="Latin Modern Math", fontsize=26)
    axis = Axis(fig[1,1],
                xscale=CairoMakie.log10,
                yscale=CairoMakie.log10,
                xlabel="number of genomes",
                ylabel="time(sec)",
                xticks=([10, 100, 1000],    [L"10^1", L"10^2", L"10^3"]),
                yticks=([1, 10, 100, 1000], [L"10^0", L"10^1", L"10^2", L"10^3"]),
    )
    hidespines!(axis, :t, :r)
    for (i,l) in enumerate(L)
        data = [seconds.(benchmark[Params(n,l)]) for n in N]
        μ, σ = mean.(data), std.(data)
        lines!(axis, N, μ;
            color=colors[i],
            linewidth=2,
            label="$(l)"
        )

        band!(axis, N, μ-σ, μ+σ;
            color=(colors[i],0.25),
        )
    end
    lines!(N, .5*N, color=:black, linewidth=3, linestyle=:dash, label="linear trend")
    axislegend("Avg. genome length", position = :lt, nbanks=2)
    fig
end

function main(path, dest)
    plt = open(benchmarks, path) |> plots
    save(dest, plt, px_per_unit=2)
end

if abspath(PROGRAM_FILE) == @__FILE__
    main(ARGS[1], ARGS[2])
end

max_stars_repo_name: neherlab/pangraph

max_issues_repo_name: neherlab/pangraph

max_forks_repo_name: neherlab/pangraph